1. 从原始数据到可比数值：为什么我们不能直接用Counts？

刚接触RNA-seq数据分析的朋友，拿到一个基因表达矩阵，第一眼看到的往往是那些密密麻麻的数字，我们称之为“原始计数”（Raw Counts）。这些数字很直接，比如某个基因在样本A里被“数”到了1000次，在样本B里被“数”到了2000次。新手很容易就会想：“哇，这个基因在B样本里的表达量是A样本的两倍！” 但如果你真这么下结论，大概率会掉进坑里。我刚开始做分析的时候，就因为这个踩过雷，结果被审稿人问得哑口无言。

那么，为什么不能直接用Counts来比较基因表达的高低呢？这背后有三个“拦路虎”。第一个是测序深度。想象一下，你拿两个桶去河里舀水，一个大桶（测序深度深），一个小桶（测序深度浅）。大桶舀上来的鱼（Reads）肯定比小桶多。但这能说明大桶对应的那条河鱼更多吗？不能，这仅仅是因为你舀的水量不同。在RNA-seq里，样本A测了5千万条Reads，样本B测了1亿条Reads，即使两个样本中某个基因的真实表达水平一模一样，B样本的Counts值也可能会是A的两倍。第二个是基因长度。基因就像不同长度的“磁铁”，长的基因（比如有10个外显子）在测序过程中，被随机打断并捕获到的片段机会自然比短的基因（只有2个外显子）要多。一个长基因可能仅仅因为“个头大”，就获得了比一个高表达但很短的基因更多的Counts，这显然不公平。第三个是样本间比较。我们做差异表达分析，终极目标是找出不同条件（比如疾病 vs 健康）下，哪些基因的表达发生了真实的变化。如果不剔除测序深度和基因长度的干扰，我们看到的差异很可能只是技术噪音，而非生物学信号。

所以，我们需要对原始的Counts进行“标准化”（Normalization），把它变成一个可以跨基因、跨样本公平比较的数值。这就引出了我们今天的主角：FPKM和TPM。它们都是为了解决上述问题而诞生的“标尺”。简单理解，它们就像是把Counts这个“原始身高”数据，考虑了个人的“腿长”（基因长度）和“测量时的地面平整度”（测序深度），最终换算成了可以公平比较的“标准身高”。接下来，我们就一步步拆解，看看这两把标尺是怎么做出来的，以及它们之间如何自由转换。

2. FPKM详解：它的计算、意义与局限性

2.1 FPKM到底是什么？一个生活化的比喻

FPKM，全称是Fragments Per Kilobase per Million mapped fragments。这个名字听起来很唬人，我们把它拆开，用个生活化的例子来理解。假设你是一个图书管理员，要统计图书馆里不同书籍的受欢迎程度（类比基因表达量）。

• Fragments（片段）：读者每次借阅时，可能会复印书中的几页带回去看。这些复印的“页”就是Fragments。在RNA-seq中，这就是我们测序仪读到的DNA/cRNA片段。对于单端测序（SE），一个Fragment产生一条Read；对于双端测序（PE），一个Fragment产生两条Read，但最终在统计时，一个Fragment（无论对应1条还是2条有效Read）只算作一个计数单位。这就是为什么原始文章里强调要区分Reads和Fragments。
• Per Kilobase（每千碱基）：书有厚有薄。一本《战争与和平》和一本《小王子》，被复印的“页数”直接比较不公平。我们需要除以书的“厚度”（总页数）。在基因里，就是除以基因的长度（以千碱基，Kb为单位）。这样，就消除了基因长度带来的偏差，得到了“每千碱基的片段数”。
• Per Million（每百万）：图书馆每天的总人流量不同，周一可能来了1000人，周六可能来了5000人。人流量大的那天，所有书被复印的总页数自然会多。为了消除“人流量”（即测序深度）的影响，我们需要把每个基因的“每千碱基片段数”再除以当天总复印页数（以百万为单位）。这样，就得到了一个标准化后的指标：每百万总片段中，每千碱基基因长度的片段数。

所以，FPKM的计算公式可以直观理解为： FPKM = （某个基因的片段计数 / 该基因的长度(Kb) ） / （样本总片段计数(以百万计)）

2.2 手把手计算：从Counts到FPKM的代码实战

理论懂了，我们来看看具体怎么算。这里我提供一个非常清晰、分步的R语言代码示例，你可以跟着一步步操作。假设我们有一个表达矩阵 expr，行是基因，列是样本；还有一个基因长度向量 gene_lengths，单位是碱基数。

# 1. 准备示例数据
# 假设有3个基因，2个样本
expr <- matrix(c(100, 200, 300, 400, 500, 600), nrow=3, ncol=2)
rownames(expr) <- c("GeneA", "GeneB", "GeneC")
colnames(expr) <- c("Sample1", "Sample2")
# 基因长度（单位：碱基对）
gene_lengths <- c(1500, 2000, 1000) # GeneA长1.5Kb，GeneB长2Kb，GeneC长1Kb
names(gene_lengths) <- rownames(expr)

# 2. 将基因长度转换为千碱基（Kb）
gene_lengths_kb <- gene_lengths / 1000

# 3. 计算每个基因每千碱基的计数 (Counts per Kilobase)
# 这里需要对每个基因（每一行）进行操作，所以用 sweep 函数，或者直接矩阵除法
# 方法1：矩阵除法（要求长度向量与行对应）
counts_per_kb <- expr / gene_lengths_kb

# 4. 计算每个样本的总片段数（以百万计）
total_counts_million <- colSums(expr) / 1e6

# 5. 计算FPKM
# 关键步骤：将 counts_per_kb 矩阵的每一列，除以对应样本的 total_counts_million
# 使用 sweep 函数可以优雅地实现列操作
fpkm_matrix <- sweep(counts_per_kb, 2, total_counts_million, FUN = "/")

# 打印结果
print("原始Counts矩阵：")
print(expr)
print("FPKM矩阵：")
print(fpkm_matrix)

运行这段代码，你会看到GeneC虽然原始Counts在Sample1中最低（300），但因为它的长度最短（1Kb），其FPKM值可能反而是最高的。这就是长度标准化起的作用。一个常见的坑是：很多人直接用t(t(expr)/colSums(expr)) * 10^9这个公式，其实它等价于上述分步过程，其中10^9就是(1/长度Kb)和(1/总片段数百万)两个分母合并产生的系数。理解分步过程，能让你更清楚每个步骤的意义。

2.3 FPKM的“阿喀琉斯之踵”：为什么它不适合做样本间比较？

FPKM解决了跨基因比较的问题，但它有一个致命的缺陷，这个缺陷在我早期分析中导致过严重的误判。FPKM的标准化是针对单个样本独立进行的。也就是说，每个样本除以自己的总片段数。这听起来没问题，但仔细想：如果样本间有大量表达量极高的“看家基因”，这些基因会吃掉总片段数的很大一部分。当你用每个基因的计数除以这个被“污染”了的总数时，其他基因的FPKM值就会被系统性压低。

更关键的是，FPKM的样本间总和是不固定的。样本A的所有基因FPKM加起来可能是1.2e6，样本B加起来可能是0.8e6。这意味着，你不能直接把样本A某个基因的FPKM和样本B同一个基因的FPKM相除来得到准确的倍数变化（Fold Change）。因为分母的“尺度”不一样。这就像用一把弹性尺子去量不同人的身高，量出来的数值无法直接比较谁更高。因此，在需要进行严谨的样本间差异表达分析时（例如使用DESeq2, edgeR），必须使用原始Counts数据，因为这些工具内置了更复杂、更稳健的标准化算法（如TMM, RLE），能够处理样本间的组成差异。而FPKM数据，如果非要用于差异分析，通常需要先进行log2转换，然后使用limma等基于线性模型的工具，但其统计效力通常不如专门为Counts设计的工具。

3. TPM登场：一个更合理的标准化策略

3.1 TPM vs FPKM：核心思想与计算差异

正因为看到了FPKM的短板，TPM（Transcripts Per Million）被提出并逐渐成为更受推荐的标准化单位。TPM的全称是Transcripts Per Kilobase per Million mapped reads，注意这里最后是“per Million mapped reads”，但它实际的计算逻辑和FPKM有本质区别。

最核心的区别在于标准化的顺序。我们再用图书馆的例子：

• FPKM：先对每本书做“厚度校正”（除以长度），再对每天做“人流量校正”（除以总片段数）。问题是，每天的人流量是原始复印页数，包含了书厚度的影响。
• TPM：先对每本书做“厚度校正”（除以长度），然后把所有书经过厚度校正后的“标准页数”加起来，得到这一天的一个“标准总人气”。最后，用每本书的“标准页数”除以这个“标准总人气”。这样，每天所有书的TPM值加起来，总和都是100万。

计算上，TPM可以看作是FPKM的“再标准化”。公式为： TPM = (某个基因的FPKM / 所有基因FPKM之和) * 10^6

这个简单的再标准化，带来了巨大的好处：TPM使得不同样本间的数值具有了直接可比性。因为每个样本的TPM总和都是相同的（100万），这相当于给所有样本强制统一了“尺度”。样本A中某基因TPM=100，样本B中同基因TPM=200，那么我们可以更有信心地说，在B中的表达量大约是A中的2倍。当然，这仍然是一个相对比较，但比FPKM可靠得多。

3.2 代码实现：三种路径计算TPM

在实际操作中，你可以从原始Counts直接计算TPM，也可以从已有的FPKM矩阵转换。下面我给出两种最常用方法的代码。

方法一：从Counts直接计算TPM（推荐） 这是最直接、误差最小的方法。原理是：先计算长度标准化后的计数（Counts per Kilobase），然后将每个样本的这些数值归一化到百万。

# 接上一节的 expr 和 gene_lengths_kb 数据
# 1. 计算每千碱基计数
counts_per_kb <- expr / gene_lengths_kb
# 2. 计算每个样本的“标准总计数”（即长度标准化后的计数之和）
norm_factor <- colSums(counts_per_kb)
# 3. 计算TPM：将 counts_per_kb 的每一列，除以该列的 norm_factor，再乘以 10^6
tpm_from_counts <- sweep(counts_per_kb, 2, norm_factor, FUN = "/") * 1e6

# 验证：每个样本的TPM总和应为1e6
colSums(tpm_from_counts)

方法二：从FPKM转换到TPM 如果你手头只有FPKM数据，可以用这个方式转换。

# 假设 fpkm_matrix 是上一节计算得到的FPKM矩阵
# 计算每个样本的FPKM总和
fpkm_sums <- colSums(fpkm_matrix)
# 将FPKM矩阵的每一列除以其列和，再乘以 10^6
tpm_from_fpkm <- sweep(fpkm_matrix, 2, fpkm_sums, FUN = "/") * 1e6

# 验证转换结果与方法一是否一致（应基本一致，忽略极小浮点误差）
all.equal(tpm_from_counts, tpm_from_fpkm, tolerance = 1e-9)

方法三：使用对数空间计算函数（数值更稳定） 当数据量极大时，直接在原始空间进行乘除可能会遇到数值计算问题。下面这个函数通过对数操作来提高稳定性，这也是很多成熟软件包内部采用的方式。

Counts_to_TPM <- function(counts, length_kb) {
  # counts: 基因×样本矩阵
  # length_kb: 基因长度向量（单位Kb）
  rate <- log(counts) - log(length_kb) # 等价于 log(counts/length)
  denom <- log(colSums(exp(rate))) # 对每一列计算 log( sum( counts/length ) )
  # exp(rate - denom + log(1e6)) 等价于 (counts/length) / sum(counts/length) * 1e6
  tpm <- exp( sweep(rate, 2, denom, `-`) + log(1e6) )
  return(tpm)
}

# 使用函数计算
tpm_by_function <- Counts_to_TPM(expr, gene_lengths_kb)

4. 实战指南：在不同分析场景中如何选择与转换

4.1 差异表达分析：为什么必须回归Counts？

这是新手最容易犯错的地方。当你已经辛辛苦苦算好了FPKM或TPM，准备做差异分析时，请务必刹车。几乎所有主流的差异表达分析软件（DESeq2, edgeR, limma-voom）的官方文档都会强调：请输入原始计数矩阵（Raw Counts）。

为什么呢？因为这些工具设计的统计模型（如负二项分布）是基于原始计数数据的离散特性。它们内置的标准化方法（如DESeq2的“median of ratios”， edgeR的“TMM”）不仅仅考虑了测序深度，还巧妙地估计了基因间的离散度，并校正了样本间的RNA组成差异。这些算法比简单的FPKM/TPM标准化要强大和稳健得多。如果你把已经标准化过的FPKM/TPM（尤其是经过log转换的）喂给DESeq2，就相当于破坏了数据原有的统计分布假设，得到的结果可能是不可靠的。我自己的教训是，曾经用FPKM做DESeq2分析，发现一些看家基因也出现了显著差异，这明显不符合生物学常识，复查后才发现是数据格式用错了。

那么FPKM/TPM什么时候用？ 它们主要用于：

1. 结果展示与可视化：在论文图表中，展示基因的表达水平时，使用TPM比展示原始Counts更易于读者理解。
2. 样本内基因表达丰度排序：例如，找出某个样本中表达量最高的10个基因。
3. 作为其他下游分析的输入：例如，某些共表达网络分析（WGCNA）建议使用经过合适转换的TPM数据。

4.2 数据可视化与报告：如何优雅地使用TPM

当你用DESeq2做完差异分析，得到了显著差异基因列表后，在绘制热图、小提琴图或展示表达量柱状图时，使用TPM数据会让图表更美观、信息更直观。

一个常见的流程是：

1. 用原始Counts和DESeq2进行差异分析，得到p值、校正后p值以及基于模型估计的标准化计数（如DESeq2::counts(dds, normalized=TRUE)）。注意，这个标准化计数是为了模型内部比较，其数值本身没有像TPM那样的“每百万”解释性。
2. 为了绘图，我们可以单独计算所有基因的TPM值。
3. 提取感兴趣基因（如差异基因）的TPM数据用于绘图。

# 假设差异分析已完成，我们有一个包含所有样本的原始count矩阵 all_counts
# gene_lengths_kb 是已知的

# 计算所有基因在所有样本中的TPM（用于绘图）
all_tpm <- Counts_to_TPM(all_counts, gene_lengths_kb)

# 从差异分析结果中，提取前50个显著差异基因的ID
top50_genes <- rownames(diff_results)[1:50]

# 提取这些基因的TPM数据
plot_data <- all_tpm[top50_genes, ]

# 然后就可以用 plot_data 去绘制热图了
library(pheatmap)
pheatmap(log2(plot_data + 1), # 热图通常对TPM取log2转换以改善视觉效果
         scale = "row", # 按行标准化，突出基因在不同样本间的模式
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         show_rownames = FALSE)

4.3 处理公共数据：当遇到不同标准化单位时

在利用GEO、TCGA等公共数据库数据时，你下载到的表达矩阵可能是多种多样的：RPKM、FPKM、TPM，甚至是已经log2转换后的值。这时，你需要：

1. 仔细阅读数据提交者的说明，确定数据单位。
2. 如果提供的是FPKM/RPKM，并且你需要进行跨样本的定量比较（比如自己再加一批样本做整合分析），强烈建议将其转换为TPM。转换方法就是本章第2节中“从FPKM转换到TPM”的代码。
3. 如果提供的是TPM，并且你只想做差异表达分析，那么很遗憾，你无法直接使用DESeq2/edgeR。一种退而求其次的方法是使用limma包对log2(TPM + 一个小的偏移量)进行分析。但这并非最佳实践，因为TPM失去了原始计数的离散性。更好的情况是，你能从数据库中找到或申请获取原始的Counts数据（很多项目现在都会提供）。
4. 绝对不要混合使用不同标准化方法的数据。千万不要把一部分Counts数据和另一部分FPKM数据直接合并在一起分析，这会导致严重的批次效应和错误结论。必须统一到同一尺度（最好是原始Counts）后再进行整合分析。

5. 常见问题与避坑指南

5.1 基因长度应该怎么取？这是个关键参数

计算FPKM/TPM时，基因长度是一个外部输入参数，它的定义直接影响结果。这里有几个坑需要注意：

• 长度定义：通常使用外显子总长度，而不是基因的起始到终止坐标的长度。因为测序捕获的是成熟mRNA（拼接后的），内含子区域不会被测到。你可以从GTF/GFF注释文件中提取每个基因所有外显子区间的并集长度。
• 使用哪个转录本：一个基因可能有多个转录本（亚型），它们的长度不同。常用的做法是：要么使用最长的转录本的长度来代表该基因（简单，但可能不精确），要么在转录本水平计算FPKM/TPM后再通过某种方法（如求和）聚合到基因水平。工具如StringTie、Cufflinks是在转录本水平进行定量的。
• 获取长度的方法：在R中，你可以使用GenomicFeatures包从GTF文件生成TxDb对象，然后用exonsBy和sum(width)函数来计算基因长度。

library(GenomicFeatures)
# 假设你有hg38的GTF文件
txdb <- makeTxDbFromGFF("Homo_sapiens.GRCh38.xx.gtf", format="gtf")
# 获取按基因分组的外显子
exons_by_gene <- exonsBy(txdb, by="gene")
# 计算每个基因的外显子总长度
gene_lengths <- sum(width(exons_by_gene))
# 注意：gene_lengths 是一个列表，需要转换为向量并与你的基因ID匹配

5.2 我的数据是双端测序（PE），该用FPKM还是RPKM？

原始文章里已经点明了：对于双端测序（PE），使用FPKM；对于单端测序（SE），使用RPKM或FPKM均可（此时两者数值相等）。关键在于你使用的定量软件输出的是什么。例如，Hisat2+StringTie流程输出的就是FPKM；而如果你用featureCounts统计了Reads数，然后自己计算，那么你得到的是基于Reads的计数，此时计算RPKM在数学上没问题，但为了与主流命名一致，很多人也称之为FPKM。在实际应用中，只要理解其原理，知道自己的“计数单元”是Fragments还是Reads，名称的区分并不影响计算和比较。

5.3 标准化后的数值为零或极小怎么办？

在取对数（log2）进行可视化或后续分析时，零值会导致问题（log2(0) = -Inf）。常见的处理方法是加一个伪计数（Pseudocount），通常是加1，即计算log2(TPM + 1)。但要注意，这个“1”是在TPM尺度上加的。对于Counts数据，DESeq2等工具内部会采用更复杂的方法处理零值。在你自己处理时，统一加一个小的值（如1）是一种简便做法，但要知道这会轻微扭曲低表达基因的分布。

5.4 一个完整的自查清单

在你完成计算后，可以用下面这个清单快速检查你的FPKM/TPM数据是否合理：

• [ ] 检查TPM总和：每个样本的TPM值总和是否非常接近1,000,000（允许有极小的浮点误差）？这是TPM计算正确的铁证。
• [ ] 检查高表达基因：表达量最高的基因通常是看家基因（如ACTB, GAPDH），它们的TPM是否在几百到几千的范围内？如果某个样本的看家基因TPM异常低（如<10），可能提示该样本质量有问题或标准化失败。
• [ ] 检查长度相关性：绘制基因表达量（Counts）与基因长度的散点图。在标准化前，两者应有正相关趋势；标准化（FPKM/TPM）后，这个趋势应该被大大削弱。你可以用cor.test(log10(counts+1), log10(gene_length))来检验相关性是否减弱。
• [ ] 与已知结果交叉验证：如果可能，将你的TPM值与类似已发表研究中的表达水平进行粗略比较，看数量级是否相符。

记住，生物信息学分析既是科学也是艺术。理解每个步骤背后的“为什么”，远比记住代码命令更重要。从Counts到FPKM再到TPM，这一路的核心思想就是“公平比较”。希望这篇指南能帮你理清思路，下次再看到这些矩阵时，你能胸有成竹地知道每一个数字背后的故事，并做出正确的选择。