ComBat 的定位：不是黑箱，是贝叶斯

我接触批次校正这么多年，最常被问到的问题就是：「那么多方法，为什么大家这么爱用 ComBat？」
说白了，ComBat 并不复杂。它抓住了批次效应的本质：不同批次导致的系统性偏移，既有位置上的平移，也有尺度上的缩放。你想想看，两个实验室做 RNA-seq，测序深度不同、试剂批次不同，得到的表达值分布能一样吗？肯定不一样。

ComBat 的聪明之处，在于它用经验贝叶斯框架来估计这些偏移。不是硬邦邦地「减均值除标准差」，而是从数据自身「学习」批次效应的统计特征。这带来一个好处：即使某个批次只有少量样本，先验信息也能帮忙稳住估计结果。

一句话概括：ComBat 假设每个基因在批次间的差异服从正态分布，通过经验贝叶斯「借用」所有基因的信息来稳健估计两个参数——位置参数（均值偏移）和尺度参数（方差变化）。

经验贝叶斯框架：怎么个「借用」法？

我在刚接触 ComBat 时，也被这个「经验贝叶斯」绕晕过。后来我换了个角度理解：
传统方法，比如简单均值中心化，每个基因独立校正。这在小样本批次里很容易过拟合——一个基因跑了 3 个样本，均值估计本身就不可靠，再去做校正，结果还不如不校。
ComBat 怎么做？它先把所有基因的方差估计看作一个整体分布，然后「借用」所有基因的信息来缩小每个基因的方差估计。嗯，这里要注意，它并没有假设所有基因的方差一样，只是假设它们来自同一个分布族。

别混淆：经验贝叶斯 ≠ 完全贝叶斯。ComBat 先通过数据（经验）得到先验分布的参数，再用这些先验做后验估计。它没有马尔可夫链蒙特卡洛采样，计算效率高得多。我见过一些朋友跟审稿人解释时写「完全贝叶斯」，审稿人马上指出问题。

具体步骤大致是这样：

  1. 分步估计：先拟合一个线性模型，去掉生物学变异（比如处理组 vs 对照组），留下残差。
  2. 假设分布：每个基因在残差中的批次效应，假设服从正态分布：均值 ~ Normal, 方差 ~ Inverse Gamma。
  3. 学习先验：用所有基因的方差（或均值）估计，拟合先验分布的参数
  4. 收缩估计：用贝叶斯后验公式，把每个基因的估计「拉向」整体的先验中心。
  5. 校正：用收缩后的参数重新对齐每个批次的数据。

这个「收缩」过程，就是 ComBat 的精华。我记得一个做芯片数据的老同事跟我说，他用 ComBat 之前总是要手动剔除掉方差最小的 20% 基因，因为那些基因的批次效应估计太吵。ComBat 直接把这个工作内化成了统计公式。

ComBat 经验贝叶斯核心流程 输入：表达矩阵 + 批次标签 1. 拟合线性模型 2. 计算残差 & 批次矩量 3. 经验先验估计 4. 后验收缩参数 5. 校正原始表达矩阵 关键理解 「借用」所有基因 的信息来稳定估计

参数详解：你不能乱调的开关

ComBat 的默认参数，在 90% 的情况下已经够用。但剩下 10% 的情况——比如芯片数据混了 RNA-seq、或者批次效应极度不均匀——就需要动参数了。我先列一个参数速查表，再逐个讲。

参数默认值作用我的使用建议
mean.onlyFALSE是否只校正均值，不校正方差RNA-seq 数据强烈建议设 TRUE
par.priorTRUE是否使用参数先验（Normal+InvGamma）中等样本量（每个批次>=3）用 TRUE；小样本用 FALSE
prior.plotsFALSE是否诊断检验先验拟合质量第一次建模时务必开一次，后面可以关
ebTRUE是否执行经验贝叶斯收缩几乎不动，关了就是普通位置尺度校正
survNULL保留的生存/协变量（不参与校正）有临床协变量时传一个数据框

mean.only：最容易踩坑的参数

为什么我要把 mean.only 放在第一个讲？因为我在项目中吃过亏。那时候合并两个 TCGA 批次，RNA-seq 数据，我默认跑了 ComBat，结果下游差异分析调出来一堆假基因。一查，原来是 ComBat 把方差也强行拉齐了。
RNA-seq 数据的方差通常与均值相关（高表达基因方差大，低表达基因方差小）。如果连方差一起校，可能会抹掉真正的生物学高变异基因。
所以我的习惯是：做 RNA-seq 数据（包括单细胞拟批量），第一件事先把 mean.only 设成 TRUE。 

# R 代码示例
library(sva)
adjusted_data <- ComBat(
  dat = expression_matrix,
  batch = batch_vector,
  mod = model.matrix(~ condition),
  mean.only = TRUE     # RNA-seq 推荐
)

技巧：如果你不确定是否该用 mean.only，可以快速跑两个版本，然后比较校正后数据的主成分图。如果 mean.only=TRUE 版本的 PCA 上批次之间的重合度已经很好，就不要再加方差校正了。

par.prior：什么时候用非参数？

par.prior = TRUE 使用参数化的先验分布（正态 + 逆伽马）。这需要数据量中等以上。如果你的某个批次只有 2 个样本，参数先验估计出来的方差可能不可靠。这时候换成 par.prior = FALSE，ComBat 会用非参数估计——说白了，就是直接用经验分布的分位数来估计。

我曾经在一个单细胞数据合并的项目里，某个批次只有 1 个样本（是的，你没看错，用户非要加一个外部参考）。我试了 par.prior = FALSE 才没有报错。不过说实话，这种极端情况我建议你还是别用 ComBat 了，直接用 Harmony 或者 scVI 更稳妥。

prior.plots：帮你判断要不要调 par.prior

这个参数被很多人忽略。我建议每个项目至少开一次 ──── 它会画两张图：第一张显示基因方差在对数尺度上是否接近正态，第二张显示剩余方差的正态 Q-Q 图。如果图上的点严重偏离对角线，说明你当前的先验假设不合理，可能需要换 par.prior = FALSE。

# 开启诊断图
adjusted_data <- ComBat(
  dat = expression_matrix,
  batch = batch_vector,
  prior.plots = TRUE
)
# 图片会在 R 的绘图设备中显示
# 如果 Q-Q 图整体呈 S 形，考虑 par.prior = FALSE

eb：开关的取舍

eb = TRUE 是灵魂。如果把它设成 FALSE，ComBat 就退化成了标准化的位置-尺度校正——每个基因独立估计，不再借用信息。我只有在一种情况下关掉它：直接对比「用 eb」和「不用 eb」对下游分析的影响，用来给审稿人证明经验贝叶斯是必要的。日常分析中，我没有关过 eb。

实战案例：TCGA 泛癌数据合并

我记得有个项目需要合并 TCGA 的 BRCA 和 LUAD 两个癌症数据，它们由不同测序中心产出，批次效应非常明显。我做了几步：

  1. 预处理：FPKM 转 TPM，取 log2(TPM + 1)。这一步是为了让数据更接近正态，满足 ComBat 的假设。
  2. 设置协变量：把癌症类型（BRCA vs LUAD）作为 mod 传入，防止 ComBat 把真实的生物学差异当成批次校掉。
  3. 启动校对：使用 mean.only = TRUE, prior.plots = TRUE 先跑一次，看诊断图没问题后，关掉 prior.plots 跑最终版本。
# 完整示例
library(sva)
expr <- readRDS("merged_tcga_logtpm.rds")
batch <- ifelse(grepl("BRCA", colnames(expr)), "Batch1", "Batch2")
cancer_type <- ifelse(grepl("BRCA", colnames(expr)), "BRCA", "LUAD")

# 协变量矩阵（保护生物学差异）
mod <- model.matrix(~ cancer_type)

# 第一次跑，开启诊断
adjusted <- ComBat(dat = expr, batch = batch, mod = mod,
                   mean.only = TRUE, prior.plots = TRUE)

# 诊断图通过后，关闭诊断，跑正式版本
adjusted_final <- ComBat(dat = expr, batch = batch, mod = mod,
                         mean.only = TRUE, prior.plots = FALSE)

结果：校正前 PCA 第一主成分完全按批次分开（方差占比 38%）；校正后批次不再成为主要变异来源，但 BRCA 和 LUAD 的差异依然保留。后来这篇文章投出去，审稿人对批次处理没有任何质疑。

避坑指南

  • 别把生物学变异扔进 mod：我见过有人把处理/对照的标签当成批次传进去，结果 ComBat 把处理效应一并校正掉了。mod 参数里放的是你想保留的生物学因素。
  • 小心极端小批次：如果某个批次只有一个样本，ComBat 会直接跳过这个批次的位置校正（你会在控制台看到一条 warning）。原因是它没法估计该批次的方差，算不出可靠的 shrunken estimates。建议合并最小批次到邻近批次。
  • RNA-seq 不要混入芯片数据：两者的分布假设差距太大。我曾经试过把 RNA-seq 和芯片直接合并跑 ComBat，结果芯片数据的低表达基因全部被拉成了负值。后来我改用 limma::removeBatchEffect 配合 mean centering 才勉强解决。如果非混不可，先用 quantile normalization 对齐分布，再用 ComBat 只校 mean.only。
  • 避免循环使用 ComBat：有些朋友在大型队列里分批次多次校正——先校批次 A vs B，再校 B vs C。这会造成批次效应的「传递」，我建议一次性把所有批次标签传进去，让 ComBat 同时处理所有批次。
← 上一章 📖 返回目录 下一章 →