github-actions[bot]
commented
4 years ago FPKM与TPM的故事 by 生信小知识
微信公众号：生信小知识
关注可了解更多的教程及单细胞知识。问题或建议，请公众号留言;
FPKM与TPM的故事
前言FPKM和RPKM都错了？为什么FPKM/RPKM是错的？举例理解为什么TPM的对的？举例理解为什么还在用FPKM小结RNA-seq的protoccol差异不同组织的RNA-seq数据不能直接比较mRNA的空间定位RNA-seq的链特异性不同mRNA整体水平的细胞不可比较小结
前言
数据的normalize非常的重要，我想所有懂一点生信的人都明白这个道理的，但是，似乎现在没有人好好总结下关于数据normalization的问题，或者说，没有人能让不是生信专业的人弄明白到底normalizetion是有什么问题。所以，我决定整理下关于normalization的问题。
好吧，其实主要是因为我组会要讲，哈哈哈哈哈哈哈哈！！！
这里先解读一篇关于RPKM和TPM的文献吧~
FPKM和RPKM都错了？
这里其实有很多人都讲过什么是FPKM/RPKM，因为已经有很多很不错的说明了，我就不多说了，具体可以去看非常有深度，但是又非常容易理解的推文：
这个文章提出的观点让我觉得有点有趣，于是就pubmed了一下，结果就又找到了一篇文章：
http://www.rnajournal.org/cgi/doi/10.1261/rna
经过阅读这篇推文和文献，我总结下我的体会。
为什么FPKM/RPKM是错的？
我们可以假定，对于样本X，其有一个基因g被转录了mRNA_g次，同时样本X中所有基因的转录总次数假定是mRNA_total，那么正确描述基因g转录丰度的值应该是：
同时，假设这个样本X中一共有n种基因，那么其他基因的转录丰度的计算也应该和上述公式类似：
除了要把分子mRNA_g换为其他基因对应的转录次数之外，分母mRNAtotal都一样。于是有趣的事情就是，所有基因转录本丰度的均值rmean将是一个恒定不变的数，由以上定义这个数就是：
这个期望值和测序状态无关！仅仅由样本中基因的总数n所决定的。也就是说，对于同一个物种，不管它的样本是哪种组织（正常的或病变的），也不管有多少个不同的样本，只要它们都拥有相同数量的基因，那么它们的r_mean都将是一致的。这是一个在进行比较分析的时候，非常有意义的恒等关系。
但在实际的操作中，我们是无法直接计算这些转录丰度r值的。好在只要能够保证模型本身的自洽性，我们是能通过自建一些统计量来对r值进行间接描述的。比如FPKM和RPKM，实际上它们的目的就是为了描述r。
但这里的大前提就是，所有这些要用来描述转录本丰度的统计量，都应该能等价表示出这一恒等关系。即：
不管如何，这些统计量所描述出来的转录本丰度应该得是其真实丰度rg的m倍（m必须是一个根据模型定出的不变值），均值也将是rmean的m倍，只有这样才是得出有意义的结果。
现在，我们回过头来看看FPKM和RPKM的计算式，就会发现它们根本做不到！
下面举例说明下。
举例理解
先举个例子来说明（以FPKM的计算为例），我们假定有两个来自同一个个体不同组织的样本X和Y，这个个体只有5个基因，分别为A、B、C、D和E它们的长度分别如下：
| gene | gene length (Kb) |
|---|---|
| A | 100 |
| B | 50 |
| C | 25 |
| D | 5 |
| E | 1 |
由此，我们可以得到，样本X和Y的转录本的不变量，rmean值是1/5。
如果FPKM或RPKM是一个合适的统计量的话，那么，样本X和Y的平均FPKM（或RPKM）应该相等。
我们以FPKM的计算的为例子，以下这个表格列出的分别是样本X和Y在这5个基因分别比对上的fregment数和各自总的fregment数量：
| X | Y | |
|---|---|---|
| A(100) | 80 | 20 |
| B(50) | 10 | 20 |
| C(25) | 6 | 50 |
| D(5) | 3 | 10 |
| E(1) | 1 | 400 |
| Total counts | 100 | 500 |
于是，按照以上公式我们可以得到样本X和Y在这5个基因上的FPKM值分别为：
| X | Y | |
|---|---|---|
| A | 8000 | 400 |
| B | 2000 | 800 |
| C | 2400 | 4000 |
| D | 6000 | 4000 |
| E | 10000 | 800000 |
| Total counts | 100 | 500 |
接下来我们就可以计算FPKM的均值了。我们得到：
样本X在这5个基因上的FPKM均值FPKMmean = 5,680;
样本Y在这5个基因上的FPKM均值FPKMmean = 161,840
虽然理论上来说FPKM是为了描述转录丰度值r，但是不同样本间FPKMmean却根本不同，而且差距相当大，那么究竟为什么会有如此之大的差异？
这是完全是其数学公式所决定的。
我们先来看FPKM的计算公式：
其中：
nf——是比对至目标基因的fregment数量
L——是目标基因的外显子长度之和除以1000（即单位是kb而非bp）
N——是总有效比对至基因组的read数量
首先，我们可以把FPKM的计算式拆分成如下两部分。
第一部分的统计量是对一个基因转录本数量的一个等价描述：nf/L
第二部分的统计量是测序获得的总有效Fregment数量的百万分之一：N/106
FPKM其实就只是这两部分的商！这有道理吗？分开来看它们似乎都有点道理，但是合起来的时候其实很没逻辑。
尤其是第二部分（N/106），本来式子的第一部分是为了描述一个基因的转录本数量，那么正常来讲，第二部分就应该是样本的转录本总数量（或至少是其总数量的等价描述）才能形成合理的比例关系，而且可以看出来FPKM/RPMK是有此意的，这本来就是这个统计量的目的。
但是！！N/10^6的大小其实是由RNA-seq测序深度所决定的，并且是一个和总转录本数量无直接线性关系的统计量——N与总转录本数量之间的关系还受转录本的长度分布所决定，而这个分布往往在不同样本中是有差异的！
这里为了方便后续内容的理解，先在此总结下上面的推断：
将FPKM理解为（nf/L）/（N/106）
nf/L可以理解为对一个基因转录本数量的一个等价描述
理论上来说FPKM应该是用来描述转录丰度的，那么当分子nf/L是对一个基因转录本数量的一个等价描述，分母N/106应该是对总转录本数量的一个等价描述。
总转录本数量的一个等价描述应该要考虑到测序深度及基因长度的影响。但是这里的分母并没有考虑到基因长度的影响。
N必须将各个被转录的基因的长度考虑进去才能正确描述总体的转录本数！而FPKM/RPKM显然没有做到这一点，这便是FPKM和RPKM出错的内在原因。
为什么TPM的对的？
其中：
nf——是比对至目标基因的fregment数量
L——是目标基因的外显子长度之和除以1000（即单位是kb而非bp）
这样的话，当我们计算TPMmean时就可以得到：
这样子，整个统计量就很好理解了，这个值刚好是rmean的106倍，满足等价描述的关系。
下面还是举个实际例子来说明
举例理解
还是用前面的数据：
| X | Y | |
|---|---|---|
| A(100) | 80 | 20 |
| B(50) | 10 | 20 |
| C(25) | 6 | 50 |
| D(5) | 3 | 10 |
| E(1) | 1 | 400 |
| Total counts | 100 | 500 |
我们计算下TPM值：
| X | Y | |
|---|---|---|
| A(100) | 281690.14 | 494.32 |
| B(50) | 70422.54 | 988.63 |
| C(25) | 84507.04 | 4943.15 |
| D(5) | 211267.61 | 4943.15 |
| E(1) | 352112.68 | 988630.75 |
| Total counts | 100 | 500 |
接着，我们可以分别计算样本X和Y的TPM_mean,并且很明显它们都是——200,000
而106/5=200,000
经过这样的标准化之后，X和Y就处于同样的一个标准上了，此刻，在同一个坐标系下的差异比较才是真正有意义的。
为什么还在用FPKM
既然FPKM有错误，那为什么大家还在用？为什么还真找到了（能被实验所验证）的有价值结果呢？
下面内容来自碱基旷工。
PKM/RPKM之所以看起来会是一个合适的值，我想主要原因有三：
其一，它们和TPM之间存在一定的正比关系，如下。通过比较它们各自的数学计算式就可以看出来（以基因A的TPM/FPKM计算为例）：
可以看到，FPKM转TPM的话，对于一个个体来说实际上是一个定值：所有基因的（reads/基因长度）之和，再除以总reads数。
于是我们可以算出：
样本X，RPKM/TPM=0.0284
样本Y，RPKM/TPM=0.8092
本质的原因是，FPKM这个转换会对本来已经正确标准化了的结果，再次做了一次无意义的不等变换，最终导致了结果不可解释。如何理解呢？后文会有补充，这里先简单说一下：这个数学转换式子仅是告诉了我们这样子来计算是可行的，但是在RNA-seq的实际应用场景中，它其实是无生物（或物理）意义的
其二，实验验证的精度是有限的，常用的qPCR也只能给出定性的比较结果，而且实验验证也未必总能成功；
其三，习惯了，形成了路径依赖，大家都这样做，所以也这样跟随，这应该是最不可取的一个情况。
小结
现在回过头来总结一下。事实上，FPKM/RPKM最大的问题就在于其无意义性。我们所要表达出来的任何统计量，它的变化都应该能对应到物理或生物过程中的变化，如果做不到这一点，那么这个统计量往往都是无意义的，用它得到的结果就算看起来符合预期也只不过是数值上的巧合，本质上是不可解释的。
FPKM/RPKM的分母(N/10^6)并不具有任何形式的生物意义，它所能表达出来的这个量，只能代表测序深度的变化，而无法作为对生物基因表达量的描述，即无法代表（等价代表）样本中转录本的总量。
RNA-seq的protoccol差异
主要来自一篇综述的整理内容：
Misuse of RPKM or TPM normalization when comparing across samples and sequencing protocols
我们知道在进行RNA-seq建库测序的时候，我们会有许多种建库方法。但是，因为在人体内所有RNA中，核糖体RNA（rRNA）占到了80-90%的总RNA比例，而我们的测序结果中并不关心rRNA的情况，所以在测序之前，必须从总RNA中去除高度丰富的rRNA。
这就要讲到关于RNA-seq建库的方法，标准方法包括
①使用oligo dT去钓带有polyA尾巴的mRNA
②通过杂交捕获和磁珠分离来消耗高丰度的 rRNA。
这两种方法各有优缺点：
polyA 方法用于大多数转录组研究，因为当焦点主要放在蛋白编码基因上时，所需的测序深度相对较低。
最近的几项研究指出，尽管rRNA耗竭的方法更贵，但是可以提供更多转录的信息。且还有一个优势在于当测定存在降解的mRNA时，这种方法更好。
通过阅读另外一篇文章：
https://www.nature.com/articles/s41598-018-23226-4
作者比较2种方法对人体组织样本和血液样本进行建库，得出以下结论：
rRNA耗竭的方法可以捕获到更多的基因feature。捕获到的read很多会mapping到内含子区域，导致很多reads是无用的，且引起计数的偏差。
polyA方法在基因CDS区的富集程度更高，外显子区覆盖度更高，计数更准。
达到相同的外显子覆盖度，用polyA方法只需要50%的reads，而用rRNA耗竭的方法则需要220%的reads。
作者推荐：如果实验目的是为了做定量（后续做差异分析），更推荐使用polyA方法建库。
其实这篇综述：
http://www.rnajournal.org/cgi/doi/10.1261/rna
也是在总结上面文章的内容，只不过画的图非常直观，下面数据来自同一个体，血液和结肠组织，经过polyA方法和rRNA耗竭方法建库后的RNA-seq数据：
A图：polyA方法测到更多的是protein-coding基因，而rRNA耗竭方法测到更多是small RNA。
B图：不同方法测到表达量最高的TOP3基因占总转录本的比例关系。
https://mp.weixin.qq.com/s/ckoCC2VE_DXwIc-DtwUmVQ