批次效应识别与去除

跟完这篇，您能用 PCA 一眼看出数据里有没有批次效应，判断批次和分组是不是被混在了一起，并在三种做法之间选对：用 removeBatchEffect 或 ComBat 得到「去批次矩阵」只为看图，用 ComBat-seq 得到校正后的 counts，还是更稳的——把批次直接放进统计模型当协变量。

什么是批次效应

同一批样本，如果分几次建库、几台仪器、几个日期、几个操作者去测，非生物学因素也会在数据里留下系统性的差异，这就是批次效应。它的麻烦在于：会和您真正关心的分组差异（如给药 vs 对照）混在一起，让本不该显著的基因显著、该显著的被淹没。

处理批次效应分两步，顺序不能反：先识别、判断能不能救，再决定怎么处理。识别没做清楚就急着「去批次」，很容易把真实的生物学差异一起抹掉。

第一步：用 PCA 把批次看出来

PCA 把上万个基因压到两三个主成分上，让您在一张散点图里看样本怎么聚堆。方法很简单：同一张 PCA 图，分别按「批次」和按「分组」上色，对比着看。

一个关键前提：不要对原始 counts 直接跑 PCA。 counts 的方差随均值变大，高表达基因会主导结果。先做方差稳定变换（VST）或 log 变换，再进 PCA。

library(DESeq2)

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = colData,        # 需含 batch、condition 两列
  design    = ~ batch + condition
)

# 方差稳定变换；blind=TRUE 表示不看分组、纯做质控
vsd <- vst(dds, blind = TRUE)

# 同一份数据，两种上色方式对比看
plotPCA(vsd, intgroup = "batch")       # 按批次上色
plotPCA(vsd, intgroup = "condition")   # 按分组上色

（plotPCA 默认取变异最大的 500 个基因，可用 ntop 参数调整。）

怎么读这两张图：

PCA 上的样子	说明	该怎么办
按分组上色分得开，按批次上色混着	分组差异是主导，批次影响小	通常无需专门处理，或把 batch 放模型里保险
按批次上色分得开，按分组上色混着	批次效应盖过了生物学差异	需要处理批次（见下文）
两种上色都分得开、且方向一致	批次和分组可能被混淆了	先查混杂，别急着去批次

第三种情况最危险，下一节专门讲。

判断：批次和分组是不是混杂了

去批次之前，必须先确认一件事：每个批次里是不是各组样本都有？ 如果不是——比如「批次 1 全是对照、批次 2 全是给药」——那么批次和分组就**完全混杂（confounded）**了，此时任何算法都无法把「批次导致的差异」和「给药导致的差异」分开，因为它们在数据里根本是同一个东西。

用一张交叉表就能看清楚：

table(colData$batch, colData$condition)

对着结果表判断：

交叉表的样子	混杂程度	结论
每个批次里各组都有、数量接近	无混杂	可放心处理，各法皆可
每个批次里各组都有，但数量不均	部分混杂	可校正，但检出力会打折
某些批次只含某一组	完全混杂	无法区分，任何去批次都是自欺

完全混杂是设计层面的缺陷，不是分析能补救的。这时硬跑 ComBat 或往模型里塞 batch，会遇到「模型矩阵不满秩（not full rank）」的报错，或者悄悄把真实差异当成批次抹掉。诚实的做法是如实说明这一局限；补做实验时，务必让每个批次都覆盖到各个分组。

为可视化去批次：removeBatchEffect

确认不是完全混杂、且批次已知，就可以处理了。如果您的目标只是看图——画 PCA、热图、做无监督聚类——最轻便的工具是 limma 的 removeBatchEffect。它在 log 尺度的表达矩阵上，把批次带来的系统性偏移扣掉。

library(limma)

mat    <- assay(vsd)                              # VST 后的 log 尺度矩阵
design <- model.matrix(~ condition, colData)      # 要「保留」的生物学变量

mat_nobatch <- removeBatchEffect(
  mat,
  batch  = colData$batch,
  design = design            # 传它，才不会把分组差异也一起去掉
)

# 用去批次后的矩阵重画 PCA，检查批次是否被压平

这里有两条红线，请务必记住：

一定要传 design（要保留的生物学变量）。 不传的话，若批次和分组有部分重叠，函数会连真实的分组差异一起抹掉。
removeBatchEffect 的输出只能用来「看」，不能回喂差异分析。 得到的 mat_nobatch 拿去画图、聚类没问题；但不要把它当作 DESeq2、edgeR、limma 差异检验的输入。原因见后面「更稳的做法」一节。

热图怎么用这份去批次矩阵，见表达热图怎么画。

ComBat 与 ComBat-seq：经验贝叶斯校正

removeBatchEffect 之外，另一个广泛使用的工具是 sva 包里的 ComBat。它用经验贝叶斯的方法在批次间「借力」，在批次数较多、每批样本偏少时通常比简单相减更稳。ComBat 有两个版本，用错版本是最常见的坑，先分清：

函数	输入	输出	用在哪
`ComBat`	log 尺度的归一化数据（如芯片、logCPM、VST）	校正后的 log 矩阵	可视化、聚类
`ComBat_seq`	原始 counts（整数）	校正后的 counts	可回喂 DESeq2/edgeR

RNA-seq 计数数据，用 ComBat_seq。 它专为计数设计，输出仍是整数样的 counts，可以继续走标准的差异分析流程：

library(sva)

adjusted <- ComBat_seq(
  as.matrix(counts),
  batch = colData$batch,
  group = colData$condition   # 告诉它哪些差异要保留
)
# adjusted 是校正后的 counts，可继续喂给 DESeq2 / edgeR

已经是 log 尺度的连续数据（如芯片），用 ComBat：

mod <- model.matrix(~ condition, colData)   # 要保留的生物学变量

combat_log <- ComBat(
  dat   = logexpr,            # log 尺度的归一化表达矩阵
  batch = colData$batch,
  mod   = mod
)

和 removeBatchEffect 同理：group／mod 参数一定要传，否则会把生物学差异一并校正掉。另外，ComBat 的经验贝叶斯在每批样本太少、或只有一两个批次时并不稳，此时结果要谨慎看待。

更稳的做法：把批次放进模型当协变量

到这里要讲最关键的一点。前面几种是「先把数据改了，再拿改过的数据去分析」。但对差异表达分析来说，更被推荐、也更稳妥的做法是反过来——不去动数据，而是把批次作为协变量写进统计模型，让模型在估计分组效应时自动扣除批次的影响。

在 DESeq2 里，就是把 batch 加进设计公式，且关注变量放最后：

dds <- DESeqDataSetFromMatrix(
  counts, colData,
  design = ~ batch + condition       # 关注的 condition 放最后
)
dds <- DESeq(dds)

# 默认给的就是「校正批次后」的 condition 效应
res <- results(dds, contrast = c("condition", "treated", "control"))

limma／edgeR 同理，把批次放进 model matrix 即可：

design <- model.matrix(~ batch + condition, colData)
fit    <- lmFit(logCPM, design)
fit    <- eBayes(fit)
topTable(fit, coef = "conditiontreated")   # 已校正批次

为什么这样更稳？因为**「先去批次、再做检验」会做两次校正**：去批次那一步已经用掉了数据里的自由度，后续检验却把这份数据当成原始观测，算出来的 p 值会偏乐观（假阳性偏多）。而把批次放进模型，是在同一个统计框架里一次性把批次和分组一起估计，自由度、方差都算对了，这也是 DESeq2、edgeR、limma 官方文档一致推荐的做法。

一个常被问到的补充：如果批次未知或没记录，但您怀疑存在未测量的系统性变异，可以用 sva 包估计出「替代变量（surrogate variables）」，再把它们当协变量放进模型，思路和放已知 batch 是一样的。

一张决策表

把上面的取舍浓缩成一张表，按「您要拿结果干什么」查一行：

您的场景	推荐做法	关键点
已知批次，做差异表达	放进模型：`~ batch + condition`	首选；不预先改数据
已知批次，画 PCA／热图／聚类	`removeBatchEffect` 或 `ComBat` 得去批次矩阵	只用于看，别回喂检验
下游工具只吃 counts、不吃协变量	`ComBat_seq` 得校正后 counts	计数数据专用
批次未知／未记录	sva 估计替代变量，当协变量放进模型	思路同放已知 batch
批次与分组完全混杂	无解，如实说明设计缺陷	补实验时让每批覆盖各组

几个反复出现的误区，一并记住：

别对原始 counts 跑 PCA：先做 VST 或 log 变换。
去批次矩阵别回喂差异分析：那是双重校正，p 值偏乐观；要做差异分析就把批次放模型里。
ComBat／removeBatchEffect 别忘传保留变量：不传 mod／design，生物学差异会被一起抹掉。
完全混杂别硬去：任何算法都救不回来，只会制造假象。

本文用到的 PCA 判读、混杂判断，都是定性的常规做法、非死板规则：批次效应有多严重、要不要处理，要结合样本量、实验设计和这批数据的实际分离程度来看，没有一条能套所有情况的硬阈值。

在平台上跑一遍

不想在 R 里配环境、也不想为「该放模型还是该 ComBat」纠结，可以在百沐一下直接做：上传表达矩阵和含批次、分组的样本信息表，用一句话说清「按批次和分组各画一张 PCA」「差异分析时校正批次」，平台会零代码替您画出 PCA、提示是否存在混杂，并在差异分析里自动把批次放进模型。您要做的只是核对批次列有没有标对、混杂判断合不合理。

完整流程见组学分析智能体使用指南，或直接开始：https://app.baimuyixia.com