两周前，我接触了一个RNA-seq的项目，做完之后，我重新思考了FPKM和RPKM的计算，觉得它们很可能是不对的，后来查阅了一些文献终于验证了我的想法。现在我重新将这个过程记录了下来：
FPKM和RPKM分别是什么
RPKM是 Reads Per Kilobase Per Million的缩写，它的计算方程非常简单：
$$RPKM = \frac{10 ^ 6 \times n_r}{L × N}$$
其中，$n_r$ 是比对至某一个基因的read数量；$L$是该基因的外显子长度之和除以1000，因此，要注意这里的$L$单位是kb，不是bp；$N$ 是 有效比对至基因组的总read数量。 FPKM是 Fregments Per Kilobase Per Million的缩写，它的计算与RPKM极为类似，如下:
$$FPKM = \frac{10 ^ 6 \times n_f}{L \times N}$$
其中，$n_f$是比对至目标基因的fregment数量。FPKM与RPKM唯一的区别是：F代表fragments，R代表reads。如果是Pair-end测序，每个fragments会由这两个成对的reads构成，因此FPKM只计算两个reads能比对到同一个转录本的fragments数量；而RPKM计算的是可以比对到转录本的reads数量（不管Pair-end的两个reads是否能比对到同一个转录本上）。如果是single-end测序，那么FPKM和RPKM计算的结果将是一致的。
以上是这两个量的计算方式。这样计算的目是为了解决在计算RNA-seq转录本丰度时的两个bias：
（1）相同表达丰度的转录本，往往会由于其基因长度上的差异，导致测序获得的Read（Fregment）数不同。总的来说，越长的转录本，测得的Read（Fregment）数越多，但这并不代表表达量就真的多。
（2）由测序文库的不同大小而引来的差异。即同一个转录本，其测序深度越深，通过测序获得的Read（Fregment）数就越多。
FPKM和RPKM通过同时除以L（转录本长度）和除以N（有效比对的Read（Fregment）总数）的办法，最终将不同样本（或者同个样本在不同条件下）的转录本丰度归一化到一个能够进行量化比较的标准上。
上面的式子看起来似乎合情合理，但是它们却都做错了。
为什么FPKM/RPKM是错的
要回答这个问题，我们需要先撇开所有形式上的计算，重新思考这个问题——到底什么是RNA转录本的表达丰度？事实上，对于任何一个取得的样本，它上面任何一个基因的表达量（或者说丰度），都将已是一个客观存在的值，这个值是不管你改变了多少测序环境都不会变的。而且细胞中此刻总共有多少个基因在表达，实际上也已经是客观定好了的。一旦我们开始以这样一种“先知”的形式来理解的时候，有趣的事情就开始出现了。
此时，我们可以假定，对于样本X，其中有一个基因G被转录了g次，同时样本X中所有基因的转录总次数假定是total，那么正确描述基因G转录丰度的值应该是：
$$r_{g}=\frac{mRNA_g}{mRNA_{total}}$$
没毛病！而且与此同时，样本X中其他基因转录丰度的计算也和以上式子类似，除了要把分子换为其他基因对应的转录次数之外，分母都一样。于是这个有趣的事情就是，所有基因转录本丰度的均值$r_{mean}$将是一个恒定不变的数，由以上定义这个数就是：
$$r_{mean} = \frac{1}{g_{total}}\sum_{g}^{G}{r_g} = \frac{1}{g_{total}}\frac{\sum_{g}^{G}{mRNA_g}}{mRNA_{total}}$$
而
$$\sum_{g}^{G}{mRNA_g} = mRNA_{total}$$
所以
$$r_{mean} = \frac{1}{g_{total}}$$
这个值是由基因的总数决定的，也就是说，对于同一个物种，不管它的样本是哪种组织（正常的或病变的等），也不管有多少个不同的样本，只要它们都拥有相同数量的基因，那么它们的$r_{mean}$都将是一致的。这是一个在进行比较分析的时候，非常有意义的恒等关系。
但在实际的操作中，我们是难以直接计算这些r值的。好在只要能够保证模型的自洽性，我们是能通过自建一些统计量来对r值进行间接描述的，比如FPKM和RPKM。本质上它们的目的就是为了描述r。虽然如此，但我们也要注意，所有这些要用来描述转录本丰度的统计量，都应该能等价描述这一恒等关系。也就是说，不管我们使用了什么统计量，它所描述出来的转录本丰度应该且必须是真实丰度$r_g$的m倍（m必须是一个根据模型定出的不变值），它的均值也将是$r_{mean}$的m倍，至少这样才是得到有意义结果的前提！
（那么）现在，我们回过头来看看FPKM和RPKM的计算式，就会发现它们根本做不到。
举个例子来说明（以FPKM的计算为例），我们假定有两个来自同一个个体不同组织的样本X和Y，这个个体只有5个基因，分别为A、B、C、D和E，它们的长度分别如下：
由此，我们可以得到，样本X和Y的转录本的不变量，$r_{mean}$值都是$r_{mean} = \frac{1}{5} = 0.2$。如果FPKM或RPKM是一个合适的统计量的话，那么至少，样本X和Y的平均FPKM（或RPKM）值应该相等。
我们以FPKM的计算的为例子，以下这个表格列出的分别是样本X和Y在这5个基因中比对上的fregment数和各自总的fregment数量：
于是，按照以上公式我们可以得到样本X和Y在这5个基因上的FPKM值分别为：
接下来就可以计算FPKM的均值了。我们得到，样本X在这5个基因上的FPKM均值$FPKM_{mean} = 5,680$；而样本Y的FPKM均值却是$FPKM_{mean} = 161,840$!! 它们根本不同，而且差距相当大，**那么究竟为什么会有如此之大的差异？**难道这是我故意构造出来的例子所造成的吗？当然不是，**这是由其数学计算上的缺陷所导致的**。
首先，我们可以把FPKM的计算式拆分成两个部分：（1）等价（其实严格来讲也没那么等价）描述某个基因转录本数量的统计量（$\frac{n_f}{L}$） 和（2）测序获得的总有效Fregment数量的百万分之一（$\frac{N}{10 ^ 6}$）；看，FPKM便是这两部分的商！分开来看它们貌似都有点道理，但是合起来的时候其实很没逻辑，尤其是第二部分$\frac{N}{10 ^ 6}$，本来式子的第一部分 （$\frac{n_f}{L}$）就已经是描述某个基因的转录本数量，那么正常来讲，第二部分就应该是描述样本总体的转录本数量（或至少是其等价描述）才能说得通，而且可以看得出FPKM(RPMK)是有此意的，因为这本身就是这一统计量的目的。然并卵，它失败了！$\frac{N}{10 ^ 6}$的大小其实是由RNA-seq的测序深度所决定的，并且是一个和总转录本数量无直接线性关系的统计量——N与总转录本数量之间的关系还受转录本的长度分布所决定，而这个分布往往在不同样本中是有差异的！比如，有些基因，虽然有效比对到它们身上的Fregment数目是相等的，但很明显，长度越长的基因，其被转录的次数就越少。也就是说，N必须将各个被转录的基因的长度考虑进去才能正确描述总体的转录本数！而FPKM（RPKM）显然没有做到这一点，这便是FPKM（RPKM）出错的内在原因。
那么应该是用什么样统计量才合适
其实，通过以上分析，我们已经可以确定一个更加合理的统计量来描述RNA转录本的丰度了。我意外地发现，这个统计量其实在2012年所发表的一篇关于讨论RPKM的文章（RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012.）中就已被提到过了，它被称之为TPM —— Transcripts Per Million，它的计算是：
$$TPM = \frac{\frac{n_r \times read_l}{ g_l} \times {10}^{6}} {T} = \frac{n_r \times read_l \times {10} ^ {6} } {g_l \times T}$$
$$T =\sum_{g=i}^{G}{ (\frac{n_r \times read_l}{g_l})_i }$$
其中，$read_l$是比对至基因G的平均read长度，$g_l$ 是基因G的外显子长度之和（这里无需将其除以1000了）。在不考虑比对剪切的情况下，$read_l$这个值往往都是一个固定值（如100bp或者150bp等），因此我们也可以将$read_l$统一约掉，那么分子就会蜕变成RPKM计算式的第一部分，但把$read_l$留着会更合理。这样，整个统计量就很好理解了，分子是基因G的转录本数（等价描述），分母则为样本中总转录本的数量，两者的比值TPM——便是正确描述基因G的转录本丰度！并且，简单计算我们就可以知道TPM的均值是一个独立于样本之外的恒定值：
$$TPM_{mean} = \frac{10^6}{N}$$
这个值也刚好是$r_{mean}$的$10^6$倍，满足上述等价描述的关系。我们仍然通过上面的例子来进作说明，为简单起见我们只把fregment换为read，其他数字都一样，并且统一假设$read_l$都是一样的：
接着，我们可以分别计算样本X和Y的TPM_mean,并且很明显它们都是$200000 = 10^6 / 5$. 而且，经过这样的标准化之后，X和Y就处于同样的一个标准上了，此刻，彼此之间的比较分析才是真正有意义的。
既然FPKM/RPKM是错的，那为什么大家直到现在都还在用，而且还真找到了（能被实验所验证）有价值的结果呢？
关于对于这个问题，我也思考过。而且我们都知道2008那篇关于RPKM的文章更是用实验结果证明了，RPKM是一个合适的统计量，符合qPCR的验证结果。但归根到底，我觉得眼见未必为实，很多实验其实是表象的，我们更应该从其本质意义和原理上去考虑。FPKM/RPKM之所以看起来会是一个合适的值，我想主要原因有二：
其一，它们和TPM之间存在一定的正比关系。这可通过它们各自的数学计算方程式看出来（以RPKM的计算为例）：
$$RPKM = \frac{T \times 10^3}{ N \times read_l } \times TPM$$
而且在同一个样本内部由于T，N和$read_l$实际上都是定值，因此同个样本内的RPKM和TPM是可以恒等转换的。然而在样本与样本之间就不行，因为不同样本T和N是不同的（假定测序长度$read_l$都一样），这就导致它们之间的转换因子大小不一样！
如以上例子，对于样本X，TPM转换到RPKM的转换因子为：0.0284，但在样本Y中，它的转换因子却是：0.8092。而由于这个基础标准的改变，导致其原本所要描述的“转录本丰度”变得不可比较。然，这其实不是最根本的原因，更本质的原因是，这个转换会对本来已经正确标准化了的结果——TPM，再次做了一次无意义的不等变换，最终导致了结果不可解释。如何理解呢，后文会有补充，这里先简单说一下：这个数学转换式子仅是告诉了我们这样子来计算是可行的，但是在RNA-seq的实际应用场景中，它其实是无生物（或物理）意义的；
其二，实验验证的精度是有限的，常用的qPCR也只能给出定性的比较结果，而且实验验证也未必总能成功。
总结
现在回过头来总结一下。事实上，FPKM/RPKM最大的问题就在于其无意义性。我们所要表达出来的任何统计量，它的变化都应该要能对应到物理或生物过程中的变化，如果做不到这一点，那么这个统计量往往都是无意义的，用它得到的结果就算看起来符合预期也只不过是数值上的巧合，本质上是不可解释的。FPKM/RPKM的分母($N/10^6$)并不具有任何形式的生物意义，它所能表达出来的这个量，只能代表测序深度的变化，而无法作为表达生物过程的量，比如无法代表（等价代表）样本中转录本的总量。
一个统计量该如何计算，说到底都只是一个“术”的问题，而我们应该尽可能在接近其本质意义的地方去确定。
FPKM/RPKM和TPM存在一定的正比关系，因此我们在使用FPKM/RPK时，有些时候确实也能获得可以被实验所验证的“好”结果，但其实它是一个橡皮筋，它的单位刻度是会随着样本的不同而改变的。到头来，样本之间的差异比较实际上也只是在不同的标准下进行的，这样的比较就算得到了所谓的“好”结果，那又有什么意义，根本就是个错误的东东。想想就是由于这种统计量，我们一定已经获得了许多的假阳性结果，同时也肯定错过了许多本来真正有意义的差异，真是弯路走尽也不知，而且还浪费了大堆的心情和时间。
这篇文章：A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Briefings in Bioinformatics.10.1093/bib/bbs046. 对7种主要的RNA-seq标准化方法（但不包含本文提到的TPM）做了一个详细的比较，它用实际结果进行比较（不同于本文所用的数学方式）也得出了RPKM/FPKM这些统计量应该被摒弃的结论，因为它所描述出来的结果是最不合理的，其实所有类似于RPKM/FPKM的统计量在描述转录本丰度的时候都应该被摒弃。
欢迎关注我的个人公众号：helixminer（碱基矿工）
