单细胞normalization

努力上进的三心草

已于 2022-08-05 09:42:31 修改

阅读量937

点赞数 1

文章标签： sql

于 2022-07-10 16:33:44 首次发布

本文链接：https://blog.csdn.net/ruiyan_hou/article/details/125707453

版权

基本概念

RPKM （reads per kilobase million)
FPKM (fragments per kilobase million)
TPM (transcripts per million)
CPM (counts per million)
对于这几个概念的阐述，我感觉https://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/ 这篇博客里面阐述的比较清楚。
我觉得可以按照时间的发展顺序来考虑前三个概念：RPKM和FPKM首先被提出，接着TPM对前面的概念做了改进。
FPKM和RPKM可以放在一起进行，只不过由于pair-end sequence 有两条reads,它就会比对到基因的两个部分，这样我们不能再以比对上的reads数作为评价标准了，就以能够比对上的fragments数目作为评价标准。所以说，对于单端测序用RPKM，对于双端测序用FPKM，接下来我们就以展示RPKM为主。

对于RPKM，我们来想一个问题，我们平时说，一个基因，它的表达高，或者它的表达低，是不是看有多少条reads能够比对到这个基因上。在这个过程中，我们想想，有什么因素在影响比对的这个reads数目。
1、基因的长度。我们最终可能得出的结论是基因A在某一个样本中的表达要高于基因B，那么就是直接把基因A和基因B放在一起进行比较的。那么这种情况下，必然要保证消除它们的长度对这个基因的影响。
2、测序的深度。我们要保证测定的总的reads数目在需要比较的几组之间应该是相同的，如果不能保证怎么办，那就除以这个总的reads的数目，消除它的影响。
基于以上的想法，人们提出了RPKM。但是要注意，虽然RPKM的计算公式是以上那样的，但是具体的计算过程中，是通过以下步骤计算的
1）把一个sample里面所有的reads数加起来,然后除以1 million，最后得到的这个数就称为scaling factor (这也就是名称里面的per million （但是在这里，这个per million可能就是一个合适的数，在上面的公式里面，它就用的是10^9)
2)然后用上面得到的结果/transcript length，最后得到的结果就称为 per kilobase ( 因为基因的长度的单位一般都是kb)

解释完RPKM之后，我们就来看一看TPM。首先明确TPM是对RPKM的改进，因此它有它的优势。

TPM的计算步骤和RPKM相反
１）首先用能够比对到每一个transcript上的reads的数目/transcript的长度。
2）把上一步除过长度之后的reads数目加起来,然后除以10^6,这个得到的数就是scaling factor.
3）用第一步的结果/第二步的结果，得到的结果就是per million 。

然后：通过在相除的时候抵消一些东西，TPM可以由RPKM计算得到。

最后我们来看一看CPM吧。CPM不涉及transcript　length
CPM= number of reads mapped to gene * 10^6 /Total number of mapped reads
它只考虑了测序深度的影响。

实际应用

应用一（10x和smart-seq2的normalization）

在单细胞测序，常用的单细胞测序方法主要有两种，一种是smart-seq2测序，另外一种是10x genomics测序。前者是测定基因的全长，后者是测定3‘端，或者5’端，仅仅测定一部分固定的长度。根据这个测序原理的话，对于smart-seq2测序技术，就应该以TPM作为单位进行标准化，对于10x genomic就应该以CPM为单位进行标准化。
因为对于10x genomic的数据，我们常用的scanpy来说，它有一步normalization
sc.pp.normalize_total(adata, target_sum=1e4)
https://scanpy.readthedocs.io/en/stable/generated/scanpy.pp.normalize_total.html
可以仔细参考它的说明文档，一般我们用的时候就根据它的tutorial把target_sum设置成为1e4, 当把target_sum设置成为1e6的时候，最后得到的值就是CPM。（这个1e4和1e6的关系就像上面说的那个10^6和109的关系一样。
但是对于smart-seq2数据就不能直接用上面的normalization的方法了。一个经验教训就是在使用scVelo的时候直接套用10x的方法，最终出错。对于使用scVelo，如果是smart-seq2数据的话，就应该把X，layers里面的spliced, unspliced都进行TPM标准化才可以。

应用二（pseudo bulk)

参考：https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-022-01021-1#Sec2
在这篇文章中，作者使用hierarchical cluster来阐明不同COVID-19严重状况下，不同cell type的transcription profile的相似性。详细来说就是，现在又四个COVID-19的状况：正常，轻度，中度和重度。对于每一种情况都会有不同的细胞type,现在作者想看看是不同的COVID-19的状况占主导，还是不同的cell type占主导，如果是前者占主导的话，那么同一个COVID-19的情况下的不同的cell type会聚在一起。如果是后者占主导的话，那么同一个cell type会聚在一起，即使是不同的COVID-19的情况。最后展示出来的是这样一张图

假设总共有四种状况，A B C D，每一种状况下面又有5个cell type，那么总共就有20行，20列。每一个格子里面的数值是每一个状况下的细胞类型和其他的细胞类型的相关性。也叫做Pearson correlation coefficient, 也就是PPC。同理，这个PPC其实也可以是JSD （Jensen-Shannon divergence，也就是不同的cell type之间的差异）。
那么如何计算出来这个correlation coefficient呢？这时候就只能使用pseudo bulk了，也就是把这个cluster 下面的所有细胞当成是一个sample。然后用能够比对到这个gene上的reads count/ total counts , 然后再x 10^6,只有这样进行normalization之后才能去除测序深度的影响，使不同的cell type之间具有可比性。举了下面的例子。

另外一篇文献中展示的pseudo bulk的方法。先把同一个cell type 中的所有的细胞加起来，然后再用TMM进行normalized
在这里插入图片描述