4. 批次效应校正

好，咱们今天聊聊批次效应。这东西，说白了就是多组学数据分析里最让人头疼的问题之一。

你想想看，我做过的Pancancer项目，样本来自全球几十个实验室。不同机器、不同试剂批次、不同操作员，数据一合并，差异就来了。有些差异是生物学意义，有些纯粹是技术误差。

嗯，我今天就带你把ComBat、Limma、Harmony这三把刷子过一遍。我会讲清楚它们各自的原理、适用场景，还有我自己在实际项目中踩过的坑。

为什么需要批次效应校正？

先看一张图，帮大家理解批次效应到底是个什么东西。

看懂了这张图，你就理解了为什么必须做。不做校正，你今天发现的一个"标志物"，很可能只是不同批次之间的测量误差而已。

ComBat：贝叶斯框架下的稳健选择

原理要点

ComBat是我个人用得最多的方法，尤其在转录组数据中。它基于一个假设——批次效应是系统性偏移加上特征特异性的缩放。

数学上，它用经验贝叶斯来估计每个批次的均值和方差调整参数。说白了，就是"借力"——用所有特征的信息来稳定单个特征的估计。这在样本量小的批次里特别有效。

# R语言示例：sva包的ComBat
library(sva)
# 假设 expr 是基因表达矩阵（基因x样本）
# batch 是批次向量
combat_expr <- ComBat(
  dat = expr, 
  batch = batch, 
  mod = model.matrix(~condition),  # 保留生物学变量
  par.prior = TRUE,                # 使用参数先验
  prior.plots = FALSE
)

适用场景

• 转录组微阵列/RNA-seq：默认首选
• 小批次样本：经验贝叶斯能稳定估计
• 已知批次标签：必须有明确的批次分配

Limma：线性模型的优雅解法

原理要点

Limma本来是做差异表达分析的。但它内置的removeBatchEffect函数，也能做批次校正。

做法很直接：把批次作为协变量放进线性模型，拟合完以后，把批次贡献的部分扣掉，保留残差。这个方法好处是透明——你知道模型里包含了什么，排除了什么。

# R语言示例：limma包的removeBatchEffect
library(limma)
# 设计矩阵包含批次和生物学变量
design <- model.matrix(~ condition + batch)
# 拟合线性模型
fit <- lmFit(expr, design)
# 去除批次效应，保留condition的贡献
cleaned_expr <- removeBatchEffect(
  expr, 
  batch = batch, 
  design = design[, 1:2]  # 只保留生物学变量列
)

适用场景

• 基因表达数据：如果批次效应比较"干净"，线性假设就够了
• 甲基化芯片/蛋白组：特征维度没那么高时效果不错
• 需要保留完整自由度：ComBat可能损失部分方差，Limma更保守

Harmony：单细胞数据的利器

原理要点

Harmony是这几年火起来的方法，特别针对单细胞数据。它的思路和上面两个完全不同。

它不做"校正"后再分析。而是在降维的时候，迭代地调整数据，让来自不同批次的细胞在低维空间中"靠拢"。每次迭代，它会计算每个细胞的"模糊聚类分配"，然后修正主成分。

# Python + scanpy 示例
import scanpy as sc
import harmonypy

# 假设 adata 已做过PCA
adata.obsm['X_pca'] = pca_result  # 50维主成分

# 运行Harmony
ho = harmonypy.run_harmony(
    data_matrix = adata.obsm['X_pca'],
    meta_data = adata.obs,
    vars_use = ['batch'],
    max_iter_harmony = 10,
    do_pca = False
)

# 用Harmony校正后的嵌入
adata.obsm['X_harmony'] = ho.Z_corr.T

适用场景

• 单细胞转录组/scATAC-seq：数据稀疏、高度异质，Harmony很擅长
• 其他方法校正后仍有残留批次：可以尝试二次校正
• 可视化/聚类前处理：Harmony嵌入可以直接用于UMAP和聚类

方法横向对比

为了方便选择，我做了一张对比表：

实战中的选择策略

我通常会按以下步骤来决策：

1. 先看数据规模：样本数 < 50，优先ComBat；样本数 > 1000且是单细胞，优先Harmony
2. 批次数量：批次太多，Harmony优势明显；2-3个批次，ComBat和Limma很够用
3. 下游分析需求：需要完整表达矩阵做差异分析，不要用Harmony；只管聚类和可视化，Harmony很适合
4. 做对比试验：我的习惯是先分别跑一遍，用PCA/UMAP可视化看看批次的分离情况

# 一个简单的评估代码片段（R语言）
library(ggplot2)
library(uwot)

# 原始数据PCA
pca_raw <- prcomp(t(expr), rank. = 20)$x
# ComBat校正后PCA
pca_combat <- prcomp(t(combat_expr), rank. = 20)$x

# 绘制比较
df_plot <- rbind(
  data.frame(PC1 = pca_raw[,1], PC2 = pca_raw[,2], Method="Raw", Batch=batch),
  data.frame(PC1 = pca_combat[,1], PC2 = pca_combat[,2], Method="ComBat", Batch=batch)
)
ggplot(df_plot, aes(x=PC1, y=PC2, color=Batch)) +
  geom_point(size=1, alpha=0.7) +
  facet_wrap(~Method, scales="free") +
  theme_minimal()

避坑指南

我曾经在甲基化数据处理上吃过亏，分享一下我的教训：

• 别在训练集和测试集上分开校正：一定要先合并所有队列，再统一做批次校正。分开校正导致特征分布不一致，模型完全不可用。
• 校正前检查数据分布：如果某个批次全是NA或异常高表达，先做质控。我有一次没检查，一个批次有50%的缺失值，ComBat直接报错。
• 保留原始数据：校正后的数据不能替代原始数据。我的做法是，原始数据永不覆盖，校正后数据单独保存。
• 可视化验证：校正前后一定要看PCA/UMAP图。如果校正后不同批次的细胞完全混在一起，但已知的不同生物学分组也被抹平了，那就是过度校正。