为何limma::removeBatchEffect()无法消除PCA中的批次聚类？如何优化？

在对RNA-seq数据进行方差稳定变换（VST）后，使用limma::removeBatchEffect()校正批次效应，但PCA结果仍显示明显的批次聚类。

极简示例代码：

library(DESeq2)
library(limma)

# Example VST matrix from DESeq2
vsd_mat <- assay(vsd)  # genes × samples matrix

# Metadata with batch and biological group (e.g. cancer subtype)
design_mat <- model.matrix(~ category, data = metadata)

# Attempt batch correction
corrected <- removeBatchEffect(vsd_mat,
                               batch = metadata$batch,
                               design = design_mat)

校正后运行PCA：

pca <- prcomp(t(corrected))
plot(pca$x[,1:2], col = metadata$batch)  # Still clusters by batch

请问该函数无法消除批次结构的原因是什么？该如何优化处理？

原因分析 #

• 批次与生物分组严重混杂：如果某批次内几乎全是单一生物分组样本（比如batch1全是A亚型，batch2全是B亚型），removeBatchEffect无法区分批次效应和真实生物差异，自然无法有效去除批次信号。
• 设计矩阵未适配校正逻辑：当前design_mat仅包含生物分组，但removeBatchEffect的design参数定义的是需要保留的效应。若批次与分组存在部分重叠，函数会避免过度校正以保护真实生物差异，导致批次信号残留。
• VST后仍残留技术噪音：VST虽能稳定方差，但低表达基因、极端离群值等仍可能携带与批次相关的技术噪音，这类噪音会主导PCA结果的批次聚类。
• 批次为全局强效应：若批次带来的表达偏移远大于生物差异，单纯线性校正不足以完全抵消这种系统性偏差。

优化处理方案 #

1. 先排查批次与分组的混杂情况 #

先统计各批次内的生物分组分布，确认是否存在严重混杂：

table(metadata$batch, metadata$category)

若存在严重混杂，优先剔除混杂批次/样本，或补充样本实现分组在批次间的平衡。

2. 调整removeBatchEffect参数设定 #

• 若批次与分组无严重混杂，可在设计矩阵中纳入批次，明确保留生物分组效应：

# 设计矩阵同时纳入批次和生物分组
design_mat <- model.matrix(~ batch + category, data = metadata)
# 校正时指定保留category效应，去除batch效应
corrected <- removeBatchEffect(vsd_mat,
                               design = design_mat,
                               covariates = metadata$batch)

• 若批次与分组存在交互作用，需在设计矩阵中加入交互项：

design_mat <- model.matrix(~ category + batch:category, data = metadata)
corrected <- removeBatchEffect(vsd_mat,
                               batch = metadata$batch,
                               design = design_mat)

3. 先过滤低质量基因再校正 #

过滤低表达、低变异的基因，减少技术噪音干扰：

# 过滤行均值低于1的基因（阈值可根据数据调整）
keep <- rowMeans(vsd_mat) > 1
vsd_filtered <- vsd_mat[keep, ]
# 基于过滤后的矩阵做批次校正
corrected <- removeBatchEffect(vsd_filtered,
                               batch = metadata$batch,
                               design = design_mat)

4. 更换更适配的批次校正工具 #

• 使用sva包的ComBat：专门针对批次-分组混杂场景优化：

library(sva)
mod <- model.matrix(~ category, data = metadata)
batch <- metadata$batch
corrected_combat <- ComBat(dat = vsd_mat, batch = batch, mod = mod)

• DESeq2内置批次校正：在构建DESeq2对象时直接将批次作为协变量，再生成校正后的VST矩阵：

dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = metadata,
                              design = ~ batch + category)
dds <- DESeq(dds)
# 生成已校正批次效应的VST矩阵
vsd_corrected <- vst(dds, blind = FALSE)
corrected_mat <- assay(vsd_corrected)

5. 调整PCA分析逻辑 #

若PCA批次聚类由少数高变异基因主导，可仅用高变异基因（HVGs）做PCA：

# 计算基因变异度，取前1000个高变异基因
vars <- rowVars(vsd_mat)
hvgs <- names(sort(vars, decreasing = TRUE)[1:1000])
# 基于高变异基因做PCA
pca <- prcomp(t(corrected[hvgs, ]))
plot(pca$x[,1:2], col = metadata$batch)

内容的提问来源于stack exchange，提问作者Timothy Adegbola