4. 表达量归一化：RPKM/FPKM、TPM、CPM 的原理与 Python 实现

拿到 counts 矩阵之后，你第一反应肯定是算差异基因。等等，别急。counts 直接比较？

我给你举个例子。一个基因长 2000 bp，另一个长 200 bp，同样的 100 条 reads 落上去，每 bp 覆盖密度差 10 倍。这不公平吧？还有，不同样品测序深度不一样，1000 万 reads 跟 5000 万 reads 没法直接比数字。

所以，归一化是 RNA-seq 分析里绕不过的坎。我入行时也踩过坑——直接用 raw counts 画热图，结果长的基因全是“高表达”，吓得我赶紧重做。

为什么需要归一化？两个核心维度的校正

• 基因长度：长基因比短基因更容易被读到，必须除以长度。
• 测序深度：文库大小不同，要缩放到统一深度基准。

你想想看，如果不校正，结论很可能只是“长的基因表达高”或“深度的测序样品看起来全部上调”。

下面三种方法分别用不同的方式处理这两个维度。

RPKM / FPKM — 老派方法，但我建议你慎用

RPKM (Reads Per Kilobase per Million mapped reads)

公式很简单：

RPKM = (raw_counts * 1e9) / (gene_length_in_bp * total_mapped_reads)

分子乘 1e9 是为了让结果好看一点（每百万 reads、每千碱基）。

实际上流程是：先除以基因长度（kb），再除以 reads 总数（百万）。

FPKM 是 RPKM 的兄弟版，区别在于 FPKM 用 fragment（双端测序的配对片段）计数代替 read 计数。一个 fragment 两个 reads，FPKM = (fragment_counts * 1e9) / (gene_length * total_fragments)。

TPM — 我的首选，我觉得更合理

TPM (Transcripts Per Million)

TPM 先把 counts 除以基因长度（得到每 bp 密度），再归一化到百万。这样做的好处是：所有基因的 TPM 之和是常数（百万），样本间比较时不受 outlier 基因影响。

公式分两步：

# 第一步：计算每 kb 的 reads 数
reads_per_kb = raw_counts / (gene_length / 1000)

# 第二步：除以总数，乘以 1e6
TPM = reads_per_kb / (sum(reads_per_kb) / 1e6)

或者一个公式：

TPM = (raw_counts * 1e6 * 1000) / (gene_length * sum(reads_per_kb))

注意这个 “sum(reads_per_kb)” 是样本内所有基因的 reads_per_kb 之和，不是 reads 总数。

CPM — 最简单的归一化，只考虑深度

CPM (Counts Per Million)

CPM = (raw_counts * 1e6) / total_mapped_reads

不做基因长度校正。什么场景用？

• 同一个基因在不同样品之间比较（自身长度不变）
• 技术重复 / 批次一致性检查
• 差异分析工具（如 edgeR 的 cpm 函数）内部会再配合其他归一化步骤

但跨基因比较？不行，因为长基因天然 CPM 高。

三种方法速查表

方法	长度校正	深度校正	总和特性	推荐场景
RPKM	是（千碱基）	是（百万 reads）	不固定	已过时，尽量不用
FPKM	是（千碱基，fragment）	是（百万 fragments）	不固定	同 RPKM，已过时
TPM	是（千碱基）	是（百万 TPM 单位）	固定为 1e6	样品间比较、PCA、热图
CPM	否	是（百万 counts）	固定为 1e6	技术重复、批次检查、差异分析预处理

Python 实现：一行代码都不浪费

假设你已经有 counts 矩阵（基因×样品），以及基因长度字典 gene_lengths。

import pandas as pd
import numpy as np

def calculate_tpm(counts_df, gene_lengths):
    """
    counts_df: DataFrame, index=gene_id, columns=sample_id, values=raw counts
    gene_lengths: dict or Series, {gene_id: length_in_bp}
    return: TPM DataFrame
    """
    # 确保长度匹配
    lengths = pd.Series(gene_lengths).reindex(counts_df.index)
    if lengths.isnull().any():
        raise ValueError("有基因缺少长度信息")
    
    # 每 kb 的 reads 数
    reads_per_kb = counts_df.div(lengths / 1000, axis=0)
    # 缩放至百万
    scaling_factor = reads_per_kb.sum(axis=0) / 1e6
    tpm = reads_per_kb.div(scaling_factor, axis=1)
    return tpm

def calculate_cpm(counts_df):
    """
    仅深度归一化
    """
    total_counts = counts_df.sum(axis=0)
    cpm = counts_df.div(total_counts / 1e6, axis=1)
    return cpm

# 示例使用
gene_len = {'BRCA1': 80000, 'TP53': 26000, 'MYC': 15000}
counts = pd.DataFrame({
    'sample1': [100, 200, 50],
    'sample2': [150, 180, 60]
}, index=['BRCA1', 'TP53', 'MYC'])

tpm_df = calculate_tpm(counts, gene_len)
cpm_df = calculate_cpm(counts)
print(tpm_df.round(2))
print(cpm_df.round(2))

核心逻辑一览图

下面这张 SVG 画出从 raw counts 出发，根据你的分析目标选择归一化路径：

避坑总结

• 基因长度从哪里拿？ 绝对不要用基因组坐标跨度的长度（内含子很大）。去 GTF 文件里提取外显子合并后的转录本长度，或者用生物信息工具（如 pyensembl、biomart）获取。
• 零 counts 怎么办？ 按公式算出来是 0，没问题。但如果总 reads 里异常高表达基因占比极大，会导致其他基因 TPM 偏低。这时考虑过滤掉核糖体 RNA 等干扰基因。
• 是否要对数转换？ 我建议 PCA 或热图时做 log2(TPM + 1) 或 log2(CPM + 1)，让分布接近正态，不然高表达基因把低表达的压到看不见。
• 批次效应能靠归一化消除吗？ 不能。归一化只能消除技术差异（长度、深度），批次效应需要专门的工具，比如 ComBat 或 RUVseq，那是后几章的内容。