4. 经典方法：ComBat（参数贝叶斯）原理、假设与适用场景

好，聊完了批次效应的“长相”和“危害”，咱们得动真格的了。

这一节，我来聊聊批次校正领域的“老兵”——ComBat。全称是 Combating Batch Effects When Combining Batches of Gene Expression Microarray Data。

名字挺长，但核心就一句话：用贝叶斯统计，把不同批次的系统偏差“算”出来，然后减掉。

我在做多中心临床数据整合时，几乎把能试的方法都试了一遍。最后发现，ComBat 虽然“老”，但特别稳。很多场景下，它就是那个“保底”的选择。

一句话记住 ComBat：它假设每个批次都有自己的“系统偏差”（位置和尺度），然后通过经验贝叶斯把这些偏差估计出来，再对原始数据进行“对齐”。

4.1 ComBat 的核心思想：位置与尺度调整

说白了，ComBat 做了一件很聪明的事。它认为，批次效应主要影响两个东西：

  1. 位置（Location）：整个批次的均值偏移了。比如批次1的所有表达值整体偏高 0.5。
  2. 尺度（Scale）：整个批次的方差变了。比如批次2的数据分布得更“散”了。

那怎么办？

ComBat 为每个基因建立了一个模型。嗯，你没听错，是每个基因一个模型。它先假设数据长这样：

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

别被公式吓到。我们来拆开看：

  • Y_{ijg}：第 i 个批次、第 j 个样本、第 g 个基因的表达值。
  • α_g：这个基因的“本真”表达水平。就是没有批次效应时的基线。
  • Xβ_g：生物学变量效应。比如肿瘤 vs 正常，吃药前 vs 吃药后。这部分是我们要保留的。
  • γ_{ig}：批次 i 对基因 g 的“位置偏移”。这是 ComBat 要估计的主要参数。
  • δ_{ig}：批次 i 对基因 g 的“尺度缩放”。这是 ComBat 要估计的第二个参数。
  • ε_{ijg}：随机误差。假设服从标准正态分布。

你看，这个模型很“干净”。它把生物学差异（Xβ_g）和批次技术差异（γ_{ig} + δ_{ig}）彻底分开了。

4.2 经验贝叶斯：聪明的“借力”

这里有个关键问题：虽然我们假设每个基因都有自己独立的 γ 和 δ，但现实中，有些基因的样本量太少，根本算不准。

ComBat 的杀手锏来了——经验贝叶斯（Empirical Bayes）。

它怎么做呢？

  1. 先求“共性”：把所有基因的批次效应放在一起看，算出一个“先验分布”。比如，大多数基因的批次偏移 γ 都集中在 0 附近，只有少数会特别大。
  2. 再调“个性”：然后，对于那些样本量小的基因，ComBat 会说：“我看你的数据不太靠谱，所以我参考了兄弟们的经验（先验分布），把你的参数‘压’一下。”

这就是“收缩”的精髓。不会因为一个基因的数据波动，就把整个批次调得面目全非。

我在处理一个只有 3 个样本的小批次时深有体会。如果用普通的线性模型去逐个基因硬算，结果惨不忍睹。但 Combat 因为有收缩机制，愣是把那 3 个样本“救”了回来，校正后的分布跟大样本批次非常一致。

我的习惯：如果数据集中有样本量极小的批次（比如 n < 5），ComBat 几乎是唯一能稳定运行的方法。其他的方法要么报错，要么给出离谱的校正结果。

4.3 假设条件：你得知道它“吃”什么

任何模型都有前提假设。ComBat 也不例外。用之前，建议你对照检查一下。

假设条件 具体内容 我踩过的坑
1. 批次效应是加法/乘法模型 批次效应主要通过均值平移和方差缩放来体现。如果批次效应是非线性的，ComBat 搞不定。 在一次 RNA-seq 数据分析中，两个批次的测序深度差异导致计数分布极度扭曲（负二项 vs 正态），ComBat 直接“翻车”。
2. 生物学效应与批次效应可加 模型假设 Xβ_g 和批次效应是相加的，而不是相乘或交互。如果某个生物学效应只在特定批次出现，ComBat 会把它也当成批次效应抹掉。 这就是“过度校正”的常见原因。我曾经因为批次和分组完全混杂，把真实的治疗差异给校正没了，后来只能重做实验。
3. 各基因共享批次参数的先验分布 经验贝叶斯的根基就是“基因之间可以相互借力”。如果数据由完全异质的分子类型（如 RNA 和 DNA 甲基化数据）合并而成，这个假设就失效了。 我建议：分类校正。把 RNA 和蛋白数据分开用 ComBat，别混在一起。
4. 样本量足够估计先验 需要足够多的基因（通常 > 1000）来形成稳定的先验分布。 有一次我只测了 100 个靶向基因，ComBat 的收缩效果大打折扣，还不如直接做 Z-score 归一化。
警告：如果批次和分组变量完全共线（比如批次1全是对照组，批次2全是处理组），千万不要用 ComBat！它会把你关心的处理效应当成批次效应直接干掉。遇到这种情况，我建议先做实验设计去混杂，或者放弃分析。

4.4 适用场景：什么时候该请出 ComBat？

根据上面的假设，我总结了几个典型的适用场景。

  • 高维组学数据：基因表达芯片、RNA-seq 的 Count 或 FPKM 值、蛋白组学、代谢组学。这些数据特征多（成千上万），符合经验贝叶斯的“借力”前提。
  • 批次数量适中：2~10 个批次。批次太少（只有2个），先验分布估计不准；批次太多（>50），计算会变慢，但依然能用。
  • 生物学效应较强且明确：比如疾病 vs 健康、药物处理 vs 对照。ComBat 的假设是生物学效应远大于噪声，它擅长在强信号背景下“去噪”。
  • 需要保留原始数据的秩（Rank）：ComBat 是线性变换，校正后的数据在各批次内的秩不会改变，这一点对于某些排序类分析比较友好。

4.5 快速上手：Python 代码示例

废话少说，直接上代码。Python 的 combat 库封装得很好，核心就几行。

# 安装：pip install combat
from combat.pycombat import pycombat

# 假设你的数据长这样：
# data: pandas DataFrame，行是基因，列是样本
# batch: 一个 list 或 Series，长度等于样本数，内容为批次标签
# 注意：数据必须为数值型，不能有缺失值

# 标准化校正
corrected_data = pycombat(data, batch)

# 查看校正前后对比
print("校正前批次均值差异：",
      data.loc[:, batch == 'Batch1'].mean(axis=1).head(5))
print("校正后批次均值差异：",
      corrected_data.loc[:, batch == 'Batch1'].mean(axis=1).head(5))

嗯，就这么简单。但你得注意几个前置条件：

  • 数据格式：行是特征（基因/蛋白/代谢物），列是样本。很多新手把行列搞反了，直接报错。
  • 缺失值：ComBat 不支持缺失值。我习惯先用 KNN Imputer 或直接过滤掉缺失率 > 20% 的特征。
  • 数据分布：虽然 ComBat 本身对正态性不敏感，但如果是极度偏态的数据（如原始 Count 值），建议先做 log2 或 sqrt 变换。
我个人的代码习惯： 
# 预处理三步走
import numpy as np
import pandas as pd
from combat.pycombat import pycombat

# Step 1: 过滤低表达/高缺失特征
data = data[data.isnull().sum(axis=1) < 0.2 * data.shape[1]]
# Step 2: 对数变换（如果数据是 Count 或强度值）
data = np.log2(data + 1)
# Step 3: 运行 ComBat
corrected = pycombat(data, batch)

4.6 本章知识结构图

来，我画了一张图，帮你把 Combat 的“骨架”理清楚。它从输入到输出，中间经过了几个关键步骤。

ComBat批次校正方法核心流程 输入：多批次高维数据 步骤1：逐基因建立线性模型 Y = α + Xβ + γ(批次) + δ(批次) × ε 步骤2：估计批次参数(γ, δ) 计算基因级别的γ和δ，但依赖样本量 步骤3：经验贝叶斯收缩 对所有基因的γ和δ进行“共性强、个性稳”的收缩估计 步骤4：应用校正，输出数据 输出：与原始数据同维度的校正矩阵

图里梳理得很清楚了。从输入开始，经过 逐基因建模 → 参数估计 → 经验贝叶斯收缩 → 最终校正。整个 pipeline 的核心就是那个“收缩”步骤，让 ComBat 在数据量和稳定性之间找到了平衡。

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