常见批次校正方法概览:ComBat、RUV、limma 归纳总结
做组学数据分析的人,谁没被批次效应折磨过?我早期做转录组时就吃过这亏——同一批样本分两次测序,差异基因列表完全不一样。后来才知道,那不是生物学差异,是批次效应在捣乱。所以今天我把三个最常用的校正方法理了一遍,方便大家选型时参考。
1. ComBat——参数化贝叶斯的代表
ComBat 是 Johnson & Li 在 2007 年提出的,核心思路是把每个基因在每个批次里的均值和方差都做一次贝叶斯收缩。它假设基因表达服从正态分布,但允许不同批次有不同的均值和方差。我习惯把它称为「全参数化选手」——模型假设比较强,但只要假设合理,效果非常稳。
import combat
import pandas as pd
# 表达矩阵:基因×样本
expr = pd.read_csv("expression.csv", index_col=0)
batch = pd.Series(["batch1", "batch2", "batch3"] * 5) # 15 个样本
corrected = combat.combat(expr, batch)
print(corrected.shape)
2. RUV——无监督因子分析
RUV 全称是 Remove Unwanted Variation,Gagnon-Bartsch 在 2012 年提出的。它跟 ComBat 思路完全不同——假设批次效应可以用若干隐变量(unwanted factors)来表示,而这些隐变量可以通过「负对照基因」(negative control genes)来估计。说白了,就是找一批理论上不受生物学条件影响的基因,用它们来估算批次结构。
我记得最早用 RUV 处理 scRNA-seq 数据时,实验设计里恰好有一些「批次内重复」的样本,正好可以作为负对照。效果还不错,但要注意负对照的选择非常关键,选错了反而引入噪声。
library(RUVSeq)
set.seed(123)
# 假设 12 个样本,其中前 6 个是 control
counts <- matrix(rnbinom(12000, mu=100, size=1/0.1), nrow=1000, ncol=12)
colnames(counts) <- paste0("sample", 1:12)
# 指定负对照基因(前 100 个)
ctrl_genes <- rownames(counts)[1:100]
# 设定批次因子(这里用随机模拟)
x <- c(rep(1,6), rep(2,6))
# RUVg 方法
ruv_out <- RUVg(counts, ctrl_genes, k=2)
corrected_counts <- ruv_out$normalizedCounts
3. limma——线性模型与经验贝叶斯
limma 不是专门的批次校正工具,它本质上是差异表达分析的线性模型框架。但它有个绝招——在模型里把「批次」作为协变量放进去,然后用经验贝叶斯方法(eBayes)来调整基因水平的方差。这招非常灵活:你可以同时控制年龄、性别、批次等多个协变量。
但要注意:limma 的效果依赖于(1)实验设计可识别;(2)批次与生物条件之间没有完全共线性。我遇到过一位同事把「医院」和「治疗组」完全混杂了,结果 limma 根本无法区分,只能换方法。
library(limma)
expr <- matrix(rnorm(1200), nrow=100, ncol=12)
batch <- factor(c(rep("A",6), rep("B",6)))
condition <- factor(rep(c("control","treatment"), each=3, times=2))
design <- model.matrix(~ condition + batch)
fit <- lmFit(expr, design)
fit2 <- eBayes(fit)
topTable(fit2, coef="conditiontreatment")
4. 其他需要了解的方法
除了上面三个主流方法,还有一些场景下我会想到的备选:
- Harmony:特别适合单细胞数据,用软聚类+迭代校正,对大数据集很友好。
- SVA(Surrogate Variable Analysis):与 RUV 类似,但用奇异值分解来估计隐变量,不需要负对照。
- ComBat-SEQ:专门为 RNA-seq 的计数数据设计的,使用负二项分布。
- Quantile normalization 加批次拟合:简单粗暴,适用于批次效应较弱时。
5. 方法对比总表
| 方法 | 核心假设 | 数据要求 | 优点 | 缺点 |
|---|---|---|---|---|
| ComBat | 正态分布,每批次均值+方差不同 | 批次标签、每批次≥5样本 | 效果稳定,易于使用 | 对分布假设敏感,不适合严重偏态数据 |
| RUV | 负对照基因估计隐因子 | 需指定负对照或重复样本 | 可以捕捉未知批次因素 | 负对照选择困难,不稳健 |
| limma | 线性模型+方差收缩 | 可识别设计矩阵 | 灵活,可同时校正多个协变量 | 批次与生物因素不可完全混杂 |
| Harmony | 软聚类+正交迭代 | ≥2批次,单细胞或高维数据 | 内存友好,无需负对照 | 需要调参(theta, lambda) |
6. 我的选择经验
说实话,没有一种方法能通吃所有场景。我一般按这个顺序快速试:
- 先看能不能找到负对照基因 —— 能的话,RUV 是首选,因为它不需要参数分布假设。
- 如果找不到,且数据近似正态(比如芯片数据),直接上 ComBat,简单又稳定。
- 如果实验设计非常复杂(比如要同时校正性别、年龄、批次),我用 limma 把它放进制模型里一起解。
- 最后别忘了画个 PCA 图看看校正效果。我吃过一次亏——用了 ComBat 后 PCA 看起来干净了,但差异基因却变少了,后来发现是过度校正。所以现在每次校正完,我都会检查一下生物学变异的保留情况。
做批次校正说白了就是「在消除噪声和保留生物学信号之间找平衡」。我建议每换一种校正方法,都保留一份原始数据副本,方便回溯。如果你觉得选型困难,不妨把 ComBat、RUV、limma 都跑一遍,然后用差异基因的重叠率和已知生物学通路富集程度来评判,这也是很多发文章的人用的策略。