4. ComBat 方法原理与适用场景：经典批次校正方法的数学基础

ComBat 的全称是 Combining Batches，由 Johnson 等人在 2007 年提出。这东西在批次校正领域地位很高，说是「标配」一点都不夸张。

我第一次用 ComBat 是在一个多中心 RNA-seq 项目里。三个实验室的样本整合到一起，PCA 图一画，批次效应直接把生物学差异给淹没了。当时试了好几种方法，最后是 ComBat 把数据拉回来的。效果确实好，但我也踩过一些坑——比如没指定协变量，结果把真实的处理效应也校掉了。

这一节，咱们把 ComBat 的数学逻辑和适用场景拆开来讲。不堆公式，但核心思想得吃透。

4.1 从直觉到数学：ComBat 在做什么

说白了，ComBat 做的是两件事：调位置 和 调尺度。

你想想看，不同批次的数据，分布中心不一样（位置偏移），离散程度也不一样（尺度差异）。ComBat 对每个基因、每个批次，分别估计出这两个偏移量，然后把它们拉回到同一个基准上。

我个人习惯把这件事理解成「对齐 + 缩放」：

  • 位置对齐：减去该批次特有的均值偏移。
  • 尺度对齐：除以该批次特有的方差缩放因子。

但 ComBat 的高明之处在于——它不是硬生生地减均值、除方差。它用了 经验贝叶斯 框架，让估计更稳定，尤其是在样本量小的批次里。为什么这很重要？我后面会讲。

一句话直觉：ComBat 像是给每个批次的表达量做了一个「标准化」，但这个标准化是借了所有基因的力来做的，不是每个基因各自为战。

4.2 核心数学模型拆解

ComBat 的模型长这样（别怕，我一句句拆）：

Y_{ijg} = α_g + Xβ_g + γ_{ig} + δ_{ig} ε_{ijg}

我来解释每个符号是什么意思：

符号 含义 说明
Y_{ijg} 表达值 批次 i、样本 j、基因 g 的表达量
α_g 基因基线均值 该基因在所有样本中的平均表达水平
Xβ_g 生物学协变量效应 如处理组 vs 对照组，这部分我们要保留
γ_{ig} 批次位置效应 批次 i 中基因 g 的均值偏移
δ_{ig} 批次尺度效应 批次 i 中基因 g 的方差缩放
ε_{ijg} 随机误差 服从 N(0, σ²)，独立同分布

核心目标就是：估计出 γ_{ig} 和 δ_{ig}，然后从原始表达值里剥离它们。

调整后的数据长这样：

Y*_{ijg} = (Y_{ijg} - α̂_g - Xβ̂_g - γ*_{ig}) / δ*_{ig} + α̂_g + Xβ̂_g

注意看：分子把批次效应和基线都去掉了，除以尺度因子，再加回基线和协变量效应。这样做的目的是——只去掉批次效应，保留生物学变异。

经验贝叶斯是怎么进来的？

这里就是 ComBat 的精髓了。它给 γ_{ig} 和 δ_{ig} 设了 先验分布：

  • γ_{ig} 服从正态分布：γ_{ig} ~ N(μ_i, τ_i²)
  • δ_{ig}² 服从逆伽马分布：δ_{ig}² ~ InvGamma(λ_i, θ_i)

然后用所有基因的数据来估计这些先验参数（这就是「经验」二字的来源）。有了先验之后，再结合每个基因自身的观测数据，计算后验均值作为 γ*_{ig} 和 δ*_{ig}。

为什么这招管用？ 假设批次 1 只有 2 个样本，基因 A 的均值偏移估计会很不靠谱。但经验贝叶斯把所有基因的信息「借」过来，相当于给这个小批次做了正则化——估计结果更稳、更不容易过拟合。

我在项目中遇到过一个小批次只有 1 个样本的情况。如果不用经验贝叶斯，那个批次的位置效应根本没法独立估计。ComBat 硬是靠着「借力」把这组数据校了回来。效果当然不是完美，但至少能用了。

4.3 适用场景：什么时候该用它

ComBat 不是万能的。我总结了下面这个对照表，帮你快速判断：

适用场景 不适用 / 需谨慎的场景
批次数量 ≥ 2，样本分布较均衡 批次与生物学分组完全混杂（如批次1全是对照，批次2全是处理）
表达矩阵是连续值（log-CPM、FPKM、TPM 等） 输入是 counts 整数矩阵（建议先做 VST 或 log 转换）
批次效应是加性 + 乘性混合 批次效应仅表现为均值偏移（此时用 limma 的 removeBatchEffect 更轻量）
你关心的是差异表达基因的准确性 你关心的是单个样本的绝对表达值（ComBat 会改变数据分布）
有明确的生物学协变量需要保留 没有协变量信息，且批次与分组高度相关

我建议你在运行之前，先画一张 PCA 图看看批次和分组的关系。如果批次和分组完全重叠——嗯，那就不是统计方法能解决的问题了，得从实验设计层面找原因。

注意：ComBat 假设各批次共享相同的生物学效应。如果某个批次里只有一种处理条件，而另一个批次里只有另一种，它会把处理效应和批次效应弄混。这时候校正完，你的差异基因就全丢了。

4.4 关键参数与调参经验

R 包 sva 里的 ComBat() 函数是最常用的。核心参数我列一下：

library(sva)

# 基本用法
corrected_expr <- ComBat(
  dat = expr_matrix,        # 基因 × 样本 的表达矩阵
  batch = batch_info,        # 批次标签向量
  mod = design_matrix,       # 协变量设计矩阵（建议包含处理组信息）
  par.prior = TRUE,          # TRUE：参数先验（默认）；FALSE：非参数先验
  prior.plots = FALSE        # 是否绘制先验诊断图
)
参数 默认值 我的建议
par.prior TRUE 常用 TRUE，速度更快。如果批次数量少（2-3个），可以先 TRUE 后 FALSE 对比一下
mod NULL 一定要传！ 把处理组、性别、年龄等协变量放进去，防止被校掉
prior.plots FALSE 首次运行建议 TRUE，看看先验分布拟合得怎么样

关于 par.prior 的选择

我个人的经验是：

  • 样本量充足（每批 ≥ 10）：par.prior = TRUE 就够了，计算快，结果稳定。
  • 小样本或批次间差异极大：试试 par.prior = FALSE，非参数先验更灵活，但计算慢一些。
  • 两个结果差异很大：说明数据不够稳健，建议检查批次分配和协变量是否完整。

我曾经在一个只有 8 个样本、3 个批次的项目里，两种先验模式给出了完全不同的校正结果。后来发现是有一个样本被标错了批次。修好之后，两种模式就一致了。所以——先验诊断图真的很有用。

4.5 实战避坑指南

最后几个我在项目里踩过的坑，希望能帮你绕过去。

坑 1：忘了加协变量
如果你不传 mod 参数，ComBat 会默认所有样本属于同一总体。如果处理组和对照组在不同批次里分布不均衡，它会把处理效应当作批次效应校掉。我曾经因为这个差点把一个阳性结果给校没了。
坑 2：对 counts 直接跑 ComBat
ComBat 假设数据是连续且近似正态分布的。原始 counts 是离散的、偏态分布。我建议先做 log-CPM 或 VST (variance-stabilizing transformation)，再输入 ComBat。直接用 counts 的话，校正效果会打折扣。
小技巧：校正完之后，画一张 PCA 图，按批次着色。如果批次还聚在一起，说明没校干净；如果样本按处理组聚在一起、批次混在一起——那就是成功了。我习惯再用 distinct 包算一下批次间的 KL 散度或者层次聚类，量化一下改善幅度。

这就是 ComBat 的数学逻辑和使用要点。你想想看，它的核心其实就两件事：用经验贝叶斯稳健地估计批次偏移，然后把偏移从数据里剥离出去。理解了这个，后面的方法学起来就顺了。

ComBat 方法核心流程 模型估计阶段 贝叶斯校正阶段 多批次表达矩阵输入 估计 α̂_g + Xβ̂_g（基线+协变量） 估计 γ̂_{ig} + δ̂_{ig}（批次偏移） 设定先验分布（正态 + 逆伽马） 计算后验批次效应 γ*_{ig} + δ*_{ig} 调整表达值，去除批次效应 输出：校正后表达矩阵 关键要点 经验贝叶斯借助所有基因的信息，稳定估计小批次的偏移参数 必须传入协变量（处理组等），防止校正掉真实的生物学信号

我用下来最大的感受是：ComBat 很强，但前提是你得用对。协变量要传、数据要预处理、批次和分组不能完全混杂。把这几点守住了，它就能帮你把多个批次的数据整合成一个可靠的整体。

← 上一章 📖 返回目录 下一章 →