理论 | 简述RPKM，FPKM，TPM
在高通量测序数据的分析中，仅仅靠raw read counts描述基因的表达量是远远不够的。受限于测序过程中的技术因素影响，read counts对于基因表达量的反映存在一定偏好（bias）。因此，Mortazavi等人提出了RPKM/FPKM的方法对read counts进行normalization以使基因表达量的比较可以在不同文库间进行。随后，Mortazavi等人更是提出了考虑转录本长度分布情况的TPM方法。本文将会简要说明为什么我们要对read counts进行normalization，以及RPKM，FPKM，TPM是什么，并通过一个简单地例子阐述为什么TPM才是被更多人认同的方法。
为什么我们要进行Normalization
我们之所以说测序得到的read counts并不是其mRNA丰度的忠实反映，是因为read counts的数量会受到多种因素的影响，例如：
-
测序深度：某些低表达量的基因只有在较高的测序深度时才能检测到。一般而言，随着测序深度的增加，基因种类以及可变剪接体的数目也会增加。同时，测序深度高的样本read counts也会较高。
Effect of sequencing depth
在上图中，样本A中的基因表达量是样本B的两倍，但这是由于测序深度引起的结果，而非真实存在的差异。
-
基因长度：由于高通量测序是将cDNA碎片化后再进行测序的，因此越长的基因产生的碎片也会更多，在测序中也会更加容易被检测到。所以对基因长度的校正也是十分有必要的。
Effect of gene length
在这个图里，基因X和Y的真实表达量是一致的，但是基因X的reads会比基因Y要多，这是由于基因X的基因长度较长所致的。
除了上述两个主要因素外，还会有其他因素对read counts的检测有所影响，例如转录组的组成，GC含量，random hexamers引起的测序偏好等等。由于上述因素的存在，导致在不同样本间使用read counts 进行比较是不太现实的，人们便提出了许多对read counts进行Normalization的方法。本文在此简单地介绍使用最为广泛却最受质疑的RPKM/FPKM，以及受到更多人认可的TPM。
RPKM/FPKM
RPKM: Reads per kilo base per million mapped reads
FPKM: Fragments per kilo base per million mapped reads
RPKM
从RPKM的公式中我们可以看到，RPKM对基因长度（gene length）以及测序深度（mapped reads from library）都进行了校正。而FPKM只是RPKM的双端测序（pair-end）版本。
TPM
TPM: Transcript per million
TPM的计算公式如下所示：
TPM
同RPKM一样，TPM对基因的长度进行了校正，计算比对到基因上的reads/基因长度得到长度校正的表达量 reads per kilobase (RPK)。再以文库中RPK之和作为Scale Factor求出TPM。
相比于RPKM使用read counts之和来作为文库校正因子，TPM使用RPK之和作为文库校正因子的好处是考虑了不同样本间的基因长度的分布。因为RPK是一个对基因长度进行校正后的表达量单位，所以RPK之和也不会再带入基因长度的bias。因此，如果需要比较的样本之间转录本分布不一致时（例如不同物种RNA-seq的比较），使用TPM是一个较佳的Normalization方案。
一个例子
以下有ABCD四个基因，并同时进行了三次重复的测序（rep1，rep2，rep3）。首先，可以看出由于基因长度的关系，B基因的read counts都是较其余的基因较高的（排除D而言）。其次，可以看出rep3的测序深度较高，得到的read counts也较多，同时还检测到其余两个重复中没有检测到的D基因（低丰度基因）。不管是测序深度还是基因长度的严重地干扰了我们对不同样本的基因表达量比较。因此，我们采用RPKM和TPM的方法进行Normalization后，再来比较其中的差异。
为了方便起见，以下的函数省略了kilo base，million base的转换
RPKM <- function(x, gene.length){
lib.size <- sum(x)
rpm <- x/lib.size
rpkm <- rpm/gene.length
return(rpkm)
}
TPM <- function(x, gene.length){
rpk <- x/gene.length
scale.factor <- sum(rpk)
tpm <- rpk/scale.factor
return(tpm)
}
expr <- data.frame(rep1=c(10,20,5,0), rep2=c(12,25,8,0), rep3=c(30,60,15,1),
row.names = c("A", "B", "C", "D"))
gl <- c(2,4,1,10)
rpkm1 <- apply(expr, 2, FUN = RPKM, gene.length = gl)
tpm1 <- apply(expr, 2, FUN = TPM, gene.length = gl)
# 对结果都乘1000，看起来直观点
rbind(rpkm1*1000, colSums(rpkm1)*1000)
## rep1 rep2 rep3
## A 142.8571 133.3333 141.5094340
## B 142.8571 138.8889 141.5094340
## C 142.8571 177.7778 141.5094340
## D 0.0000 0.0000 0.9433962
## 428.5714 450.0000 425.4716981
rbind(tpm1*1000, colSums(tpm1)*1000)
## rep1 rep2 rep3
## A 333.3333 296.2963 332.594235
## B 333.3333 308.6420 332.594235
## C 333.3333 395.0617 332.594235
## D 0.0000 0.0000 2.217295
## 1000.0000 1000.0000 1000.000000
通过以上的例子，我们可以发现RPKM和TPM最重要且直接的区别就是，TPM校正后的表达量之和在不同样本中是相等的，而在RPKM中并不是。这就导致当我们使用RPKM在样本间进行比较的时候就不能得到可靠的结果。
虽然以上的例子较为理想化，但在进行实验时，我们认同的一个前提是不同的技术重复（例如一个细胞重复测三次）所得到的基因表达量应当是一致的。在测序的过程中或许会引入一定的误差，但这些误差是可以通过Normalization的方法在一定程度上校正回来的。简而言之，RPKM和TPM这类方法就是为了使不同样本间的总体表达量趋于一致，让不同样本间的基因表达量有可比较性，而TPM在这方面表现得较为出色，故受到学界更多的青睐。
- 常用的Normalization 方法总结
| Normalization method | 描述 | 考虑因素 | 使用场景 |
|---|---|---|---|
| CPM (counts per million) | 使用read counts的总和校正counts | 测序深度 | 同一样本组重复样本之间的基因counts比较；不适用于样品内的比较或差异分析 |
| TPM (transcripts per kilobase million) | 每百万mapped reads中每kb转录本上的reads数 | 测序深度和基因长度 | 样本内或同一样本组样本之间的基因counts比较； 不适用于差异分析 |
| RPKM/FPKM (reads/fragments per kilobase of exon per million reads/fragments mapped) | 如TPM | 测序深度和基因长度 | 样本内基因间的counts比较; 不适用于样本间比较和差异分析 |
| DESeq2’s median of ratios [1] | counts除以样本特异的校正因子，该因子由基因计数相对于每个基因的几何平均值的中位数比率确定 | 测序深度及RNA组成 | 样本之间的基因counts比较以及差异分析; 不适用于样本内比较 |
| EdgeR’s trimmed mean of M values (TMM) [2] | 使用样本之间的加权截尾的对数表达量比值的均值进行TMM校正 | 测序深度，RNA组成以及基因长度 | 样品之间和样品内部的基因counts比较，适用于差异分析 |
这边补充一点是在差异表达分析中，对基因表达量的Normalization没有这么严谨，只是要求对文库大小进行校正即可，例如
edgeR包中的cpm(counts per million)只是简单地将read counts/library size。对这方面感兴趣的同学可以看我介绍
edgeR的那篇文章：简单的转录组基因差异表达分析 -- edgeR其中有更详细的介绍。
参考文章:
- https://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/
- https://www.biostars.org/p/273537/
- https://www.jianshu.com/p/1940c5954c81
- Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nature methods. 2008 Jul 1;5(7):621-8.
- https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html
完。
