RNAseq定量方法

生信喵

一、比对参考基因组还是参考基因集
       为了获取表达矩阵，可以将测序数据比对到参考基因组然后通过坐标文件 GTF（GFF 或者 BED）统计每个基因比对上的数据计算丰度，或者直接与参考基因集进行比对，直接计算每个基因覆盖深度的方法。但是两种方法之间有较大的差别：
       第一：基因集中不能保证所有基因都是准确的，例如一个基因预测结果不准确，比真实情况过长，那么长度来的部分则没有序列覆盖到；
       第二：基因集都是已知的基因，一些没有预测到的基因无法比对上，丢失了这部分信息，利用全基因组比对，可以对已知基因进行纠正；
       第三：基因集中不包括可变剪切的信息，而与全基因组序列比对，然后综合GTF 文件中已知转录本信息，则可以找到一些新的转录本；
       第四：全基因组有更长的序列，因此 reads 比对过程中更容易比对上，基因长度过短，头尾部不容易比对上，这样就会对最终计数产生影响。

二、基因水平差异还是转录本水平差异？
       真核生物由于可变剪切的存在，同一个基因通过外显子的不同组合，会表达出不同的转录本，最终翻译出不同的蛋白质序列。那么对于 RNAseq 研究，是在基因水平还是转录本水平比较差异更准确呢？
       严格来说，应该比较转录本水平的表达差异更加准确。但由于目前测序读长过短，假设有来自同一基因的三个转录本 A，B，C，长度分别是 1000bp，800bp，600bp，假设测序读长 2*100bp，最终测序出来的 reads 无法准确分配给最终的转录本 A，B 或者 C，并且还可能存在新的转录本 D，因此，当前条件下采用基因水平的差异表达计算。
       目前开发出一些工具可以用户转录本水平的定量，例如 RSEM 工具。

三、几种 RNAseq 定量方法比较
       给定两个基因，如何比较是否存在表达差异呢？如何进行量化。由于来自不同样品，并非测序全长，测序深度不同，可变剪切等因素的影响，不能直接比较覆盖 reads 数，因此需要进行均一化。由于我们通常对同一个基因进行差异表达的比较，因此，均一化的主要作用是消除测序深度的影响。目前采用的主要方法是 RPKM，FPKM，TPM，CPM 等方法，每一种方法都不是完美的，都有其特定的使用范围。
       无论采用哪种比较方法，我们都应该知道下面两个原则：
       1、RNAseq 定量都是一种相对定量的方法，不同的定量方法，结果会有所差异；
       2、一般研究只关注差异最大的一些基因，因此，无论采用哪种定量方法，差别最显著的都会凸显出来。

1. https://www.reneshbedre.com/blog/expression_units.html

复制代码

3.1 reads count
       测序完成之后，每一个基因被测序到的 reads 数目，理论上来说，基因表达量越大，基因长度越长，测序深度越大，被测序到的 reads 数目越多。

       不同基因比对得到的 Read 个数
      
       注意事项：由于可变剪切的存在，计算 reads counts 可以相对于基因，也可以相对于转录本，如果是相对于转录本就比较复杂，存在多重比对的情况。


3.2 RPKM
       RPKM 全称是（Reads Per Kb per Million reads）（Mortazavi, 2008），
       其计算公式为：

       nr 是比对至目标基因的 read 数量；
       L 是目标基因的外显子长度之和除以 1000
       N 是全部 reads 数目

       不同基因的 RPKM 值的计算结果（“M”为 100 时）
      

3.3 FPKM
       FPKM：FPKM（Fragments PerKilobase Million），FPKM 和 RPKM 的计算方法基本一致，只不过把 reads 换成了 Fragments。

      
       nf 是比对至目标基因的 fregment 数量
       对于单末端测序数据，由于 Cufflinks 软件计算的时候是将一个 read 当做一个fragment 来算的，故而 FPKM 等同于 RPKM。对于双末端测序，如果一对paired-read 都比对上了，那么这一对 paired-read 称之为一个 fragment；而如果只有一个比对上了，另外一个没有比对上，那么就将这个比对上的 read 称为一个 fragment。而计算 RPKM 时，一对 paired-read 会当成两个 read 分别计算。


3.4 TPM
       TPM（Transcripts PerKilobase Million）计算公式：

      
       下面是依照公式计算上表中基因的 TPM 值（注意：为了便于显示，同上，将公式中的 100 万改为 100）：

       不同基因的 TPM 值的计算结果（“M”为 100 时）
      

3.5 CPM
       CPM（counts per million),在 edgeR 中，提供了一种名为 CPM 的定量方式，全称为 count-per-millon。
       假定原始的表达量矩阵为 count, 计算 CPM 的代码如下
       cpm <- apply(count ,2, function(x) { x/sum(x)*1000000 })
       原始的表达量除以该样本表达量的总和，在乘以一百万就得到了 CPM 值 。从公式可以看出， CPM 其实就是相对丰度，只不过考虑到测序的 reads 总量很多，所以总的 reads 数目以百万为单位。