1. 标准化

1.1. House-keeping gene(s)

1.2. spike-in

1.3. CPM

1.4. TCS

1.5. Quantile

1.6. Median of Ration

1.7. TMM

2. 为什么说FPKM和RPKM都错了？

2.1. FPKM和RPKM分别是什么

2.2. 什么样才算好的统计量

2.3. FPKM和RPKM犯的错

2.4. TPM是一个合适的选择

1. 标准化

目前，转录组测序（RNA-seq）分析是非常成熟的研究手段，有众多分析工具和方法供大家使用，其中，对基因或转录本的读段数目（read count）进行归一化是一个非常重要的分析过程，如何对基因区域进行准确的定量和归一化，是大家十分关心的核心问题之一。

无疑，转录组测序双端数据分析中，目前FPKM是最常用的归一化方法，那FPKM归一化方法是最准确的吗？随着生物信息分析技术的快速发展，FPKM或许已经是“明日黄花”。

归一化的基本背景

总的来说，传统的转录组定量方法是相对定量，一个基因的定量结果很大程度上会受到基因的长度和测序深度的影响。基因长度越长、测序深度越高，得到该基因的read counts就越多，相对表达水平也越高。所以，

在进行下游分析的时候，例如聚类、主成分分析，如果不进行数据归一化直接使用原始read count，简直就是耍流氓。

因此，表达量归一化的精确计算需要同时考虑基因长度、测序深度等信息。

用总reads进行均一化可能最简单，其基于以下两个基本假设：

1、绝大多数的gene表达量不变；

2、高表达量的gene表达量不发生改变；

但在转录组中，通常一小部分极高丰度基因往往会贡献很多reads，如果这些“位高权重”的基因还是差异表达的，则会影响所有其它基因分配到的reads数，而且，两个样本总mRNA量完全相同的前提假设也过于理想了。那如何比较呢，各个方家使出浑身解数，有用中位数的，有用75分位数的，有用几何平均数的，有用TMM(trimmed mean of Mvalues)的等等，总之要

找一个更稳定的参考值。

除FPKM/RPKM外，标准化方法有：

1.1. House-keeping gene(s)

矫正的思路很简单，就是在变化的样本中寻找不变的量

那么在不同RNA-seq样本中，那些是不变的量呢？一个很容易想到的就是管家基因(House-keeping gene(s))

那么 Human 常用的 House-keeping gene 怎么确定？

目前大家用的比较多的一个human housekeeping gene list 来源于下面这篇文章，是2013年发表在 Cell系列的 Trends in Genetics 部分的一篇文章

1.2. spike-in

使用Housekeeping gene的办法来进行相对定量，这种办法在一定程度上能够解决我们遇到的问题。但其实这种办法有一个非常强的先验假设：housekeeping gene的表达量不怎么发生变化。其实housekeeping gene list有几千个，这几千个基因有一定程度上的变化是有可能的。

spike-in方法：在RNA-Seq建库的过程中掺入一些预先知道序列信息以及序列绝对数量的内参。这样在进行RNA-Seq测序的时候就可以通过不同样本之间内参（spike-in）的量来做一条标准曲线，就可以非常准确地对不同样本之间的表达量进行矫正。

spike-in是已知的其他物种基因组，可以对基因表达数据进行绝对值的矫正

个人认为，FPKM/RPKM/TPM 都是基因的相对表达值

DESeq2会对测序深度进行矫正，将普遍高表达的样本A看作是测序深度过高所导致，从而影响差异基因的筛选

RNA的spike-in是一组序列已知的RNA转录本，目前普遍使用的是ERCC( The External RNA Control Consortium)搞出来的一组RNA序列。当然，也可以使用序列相似度较低的物种的序列作为spike-in。如两种常见的酵母，S.pombe和S. cerevisiae，他们的序列可以作为彼此的spike-in进行后续矫正。

通过在样品制备过程中，混入指定数量的spike-in，我们就可以知道不同样本中的基因绝对比表达值。

如等细胞数的样本A和样本B，在每个样本中，我加入了等量的spike-in。最后分析发现，spike-in占样本A的1%，但是占样本B的5%。这表明样本A的RNA表达量也许普遍比样本B的表达量高五倍左右。

下面是一个简单的草图，希望可以帮助理解，左边是细胞，右边是RNA

比较常用的spike-in类型：ERCC Control RNA

ERCC = External RNA Controls Consortium

ERCC就是一个专门为了定制一套spike-in RNA而成立的组织，这个组织早在2003年的时候就已经宣告成立。主要的工作就是设计了一套非常好用的spike-in RNA，方便microarray，以及RNA-Seq进行内参定量

1.3. CPM

CPM(count-per-million)

$\text{cpm}_i=\frac{\text{read count of Gene}_i}{\text{total reads}/10^6}$

1.4. TCS (Total Count Scaling)

简单来说，就是找出多个样本中library size为中位数的样本，作为参考样本，将所有的样本的library size按比例缩放到参考样本的水平

选择一个library size为中位数的sample，以它为baseline，计算出其它sample对于baseline的normalization factor，即一个缩放因子：

$d_i=\frac{S_{baseline}}{S_i}$

然后基于该缩放因子对特定的sample中的每个基因的read count进行标准化（缩放）：

1.5. Quantile

简单来说，就是排序后求平均，然后再回序

在R里面，推荐用preprocessCore 包来做quantile normalization，不需要自己造轮子啦！

但是需要明白什么时候该用 quantile normalization，什么时候不应该用，就复杂很多了

1.6. Median of Ratio (DESeq2)

该方法基于的假设是，即使处在不同条件下的不同个样本，大多数基因的表达是不存在差异的，实际存在差异的基因只占很小的部分， 那么我们 只需要将这些稳定的部分找出来，作为标准化的内参，依据内参算出各个样本的标准化因子

参考

（1）对每个基因计算几何平均数，得到一个假设的参考样本(pseudo-reference sample)

genesampleAsampleBpseudo-reference sample

EF2A1489906sqrt(1489 * 906) = 1161.5

ABCD12213sqrt(22 * 13) = 17.7

…………

（2）对每个样本的每个基因对于参考样本计算Fold Change

genesampleAsampleBpseudo-reference sampleratio of sampleA/refratio of sampleB/ref

EF2A14899061161.51489/1161.5 = 1.28906/1161.5 = 0.78

ABCD1221316.922/16.9 = 1.3013/16.9 = 0.77

MEFV793410570.2793/570.2 = 1.39410/570.2 = 0.72

BAG1764256.576/56.5 = 1.3542/56.5 = 0.74

MOV105211196883.7521/883.7 = 0.5901196/883.7 = 1.35

………………

（3）获取每个样本中Fold Change的中位数，我们就得到了非DE基因代表的Fold Change，该基因就是我们选择的该样本的内参基因，它的Fold Change就是该样本的标准化因子

normalization_factor_sampleA <- median(c(1.28,1.3,1.39,1.35,0.59))

normalization_factor_sampleB <- median(c(0.78,0.77,0.72,0.74,1.35))

1.7. TMM (Trimmed Mean of M value, edgeR)

该方法的思想与DESeq2的Median of Ratio相同，假设前提都是：大多数基因的表达是不存在差异的。

它与DESeq2的不同之处在于对内参的选择上：

DESeq2选择一个内参基因，它的 Ratio/Fold-Change 就是标准化因子。

edgeR 选择一组内参基因集合，它们对标准化因子的计算均有贡献：加权平均。

（1）移除所有未表达基因

（2）从众多样本中找出一个数据趋势较为平均的样本作为参考样本

对所有样本求总Read数；

各样本除以各自的总Read数，得到修正Read数；

求出各自样本修正Read数的Q3值（第3个四分位数）；

所有的Q3值求平均，与平均Q3相差最小的样本即是参考样本。

（3）找出每个样本中的代表基因集，参考这些代表基因集的fold change，计算出该样本的标准化因子

寻找样本的代表基因集：依据基因的偏倚程度和Reads数大小选出——偏倚程度小、reads数居中的基因

衡量偏倚度的量：LFC (log fold change)

$LFC=\log_2\frac{\text{sample}_i}{ref}$

LFC过大或过小都表示具有偏倚性，LFC越大表示reads数在samplei中越高，即偏向samplei；LFC越小表示reads数在ref中越高，即偏向ref

衡量reads数的量：read的几何平均数 (read geometric mean, RGM)

RGM越大表示基因reads越多，RGM越小表示基因reads越少

结合偏倚度、read数画出散点图：

偏倚度小、表达量居中的那些基因落在图中的红线附近

由参考代表基因集计算样本的标准化因子：

对这些代表基因集计算加权平均数：

$\frac{\sum_i^n LFC_i\times \text{ReadCount}_i}{\sum_i^n \text{ReadCount}_i}$

该加权平均数就能代表该样本的标准化因子，只是经过了log变换，所以需要恢复为它的正值：

$\text{ScalingFactor}=2^{\text{weight-average}}$

2. 为什么说FPKM和RPKM都错了？

2.1. FPKM和RPKM分别是什么

在RNA-Seq中，最简单的定量基因表达量（gene expression）的方法就是将RNA-Seq数据比对到相应的参考序列上时，会有比对到各个基因的read数量，称为raw read counts。

但是如果要比较不同样本中基因的表达量，光有raw counts是远远不够的，因为raw cread counts受到很多因素的影响，如目标基因的转录本长度（transcript length）、总的有效比对的read数量（即测序深度 sequencing depth）以及测序的偏差（sequencing bias）等等，这些因素是如何影响raw read counts的后面会有解释。

那么为了将不同样本的基因表达量归一化到一个能够量化比较的标准上，科学家们采取的措施是将raw counts同时除以目标基因的外显子长度之和（也就是目标基因转录本长度）和总的有效比对的read总数。这就是RPKM的定义

RPKM = (10^6 * nr) / （L * N）

其中 nr 代表比对至目标基因的read数量；

L代表目标基因的外显子长度之和除以1000，单位是Kb，不是bp；

N是总的有效比对至基因组的reads数量。

注意这里的 nr：在single-end测序中，一个read就是一个read。而在pair-end测序中，若一对paired-read 都比对上了，当做两个read；若只有一个read比对上，另一个未比对上，当做一个read计算。

类似的，FPKM的定义如下

FPKM = (10^6 * nf) / （L * N）

其中 nf 代表比对至目标基因的fragment数量；

L代表目标基因的外显子长度之和除以1000，单位是Kb，不是bp；

N是总的有效比对至基因组的fragment数量。

注意这里的 nf：在single-end测序中，FPKM将read当做fragment计算，此时FPKM和RPKM是相同的。而在pair-end测序中，若一堆paired-read 都比对上了，当做一个fragment。

RPKM是Reads Per Kilobase per Million的缩写，它的计算方程非常简单：

$RPKM=\frac{10^6 \times n_r}{L \times N}$

FPKM是Fregments Per Kilobase per Million的缩写，它的计算与RPKM极为类似，如下：

$FPKM=\frac{10^6 \times n_f}{L \times N}$

与RPKM唯一的区别为：F是fragments，R是reads，

如果是PE(Pair-end)测序，每个fragments会有两个reads，FPKM只计算两个reads能比对到同一个转录本的fragments数量，而RPKM计算的是可以比对到转录本的reads数量而不管PE的两个reads是否能比对到同一个转录本上。

如果是SE(single-end)测序，那么FPKM和RPKM计算的结果将是一致的。

这两个量的计算方式的目的是为了解决计算RNA-seq转录本丰度时的两个bias：

1、相同表达丰度的转录本，往往会由于其基因长度上的差异，导致测序获得的Read（Fregment）数不同。总的来说，越长的转录本，测得的Read（Fregment）数越多；

2、由测序文库的不同大小而引来的差异。即同一个转录本，其测序深度越深，通过测序获得的Read（Fregment）数就越多。

2.2. 什么样才算好的统计量

首先，RNA转录本的表达丰度是什么？？

事实上，对于任何一个制备好的样本，它上面任何一个基因的表达量（或者说丰度），都将已是一个客观存在的值，这个值是不管你改变了多少测序环境都不会变的。而且总共有多少个不同的基因在表达，实际上也已经是客观定好了的。一旦我们开始以这样一种“先知”的形式来理解的时候，有趣的事情就开始出现了。

对于样本X，其有一个基因G被转录了mRNA_g次，同时样本X中所有基因的转录总次数假定是mRNA_total, 那么正确描述基因g转录丰度的值应该是：

$r_g=\frac{\text{mRNA}_g}{\text{mRNA}_{total}}$

则一个样本中基因表达丰度的均值为

$r_{mean}=\frac{1}{g_{total}} \sum_g^G r_g = \frac{1}{g_{total}} \frac{\sum_g^G \text{mRNA}_g}{\text{mRNA}_{total}}$

而

所以

$r_{mean}=\frac{1}{g_{total}}$

这个期望值竟然和测序状态无关！仅仅由样本中基因的总数所决定的

也就是说，对于同一个物种，不管它的样本是哪种组织（正常的或病变的），也不管有多少个不同的样本，只要它们都拥有相同数量的基因，那么它们的r_mean都将是一致的！

由于上面的结果是在理论情况下推导出来的，实际上我们无法直接计算这个r，那么我们可以尝试通过其他方法来近似估计r，只要这些近似统计量可以隐式地包含这一恒等关系即可。

（那么）现在，我们回过头来看看FPKM和RPKM的计算式，就会发现它们根本做不到。

2.3. FPKM和RPKM犯的错

实际数据来证明

假定有两个来自同一个个体但不同组织的样本（比如同一个病人的肝组织和脾组织 ）X和Y，这个个体只有5个基因，分别为A、B、C、D和E它们的长度分别如下：

值是:

$r_{mean}$

$r_{mean}=\frac 15=0.2$

如果FPKM或RPKM是一个合适的统计量的话，那么，样本X和Y的平均FPKM（或RPKM）应该相等。

以下这个表格列出的分别是样本X和Y在这5个基因分别比对上的fregment数和各自总的fregment数量：

可以算出：

样本X在这5个基因上的FPKM均值FPKM_mean = 5,680;

样本Y在这5个基因上的FPKM均值FPKM_mean = 161,840

很明显，它们根本不同，而且差距相当大

究竟为什么会有如此之大的差异？

可以从其公式上找到答案

首先，我们可以把FPKM的计算式拆分成如下两部分。

$FPKM=\frac{10^6 \times n_f}{L \times N}$

第一部分的统计量是对一个基因转录本数量的一个等价描述（虽然严格来讲也没那么等价）：

$\frac{n_f}{L}$

第二部分的统计量是测序获得的总有效Fregment数量的百万分之一：

$\frac{N}{10^6}$

这么一拆，就可以看出这个公式的问题了：逻辑上根本说不通嘛！

尤其是第二部分（N/106），本来式子的第一部分是为了描述一个基因的转录本数量，那么正常来讲，第二部分就应该是样本的转录本总数量（或至少是其总数量的等价描述）才能形成合理的比例关系，而且可以看出来FPKM/RPMK是有此意的，这本来就是这个统计量的目的。

但是N/106并不能描述样本的转录本总数量。

N/106的大小其实是由RNA-seq测序深度所决定的，并且是一个和总转录本数量无直接线性关系的统计量——N与总转录本数量之间的关系还受转录本的长度分布所决定，而这个分布往往在不同样本中是有差异的！

比如，虽然有些基因，有效比对到它们的Fregment数是相等的，但总的来讲长度越长的基因，其被转录的次数就越少。也就是说，N必须将各个被转录的基因的长度考虑进去才能正确描述总体的转录本数！而FPKM（RPKM）显然没有做到这一点，这便是FPKM（RPKM）出错的内在原因。

FPKM和RPKM通过同时除以L（转录本长度）和N（有效比对的Read（Fragment）总数）的办法，最终将不同样本（或者同个样本在不同条件下）的转录本丰度归一化到一个能够进行量化比较的标准上。

以上一切看起来都很合理

但是！！！

既然说了测序获得的read（fragment）受到基因长度的影响，RPKM和FPKM计算中也去除了目标基因长度的影响，

但是除以N时没有考虑到这个影响，N是总的有效比对的read（fragment），它同样会受到各个转录基因长度（distribution of transcript lengths）的影响。

所以FPKM/RPKM是不准确的。那么有没有一个统计量能解决这个问题呢？有！那就是TPM

2.4. TPM是一个合适的选择

这个统计量在2012年所发表的一篇讨论RPKM的文章（RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012.）中就被提出来了，称之为TPM —— Transcripts Per Million，它的计算是：

TPMi= {( nr/Lr )*10^6 } / sum( nr/Lr+……..+ nm/Lm )

nr：mapping到目标基因上的read数；

Lr：目标基因的外显子长度的总和。

在一个样本中一个基因的TPM：先对每个基因的read数用基因的长度进行校正，

之后再用校正后的这个基因read数 (nr/Lr) 与校正后的这个样本的所有校正后的read数（sum( nr/Lr+……..+ nm/Lm )）求商。

没错！TPM不是除以有效比对的read总数，而是除以经过基因长度归一化后的有效比对的read总数，即归一化后的测序深度。

因此，TPM在计算不同样本的基因表达量比较时，是更加准确的统计量。

$TPM=\frac{\frac{n_r \times read_l}{g_l}}{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$

其中，readl是比对至基因g的平均read长度，gl是基因g的外显子长度之和（这里无需将其除以1000了）。

在不考虑比对的剪切的情况下，readl这个值往往都是一个固定值（如100bp或者150bp等），因此我们也可以将readl统一约掉，那么分子就会蜕变成RPKM计算式的第一部分，但留着会更合理。

这样子，整个统计量就很好理解了，分子是某个基因g的转录本数（的等价描述），分母则为样本中总转录本的数量，两者的比值TPM——便是基因g的转录本丰度！

简单计算之后我们就可以发现 TPM的均值是一个独立于样本之外的恒定值，它等于：

$TPM_{mean}=\frac{10^6}{N}$

这个值刚好是r_mean的一百万倍，满足等价描述的关系。

TPM值就是RPKM的百分比！！！

TPM的优点是什么呢？很明显，所有基因的TPM值加起来肯定是1M，因为百分比的总和就是1嘛，与样本无关，各个样本都可以保证TPM库是一样的，这样比较更有意义！！！

最后还是贴上公式吧！

既然FPKM/RPKM是错的，那为什么大家还在用，而且还真找到了（能被实验所验证）有价值的结果呢？

。而且我们都知道2008那篇关于RPKM的文章用实验结果证明了，RPKM是一个合适的统计量，符合qPCR的验证结果。

但归根到底，眼见未必为实，实验是表象的，我们更应该从其本质意义和原理上去考虑。FPKM/RPKM之所以看起来会是一个合适的值，我想主要原因有二：

其一，它们和TPM之间存在一定的正比关系。这在通过它们各自的数学计算式就可以看出来（以RPKM的计算为例）：

而且在同一个样本内部由于T，N 和 read_l实际上都是定值，因此同个样本内，的RPKM和TPM是可以恒等转换的，只是在样本与样本之间就不行，因为它们之间的转换因子（百分比）大小不一了！

如以上例子，对于样本X，TPM转换到RPKM的转换因子（百分比）为：0.0284，但样本Y的转换因子（百分比）为：0.8092。

而由于这个标准的改变，会导致其原本所要描述的“转录本丰度”变得不可比较。

然，这其实不是最根本的原因，更本质的原因是，这个转换会对本来已经正确标准化了的结果，再次做了一次无意义的不等变换，最终导致了结果不可解释。

如何理解呢，后文会有补充，这里先简单说一下：这个数学转换式子仅是告诉了我们这样子来计算是可行的，但是在RNA-seq的实际应用场景中，它其实是无生物（或物理）意义的；

其二，实验验证的精度是有限的，常用的qPCR也只能给出定性的比较结果，而且实验验证也未必总能成功。

5. 总结

现在回过头来总结一下。事实上，FPKM/RPKM最大的问题就在于其无意义性。

我们所要表达出来的任何统计量，它的变化都应该要能对应到其物理或生物过程中的变化，如果做不到这一点，那么这个统计量往往都是无意义的，用它得到的结果就算看起来符合预期也只不过是数值上的巧合，本质上是不可解释的。

FPKM/RPKM的分母(N/106)并不具有任何形式的生物意义，它所能表达出来的这个量，只能代表测序深度的变化，而无法作为表达生物过程的量，比如无法代表（等价代表）样本中转录本的总量。

一个统计量该如何计算，说到底都只是一个“术”的问题，而我们应该尽可能在接近其本质意义的地方去确定。

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）的统计量在描述转录本丰度的时候都应该被摒弃。