【专题】基因组学技术专题(二)—— 为什么说FPKM/RPKM是错的

RNA-seq-COVER.png

【声明】本文未经许可,禁止转载!
两周前,我接触了一个RNA-seq的项目,做完之后,我重新思考了FPKM和RPKM的计算,觉得它们很可能是不对的,后来查阅了一些文献终于验证了我的想法。现在我重新将这个过程记录下来:

1. 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)总数)的办法,最终将不同样本(或者同个样本在不同条件下)的转录本丰度归一化到一个能够进行量化比较的标准上。上面的式子看起来似乎合情合理,但是它们却都做错了。
 
2. 为什么FPKM/RPKM是错的
要回答这个问题,我们需要先撇开所有形式上的计算,重新思考到底什么是RNA转录本的表达丰度这个问题。事实上,对于任何一个制备好的样本,它上面任何一个基因的表达量(或者说丰度),都将已是一个客观存在的值,这个值是不管你改变了多少测序环境都不会变的。而且总共有多少个不同的基因在表达,实际上也已经是客观定好了的。一旦我们开始以这样一种“先知”的形式来理解的时候,有趣的事情就开始出现了。

此刻,我们可以假定,对于样本X,其有一个基因g被转录了$mRNA_g$次,同时样本X中所有基因的转录总次数假定是$mRNA_{total}$, 那么正确描述基因g转录丰度的值应该是:
$$
r_g  = \frac{mRNA_g}{mRNA_{total}}
$$ 
同时,样本X中其他基因的转录丰度的计算也和以上式子类似,除了要把分子$mRNA_g$换为其他基因对应的转录次数之外,分母$mRNA_{total}$都一样。于是有趣的事情就是,所有基因转录本丰度的均值$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它们的长度分别如下:

RNA-seq-Gene.png

由此,我们可以得到,样本X和Y的转录本的不变量:$r_{mean}$值都是$r_{mean} = \frac{1}{5} = 0.2$。如果FPKM或RPKM是一个合适的统计量的话,那么至少,样本X和Y的平均FPKM(或RPKM)应该相等。

我们以FPKM的计算的为例子,以下这个表格列出的分别是样本X和Y在这5个基因分别比对上的fregment数和各自总的fregment数量:

RNA-seq-fregment.png

于是,按照以上公式我们可以得到样本X和Y在这5个基因上的FPKM值分别为:

RNA-seq-FPKM.png

接下来我们就可以计算FPKM的均值了。我们得到,样本X在这5个基因上的FPKM均值$FPKM_{mean} = 5,680$;而样本Y的FPKM均值却是$FPKM_{mean} = 161,840$!! 它们根本不同,而且差距相当大,那么究竟为什么会有如此之大的差异?难道是因为我故意构造出来的例子所造成的吗?当然不是,这是由其数学计算所决定的。

首先,我们可以把FPKM的计算式拆分成两部分:等价(其实严格来讲也没那么等价)描述某个基因转录本数量的统计量$\frac{n_f}{L}$ 和测序获得的总有效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)出错的内在原因。

3. 那么应该是用什么样统计量才合适
其实,通过以上分析,我们已经可以确定一个更加合理的统计量来描述RNA转录本的丰度了。这个统计量在12年所发表的一篇关于讨论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计算式的第一部分,但留着会更合理。这样子,整个统计量就很好理解了,分子是某个基因g的转录本数(的等价描述),分母则为样本中总转录本的数量,两者的比值TPM——便是基因g的转录本丰度!并且,简单计算我们就可以知道TPM的均值是一个独立于样本之外的恒定值:
$$TPM_{mean} = \frac{10^6}{N}$$
这个值也刚好是$r_{mean}$的$10^6$倍,满足上述等价描述的关系。我们仍然通过上面的例子来进作说明,为简单起见我们只把fregment换为read,其他数字都一样,并且统一假设read_l都是一样的:
 
RNA-seq-TPM.png

接着,我们可以分别计算样本X和Y的TPM_mean,并且很明显它们都是200000 = 10^6 / 5. 而且,经过这样的标准化之后,X和Y就处于同样的一个标准上了,此刻,彼此之间的比较分析才是真正有意义的。
 
4. 既然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是可以恒等转换的,只是在样本与样本之间就不行,因为它们之间的转换因子大小不一了!如以上例子,对于样本X,TPM转换到RPKM的转换因子为:0.0284,但样本Y的转换因子为:0.8092。而由于这个标准的改变,会导致其原本所要描述的“转录本丰度”变得不可比较。然,这其实不是最根本的原因,更本质的原因是,这个转换会对本来已经正确标准化了的结果,再次做了一次无意义的不等变换,最终导致了结果不可解释。如何理解呢,后文会有补充,这里先简单说一下:这个数学转换式子仅是告诉了我们这样子来计算是可行的,但是在RNA-seq的实际应用场景中,它其实是无生物(或物理)意义的;
 
其二,实验验证的精度是有限的,常用的qPCR也只能给出定性的比较结果,而且实验验证也未必总能成功。
 
5. 总结
现在回过头来总结一下。事实上,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这个统计量是应该被摒弃的结论,因为它所描述出来的结果是最不合理的,其实所有类似于RPKM(包括FPKM)的统计量在描述转录本丰度的时候都应该被摒弃。 

http://bib.oxfordjournals.org/content/early/2012/09/15/bib.bbs046.full
http://www.ncbi.nlm.nih.gov/pubmed/22872506


RNA-seq-circos.jpeg

16 个评论

第一次看到这种细致的RPKM/FPKM,TPM的分析,确实存在缺陷错误,楼主分析的不错!
技术贴啊
然而大家用的依旧还是这两个值
特意注册来回复楼主。
楼主我觉得,大家之所以用RPKM,是因为
1.真正的转录总数无法得到,因为你没测到并不一定说那里就没有转录。
2.主要是因为,做分析的时候有个假设,假设两个样品之间转录总次数是一样的。吧。
哈哈。。。
楼主,我觉得有些问题想和你讨论一下。
你在文章里提到的平均转录丰度(rmean)从意义上来看是基因总数的倒数,相当于是一个样品的内秉特性。在下文举的例子中,X,Y两个样本因为基因数量相同,导致有相同的平均转录丰度(rmean),因此你认为两者的平均表达量(FPKM)应该相等。
但是我觉得基因总数相等推导出两者平均表达量相同有些问题。
比如做另一个假设,样品X,Y都有5个基因。样品X中每个基因都只转录了1次,而样品Y中每个基因都转录2次。抛开FPKM的问题,楼主所说的rmean在两个样品中都相同,但按照转录次数来看,Y中基因的平均表达量的确高于X。换句话说,我认为平均转录丰度并和转录本的表达量并没有直接关系。
不知道我理解是否有误,希望楼主针对这个问题解释一下。
一个基因的表达量=该基因的转录数/总的转录数
X样品中的A基因的表达量=1/5=0.2
Y样品中的A基因的表达量=2/10=0.2
应该是这样算吧?
如果两个样品的测序深度相同,使用RPKM做差异分析应该是可以的吧?
不同组织中的mRNAtotal肯定有差异,不能忽略基因表达的时空性;对于mRNAg,在不同组织中也存在差异。如此一来,不同组织的测序文库间差异可能会很明显。如果楼主觉得这个分析在理,回复一下。
除了这里提到的问题,其实rg = mRNAg/mRNAtotal能否用来表征生物学意义上我都不敢苟同。在荧光定量PCR里,真正比较的mRNA表达差异事实上是每个细胞(至少是取材的样品中平均每个细胞)中mRNA表达量的差异,因为在均一化的过程中,选择的reference gene是被认为在每个研究的样本的细胞中表达量基本一致的。而这里,所用到的是mRNAtotal(所有基因的转录总次数)。但是不同的细胞类型,所有基因的转录总次数并不是完全一致的,表达的基因数量也不是完全一致的。这种情况下,mRNAtotal真的能有生物学上的意义吗?不知道我有没有把问题说清楚,打个比方。生长旺盛的幼嫩叶片中,各种基因表达量都很高,假设有3个基因表达,其中h为表达量非常保守的housekeeping基因,在所有细胞中都只表达1份,表达量如下(每个细胞):(a基因:5,b基因:10,h基因:1),那么转录总次数应该为16。而作为对比的衰老叶片,各种基因表达量都偏低(a基因:2,b基因:1,h基因:1),转录总次数为4。5/16 < 2/4就能说明a基因在衰老的叶片中表达量上调了吗?而在荧光定量的实验中,应该是5/1 > 2/1,反而是下调了。专门注册来回复这个问题的,希望能有回复。如果我的想法有什么问题,也请见谅,大家一起讨论: )
好比两次月考,考试题目难度完全一样,全班同学成绩整体下滑,小明在其中因为下滑得没那么多,反而在班上名次上升,能说他进步了吗?不能吧,只能说他退步没那么快。
我个人理解,在illumina建库时,会将不同样品按总量配平,就是说那些平均表达量都低的样品会被人为提高,再加上测序量基本相同,可能这样后就认为转录本总量一致吧
数学不好FPKM , TPM很久都没弄清楚。就把文中例子算了一下(因为有基因长度为1,默认所有读长都是1,其实改成75,200不影响结果),excel过来就这样了,大家可以回去excel里面自己算算看
gene length 100 50 25 5 1 avg sum
sample X 80 10 6 3 1 100
sample Y 20 20 50 10 400 500

FPKM X 8000 2000 2400 6000 10000 5680
FPKM Y 400 800 4000 4000 800000 161840

TPM X 281690.1408 70422.53521 84507.04225 211267.6056 352112.6761 200000
TPM Y 494.3153732 988.6307464 4943.153732 4943.153732 988630.7464 200000
Ti x 0.8 0.2 0.24 0.6 1 2.84
Ti y 0.2 0.4 2 2 400 404.6

那么问题来了,第3个基因根据FPKM算X样品中表达量低于Y样品的;根据TPM算法X样品中表达量高于Y样品的。那到底谁高谁低?
同时,TPM中样品基因表达总量,基因转录总次数mRNA差别远远大于reads差别数目。
而且如果把sample Y 中第二个样品的reads 改为50, 第误个样品改为370,保持总reads数目不变。这样就是同一个基因,在不同样品不同的测序深度样品中占总比例不变,为reference,那么FPKM都是2000一致,TPM 依然差别巨大。如下
gene length 100 50 25 5 1 avg sum
sample X 80 10 6 3 1 100
sample Y 20 50 50 10 370 500

FPKM X 8000 2000 2400 6000 10000 5680
FPKM Y 400 2000 4000 4000 740000 150080

TPM X 281690.1408 70422.53521 84507.04225 211267.6056 352112.6761 200000
TPM Y 533.0490405 2665.245203 5330.490405 5330.490405 986140.7249 200000
Ti x 0.8 0.2 0.24 0.6 1 2.84
Ti y 0.2 1 2 2 370 375.2

总体说来,FPKM更接近理想数据,何解?如果楼主还在,欢迎讨论。
rg = mRNAg/mRNAtotal确实不合理,只有当样品中mRNAtotal变化不大,测序量足够才可以这么用。所以现在大家也是FPKM用的最多。拿上面的数据仔细按照公式计算后跟实验中我们的计算方法也不一致,细节下面的回复。
TPMmean的公式写错了吧,下面的N代表什么???
N是基因的总数目
RPKM确实可以估计r,第一个公式没错,但第二个求r的mean值就有问题了,按照你的模型的话,加一个限制条件:所有基因除了长度不同外,其他完全相同,那么这个模型才可以成立。但是实际情况不是这样,所以RPKM的平均值无法代表r_mean,RPKM更具有实际意义。你的数学模型的第二步是错的,这是一点个人看法,不对的请指教。

要回复文章请先登录注册