基因表达定量与归一化方法：RPKM、FPKM、TPM

RNA-seq的下机数据拿到手，你准备怎么量化基因的表达量？这个问题看似简单，但选错方法，后面的分析结果就会跑偏。

我刚开始做RNA-seq分析时，也被这个问题折腾过。原始counts受两个主要因素影响：基因长度和测序深度。长基因比短基因更容易被检测到，测序深度高的样本自然计数也高。说白了，如果不做归一化，你根本分不清一个基因表达高，到底是因为它真的活跃，还是因为它长得长。

为什么需要归一化？

打个比方。你去两个医院做体检，一个医院只查了10项指标，另一个查了100项。如果只看阳性项目数，第二个医院肯定更多。但你得按比例算，对吧？RNA-seq也是一样。每个样本的测序深度不同，每个基因的外显子长度也不同。直接拿counts比较，就像拿苹果和橙子比重量。

归一化就是解决这个问题的。它把不同样本拉到同一个基准线上。我见过不少新手，拿到counts就直接做差异分析，一查结果，全是长基因。嗯，这就是没归一化闹的。

RPKM与FPKM：历史产物

RPKM（Reads Per Kilobase per Million）是最早被广泛使用的归一化方法之一。它先按基因长度校正，再按测序深度校正。公式其实不复杂：

RPKM = (reads_count / gene_length_kb) / total_mapped_reads_million

举个例子。一个基因长2kb，有100个reads，总测序深度是10M。那RPKM就是 (100/2) / 10 = 5。听起来挺合理吧？

但我得说句实话：我在项目中遇到过，RPKM对低表达基因波动很大。尤其是那些只有一个外显子的小基因，稍微有点背景噪音，RPKM就飙上去了。后来FPKM（Fragments Per Kilobase per Million）被提出来，主要是针对双端测序。双端测序产生的是fragment，不是单个read，所以FPKM用fragment counts代替reads counts。

我记得第一次手动算RPKM和FPKM的差异，差点把自己绕进去。其实核心区别就一句话：单端用RPKM，双端用FPKM。但要注意，两个方法都假设了样本间的表达分布一致。这个假设在比较不同组织类型时往往不成立。你想想看，肿瘤样本和正常样本的转录组构成差异很大，用RPKM或FPKM做跨样本比较，风险不小。

我曾经在一个项目中，用FPKM比较了50个肿瘤样本和50个正常样本。结果发现，所有正常样本的管家基因表达量都特别高，肿瘤样本的管家基因反而低。一开始我以为是生物学差异，后来发现是FPKM假设踩坑了——正常样本的mRNA比例更高，RPKM/FPKM在总reads归一化时反而压低了肿瘤样本的相对表达量。

TPM：当前标准

TPM（Transcripts Per Million）为啥成了目前推荐的标准？说白了，它的设计思路跟RPKM/FPKM正好反过来。TPM先按基因长度归一化，再按每百万总转录本归一化。

# 计算TPM的步骤
rate = counts / gene_length
tpm = rate / sum(rate) * 10^6

关键区别在这里。RPKM是先按深度再按长度，TPM是先按长度再按深度。你可能会问：顺序不一样，结果能差很多？答案是，差别很大。TPM把每个基因的表达量看作一个比例，比例的总和固定为1M。这样一来，不同样本之间的TPM可以直接比较。一个基因在样本A中是100 TPM，在样本B中是50 TPM，那它就是在A中表达更高——这个结论是成立的。

我个人习惯用TPM做后续的聚类和分类。在肿瘤亚型分类项目中，我用TPM作为输入特征，聚类结果更稳定。记得有一次，我拿到两个不同实验室的RNA-seq数据，一个用的RPKM，一个用的FPKM。为了合并数据集，我花了不少时间把它们都转成了TPM格式。转完之后，两个批次的样本终于能放在一起看了。

三种方法对比

方法 归一化思路 适用场景 特点
RPKM 深度→长度 单端测序、样本内比较 历史最久，跨样本能力弱
FPKM 深度→长度 双端测序、样本内比较 与RPKM类似，考虑fragment
TPM 长度→比例 跨样本比较、聚类、分类 当前推荐，可比性好

你可能会问：那差异表达分析为什么不直接用TPM？嗯，这里要注意。差异表达分析工具像DESeq2和edgeR，它们内部有自己的归一化方法（比如DESeq2的中位数比值法），更适合检测微小差异。如果你已经用TPM，再输入DESeq2，反而会造成统计上的问题。我个人做差异表达时，会用counts输入DESeq2，做聚类和分类时，才会用TPM。

归一化流程

从上游定量到下游归一化，我通常走这个流程：

  1. 上游定量：用kallisto或salmon直接获得counts和TPM。这两个工具速度很快，还自动处理了多映射问题。
  2. 得到counts：用于差异表达分析，输入DESeq2或edgeR。
  3. 得到TPM：用于聚类、分类和可视化。

有个小技巧：如果上游用的是STAR + featureCounts，你也可以自己算TPM。代码很简单：

# 用R计算TPM的示例
library(data.table)
dt <- fread("gene_counts.txt")
dt[, rate := counts / gene_length_kb]
dt[, tpm := rate / sum(rate) * 1e6]

我曾经在项目里试过直接用kallisto输出的TPM。一个基因有多个转录本，kallisto会给出每个转录本的TPM，然后按基因汇总。这里要注意，多个转录本的TPM不能直接加和，因为TPM是比例值。应该先把每个转录本的counts加总，再重新计算基因级别的TPM。

如果你需要合并多个基因表达矩阵，建议全部转成TPM格式再合并。我通常先用tximport把kallisto的结果导入R，然后手动转TPM。tximport自带一个 countsFromAbundance 参数，选 "lengthScaledTPM" 可以直接搞定。

实际案例：以乳腺癌为例

在我们这个肿瘤亚型分类课程中，数据来自TCGA的乳腺癌样本。上游定量用的是kallisto。每个样本得到一个TPM矩阵，行是基因，列是样本。我检查了一下数据的分布情况：

  • TPM范围从0到几万，大部分基因集中在0.1~100之间。
  • 线粒体基因的TPM普遍偏高，通常会占到总转录本的10%左右。
  • 低表达基因（TPM < 1）占了将近一半。

我记得第一次看到这个分布时，心想：这么多低表达基因，要不要过滤掉？后来实践证明，保留低表达基因反而增加了噪声。我建议在分类之前，至少过滤掉那些在80%以上样本中TPM < 1的基因。

另外，TPM非常适合做log2转换。log2(TPM + 1)可以把偏态分布拉成正态，方便后续的PCA和聚类分析。

避坑指南

  • 不要混用不同归一化方法的数据：我曾经在分析肿瘤亚型时，从两个数据库下载了表达矩阵。一个用的RPKM，一个用的TPM。我直接合并，结果聚类全按批次分开了。后来统一转成TPM才搞定。
  • TPM用于差异表达要谨慎：一个学生用TPM跑DESeq2，发现结果不稳定——因为DESeq2内部默认输入是counts，TPM的小数点破坏了它的离散度估算。
  • 基因长度变化的影响：如果你比较的是同一物种的不同条件，基因长度不变，影响不大。但比较跨物种时，基因长度差异会干扰RPKM/FPKM的结果。TPM不受这个问题影响。
一句话总结：做肿瘤亚型分类，TPM是你最好的朋友。聚类、PCA、热图、分类特征，统统用TPM。差异表达分析才用counts + DESeq2。别搞混。
← 上一章 📖 返回目录 下一章 →