企业官网建设流程全解析

做转录组项目，只要样本一多、来源一杂，批次效应就是绕不开的坎。明明是同一个处理，第一轮建库的样本和第二轮建库的样本在PCA图上硬是分成两团；明明是同一批细胞，换了个测序机器之后，差异基因列表看着就像换了物种。每次有师弟师妹拿这种图来问我，我第一句话基本都是同一个：先别急着跑差异分析，你先把批次效应处理掉。

去批次的方法很多，但大家日常问得最多、也最容易搞混的就是三个：ComBat、ComBat-seq、removeBatchEffect。这三个名字看着像一家子，原理和适用场景其实差别很大。我这段时间刚好在整理一套多批次转录组数据的分析流程，把这三个方法从原理到R代码再到使用场景完整过了一遍，也踩了几个值得记录的坑。这篇就把我的整理结果写出来，给正在跟批次效应搏斗的朋友做个参考。

1. 批次效应来自哪里，动手之前先做判断

去批次不是拿到数据就直接跑函数，第一步反而是在屏幕前把元数据看明白。批次效应的本质是技术层面的系统性偏差，它和生物学差异混在一起，让样本表达谱里既有真实信号，又有技术噪音。如果连批次是怎么产生的、和数据里的哪些因素混在一起都说不清楚，后面用再高级的方法也是一笔糊涂账。

1.1 批次效应有哪些常见的源头

先说最常见的几个来源。样本采集时间和地点不同，是最容易出问题的地方，我今天采的和上周采的，光RNA降解程度就可能不一样。接着是建库环节，试剂盒批次、反转录酶批次、加样操作人不同，都会引入系统偏差。测序环节也一样，同一批样本分了两次上机，甚至同一个flow cell的不同lane，都会带来强度差异。RNA-seq里还有个容易被忽略的因素，文库构建的总量差异导致不同样本测序深度不同，这个虽然不是严格意义上的批次，但同样会造成类似批次的系统效应。

我在实际项目里见过最典型的场景：一个课题收了三个批次的临床样本，第一批送了外送公司，后两批换了一家，最后拿回来的count矩阵合并在一起跑PCA，PC1完美地把两个公司的样本分开。这种时候你不处理批次效应，后续差异表达结果里一半的"差异基因"可能只是两个测序公司之间的技术噪音。

1.2 用什么方式确认数据里确实有批次效应

不要上来就默认有批次，也不要默认没有，用几个预处理图快速判断。我一般看三样东西：PCA图、层次聚类树、以及基因表达量与已知批次变量的相关性。

PCA是最直观的。把样本按元数据里的批次标色，如果前两三个主成分里样本清晰地按批次聚团，那就很明确有批次效应。层次聚类也有用，尤其当批次和另一个技术变量（比如RNA完整性RIN值高低）绑定的时候，聚类树经常会先按RIN值劈成两支。第三个办法是看表达矩阵和批次变量的关联强度，简单做一轮基因层面的检验或者相关性，看看有多少基因的表达量与批次显著相关，这个数量明显偏高就是信号。

library(ggplot2) library(limma) library(edgeR) # counts: 基因×样本的原始count矩阵 # meta: 样本元数据，至少包含 batch 和 group 两列 dge <- DGEList(counts = counts) dge <- calcNormFactors(dge) logcpm <- cpm(dge, log = TRUE) # PCA 快速检查 pca_res <- prcomp(t(logcpm), scale. = TRUE) pca_df <- data.frame(pca_res$x[, 1:2], meta) ggplot(pca_df, aes(x = PC1, y = PC2, color = batch)) + geom_point(size = 3) + theme_minimal()

如果这张图里不同batch的点各自抱团，那就说明批次效应已经大到影响全局结构的程度，需要认真处理了。

1.3 一个必须提前确认的问题：批次和分组是否混杂

这一步比选什么去批次方法都重要。批次效应处理的前提，是批次和生物学分组是交叉的——每个处理组里都同时存在多个批次，这样才能从数据里把批次和分组分开。如果批次和分组完全混杂，比如第一批都是对照组，第二批都是处理组，那任何算法都救不了你，因为在数学上无法区分到底表达差异来自处理还是来自批次。

遇到过类似设计的朋友应该都有体会，这时候无论跑ComBat还是ComBat-seq，出来的结果看着都对，但只要换个批次变量去标色，还是整整齐齐分成两组。唯一的解决办法是在实验设计阶段避免这种混杂，或者补充样本让批次和分组形成交叉设计。这个原则我在下面每个方法里都会反复提到，它是整个去批次流程的生死线。

2. ComBat：经验贝叶斯去批次的经典方案

ComBat是这批方法里资历最老的，来自Johnson等在2007年发表的论文，原本是给芯片数据设计的。十几年下来它成了批次效应校正的事实标准之一，R里面有现成的sva包可以直接调用。理解它的原理，对你弄明白后面两个方法为什么存在、各自解决什么问题，特别有帮助。

2.1 ComBat的原理：加性项、乘性项和收缩估计

ComBat的核心模型把基因表达量拆成这样几个部分。对于第g个基因、第j个样本（属于批次i），表达值被建模为：

Y_ijg = α_g + Xβ_g + γ_ig + δ_ig × ε_ijg

其中α_g是基因表达基线，Xβ_g是生物学协变量（也就是你想保留的处理分组等信息），γ_ig是批次i带来的加性偏移，δ_ig是批次i带来的乘性缩放。ε_ijg是随机误差项。

ComBat的聪明之处在于用经验贝叶斯来估计γ和δ。简单说，它先用数据估计出每一个基因在每个批次里的偏移和缩放量，然后做一个"向整体均值收缩"的处理。因为单个基因的估计会很飘，尤其当某个批次样本数很少的时候，单个基因的平均值可能完全不可靠。经验贝叶斯相当于把所有基因放在一起"投票"，得到一个先验分布，再把每个基因的具体估计拉回到先验均值附近。这样既保留了各基因自己的批次特征，又降低了少数离群基因对校正结果的干扰。

2.2 在R里正确使用ComBat

ComBat接收的输入是基因×样本的矩阵，要求数值已经经过归一化和log变换，比如log2(CPM+1)或者log2(TPM+1)。它的关键参数是batch和mod，batch用来指定每个样本属于哪个批次，mod用来放你想保留的生物学分组信息。

library(sva) # 输入矩阵：基因×样本，log2(CPM+1)转换后 expr_log <- log2(cpm(dge) + 1) # 模型矩阵必须包含你想要保留的生物学分组变量 mod <- model.matrix(~ group, data = meta) expr_combat <- ComBat( dat = expr_log, batch = meta$batch, mod = mod, prior.plots = FALSE )

这里最容易被忽略的是mod参数。很多教程为了省事不写这一项，直接ComBat(dat=expr_log, batch=meta$batch)。在分组和批次不相关的理想情况下，这样也能跑，但在分组和批次有部分相关性的时候，不提供mod会让ComBat把一部分生物学差异也当成分散性偏差给抹掉。我自己的原则是：只要数据里有分组信息，就一定要通过mod传进去，哪怕只是多写一行代码的事。

2.3 ComBat的局限：它不是为count数据设计的

ComBat在芯片数据上表现出色，但RNA-seq数据和芯片数据在统计分布上差别很大。RNA-seq的原始数据是整数count，方差与均值之间存在依赖关系，高表达基因的方差天然就比低表达基因大。我们习惯把count转成logCPM再用ComBat，本质上是用了一个高斯近似的假设。

在样本量充足、批次效应比较温和的场景下，这个近似工作得还算可以。可一旦某个批次样本量很少，或者文库大小差异悬殊，log变换加ComBat很容易出现两种问题：一是过度校正，把本来真实的基因表达差异也压缩了；二是校正后的矩阵是连续值，直接喂给edgeR、DESeq2这类工具时，count分布假设已经不再成立。正是这些坑，催生了专门给RNA-seq count数据设计的ComBat-seq。

3. ComBat-seq：专门喂给RNA-seq count数据的解法

ComBat-seq是Zhang等人在2020年发表在Bioinformatics上的方法，核心思路是把ComBat的经验贝叶斯框架搬到负二项分布上，直接处理原始整数count矩阵。如果你手里的数据是标准RNA-seq count，并且下游还要做差异表达分析，这个方法是目前最稳妥的选择之一。

3.1 为什么count数据不能直接套用ComBat

RNA-seq读段计数服从近似负二项分布，它的特点就是overdispersion——方差明显大于均值。你回想一下标准的泊松分布，方差等于均值，这在真实转录组数据里根本找不到。负二项分布能描述方差和均值之间的幂律关系，所以更贴近真实数据。

ComBat里的高斯假设把表达量当成均值附近对称波动，可count数据的分布是严重右偏的。你把低表达基因的count取log之后，大量0值和1值会聚成一堆，分布形态被扭曲了。这时候再做高斯假设下的经验贝叶斯校正，压缩的主要是那些本来就难以区分信号和噪音的低表达基因，结果就是校正矩阵里出现大量接近0的小数，看起来"均一化"了，实际上把count数据的离散程度抹掉了。

ComBat-seq解决的正是这个问题。它在负二项广义线性模型下估计每个基因的离散度，再对离散度做经验贝叶斯收缩，最后输出仍然保持整数特性的校正count矩阵。最大概率分布的关联，才不会被破坏。

3.2 ComBat-seq的原理和输入输出

ComBat-seq对每个基因拟合一个负二项回归模型：

log E[Y_g] = β0_g + β_group + γ_batch + log(N_j)

其中N_j是样本j的文库大小，也就是把所有样本的测序深度差异当成offset处理。模型里同时包含生物学分组项和批次项，然后通过经验贝叶斯把每个基因的离散度估计向整体水平收缩。校正时它会把批次项的贡献从预测值里去掉，保留分组项的贡献，再结合残差重新算出调整后的期望count。

这个设计带来的一个直接优势是：输出仍然是count。你可以把ComBat-seq校正后的矩阵直接拿去做edgeR或者DESeq2分析，而不需要在"已经从log尺度还原成伪count"这种绕弯操作里挣扎。另外一个优势是对稀疏基因更友好，因为它建模的时候就把低表达基因的零膨胀特性考虑进去了。

3.3 R实现与新版本sva包的注意事项

ComBat-seq的调用方式很简单，但参数细节决定了结果是否合理。

library(sva) # counts 必须是原始整数count矩阵（基因×样本） counts_adj <- ComBat_seq( counts = as.matrix(counts), batch = meta$batch, group = meta$group )

group参数的作用和ComBat里的mod类似，用来告诉算法哪些生物学分组信息需要保留。如果分析里有连续型协变量（比如年龄、临床指标），可以用covar_mod传入模型矩阵。需要特别注意的是group和covar_mod最好不要同时给，新版sva里如果两者都提供或者都不提供，函数会直接报错，要求你明确告诉它用什么结构来保留生物学信号。

另外，ComBat-seq要求输入必须是整数count。有次我图省事，直接把edgeR::cpm(dge)的结果传进去，运行到一半报错说检测到非整数值。后来养成了习惯，在调函数之前加一步类型检查：all(counts == floor(counts))，确认是整数矩阵再往下走。这个小检查帮我避开了好几次隐形错误。

4. removeBatchEffect：limma里的线性模型清理神器

相比前面两个，removeBatchEffect是最轻量、最快捷的一个。它来自limma包，本质上是线性模型框架下的残差化处理，用起来一行代码，特别适合在探索性分析阶段快速看图。但它的定位也清晰：它能帮你把数据"擦干净"去画PCA、做热图、跑聚类，但不能替代差异表达分析里把batch放进模型这一步。

4.1 稳定、快速的残差化思路

removeBatchEffect的原理不复杂。假设表达矩阵Y，我们希望保留生物学设计矩阵X对应的信号，去掉批次变量Z带来的影响。它先做一个线性拟合，估计出批次项的系数，然后从原始表达量里把批次项贡献的部分减掉，剩下的是"生物学信号+随机残差"。

这里的关键和ComBat不一样：removeBatchEffect没有经验贝叶斯收缩这一步，它直接使用最小二乘估计。这样做的优点是快、稳定、可解释；缺点是如果一个批次的样本里有极端离群值，这个离群值会直接影响该批次的均值估计，不像ComBat那样可以通过跨基因收缩来降低个别异常值的影响。

4.2 在R中实现并用于PCA等下游分析

removeBatchEffect的官方用法里有两个容易搞混的参数：design和batch。design放的是你希望保留的生物学变量，batch放的是你要去掉的批次变量。两个千万别搞反，不然你会把生物学差异抹掉、把批次差异留下，结果完全是反的。

library(limma) # logCPM矩阵 logcpm <- cpm(dge, log = TRUE) # 生物学设计矩阵，这里只放group design <- model.matrix(~ group, data = meta) # 去除batch效应，保留group差异 logcpm_adj <- removeBatchEffect( x = logcpm, batch = meta$batch, design = design ) # 用校正后的矩阵画PCA pca_adj <- prcomp(t(logcpm_adj), scale. = TRUE)

如果数据里有两个批次来源，比如建库批次和测序批次，可以把第二个传给batch2参数。还有连续型的系统来源，比如RNA的RIN值，可以用covariates参数传进去一并回归掉。我自己处理多中心数据时经常同时用到batch和covariates，比如把中心作为批次、把RIN值作为协变量。

4.3 它的定位：探索性分析工具，而非差异表达替代方案

必须说清楚：removeBatchEffect不应该是你做差异表达分析的主方法。limma的差异表达流程里面，推荐的做法是把batch作为协变量放进线性模型：

design_full <- model.matrix(~ batch + group, data = meta) fit <- lmFit(logcpm, design_full) fit <- eBayes(fit)

这才是标准的、在统计上更严谨的路线。removeBatchEffect适合的场景是那些没办法把批次变量塞进模型的下游分析——比如聚类、热图、PCA、机器学习特征筛选。在这些场景里，你需要的是一个"已经去掉批次噪音"的矩阵作为输入，所以通常先把removeBatchEffect的结果存下来，再往下走。

5. 三种方法到底怎么选，我常用的组合流程

写到这里，你大概已经看出三个方法各有所长。实际项目里很少只依赖一个方法走到底，更常见的是根据分析目标组合使用。下面先给一张对比表，再讲我日常跑的流程。

5.1 方法对比速查表

选型的核心逻辑是看两件事：你手上是count还是log值，以及你下一步要做什么。如果下一步是edgeR或DESeq2差异分析，源头又是count，那就走ComBat-seq。如果下一步是PCA、聚类、可视化，用logCPM加上removeBatchEffect就足够轻快。ComBat的位置比较微妙，它更适合那些你默认了高斯近似的场景，或者在旧流程里已经习惯用它并且结果验证过没问题的情况。

5.2 一个从counts到去批次验证的完整流程

我这套流程在多个转录组数据集上跑过，稳定性和可复现性都不错。第一步是检查批次与分组的交叉性，这一步在1.3里说过，不满足条件就直接停下，该补样本补样本，不要硬跑。第二步是做常规QC和过滤，去掉低表达基因，用edgeR::calcNormFactors做TMM归一化，同时保留原始count矩阵备用。第三步按目标分叉：要跑差异分析就用ComBat-seq处理原始count；要做无监督探索就用removeBatchEffect处理logCPM；如果课题组的旧流程里大家都用ComBat，那就明确记录参数，并且用PCA验证生物学分组没有被抹掉。

无论选了哪条路，最后都要做同一件事：验证。跑校正前后的PCA，对比两个图。理想的效果是批次分团消失，但生物学分组的分离度保住了。另外一个我常用的验证手段是挑几个课题组内公认的marker基因，看它们在两组之间的差异在校正前后是否保持方向一致。如果marker基因的差异方向都变了，说明校正过头了，需要换参数或者换方法。

# 快速输出校正前后的PCA对比 plot_pca_comparison <- function(counts_raw, counts_adj, meta) { logcpm_raw <- cpm(DGEList(counts_raw), log = TRUE) logcpm_adj <- cpm(DGEList(counts_adj), log = TRUE) p1 <- plot_pca(logcpm_raw, meta$batch, "Before correction") p2 <- plot_pca(logcpm_adj, meta$batch, "After correction") cowplot::plot_grid(p1, p2, ncol = 2) }

5.3 做机器学习或聚类时特别要注意的泄漏问题

这一条是很多临时转来做生信的人最容易忽略的：如果用去批次后的数据做机器学习分类，一定要警惕信息泄漏。假设你用全部样本做ComBat-seq校正，再用这个校正后的矩阵训练分类器、评估准确率，你的评估结果会偏乐观。因为校正过程里已经用到了测试样本的信息，测试集不再独立。

正确的做法是在交叉验证的每一折里面，只基于训练集样本估计批次校正参数，然后把参数应用到测试集样本上。ComBat和ComBat-seq本身都支持这种"用已有模型变换新数据"的做法。如果嫌麻烦不愿意在每一折里重跑，也要在论文里明确说明你是先整体校正再划分训练测试集的，并承认这是一个潜在局限。这个问题不影响肿瘤表达谱聚类的探索性结论，但会影响任何需要用准确率来支撑结论的场景。

6. 实操中踩过的坑和排查手册

方法说完了，最后整理一份我在实战里遇到的高频问题清单。这些问题在官方文档和教程里往往一句话带过，但真出问题的时候能卡住你好几天。

6.1 常见错误与现场诊断

表格里列的场景我基本都亲自踩过，其中"不传group导致生物学差异被抹掉"是最隐蔽的。它不会报错，PCA图看起来也变"干净"了，样本全部挤到原点附近，你甚至会觉得效果很好，直到发现分组之间的差异全没了才意识到出了问题。所以我的习惯是每次校正完不只画按批次标色的PCA，还会再画一张按分组标色的PCA，两张图对着看。

6.2 几个改善去批次效果的小习惯

第一个小习惯：去批次前统一基因ID格式和基因集范围。多批次合并数据经常出现同一个基因在不同批次的注释版本里用不同ID，比如Ensembl ID和Symbol混用。不统一就启动去批次，算法会把本来同一基因的几行当成不同基因处理，校正结果自然不可靠。

第二个小习惯：记录每次去批次的参数，包括软件版本。ComBat、ComBat-seq和sva包都还在迭代，不同版本之间的默认行为有差异，比如新版本里ComBat_seq对covar_mod的参数检查更严格。这一两年我吃过版本不一致的亏，同一个流程换台机器跑出来结果对不上，最后发现是sva版本不同。

第三个小习惯：把去批次这个环节当作分析流程的正式一步，写进报告。不要只在方法段落里写"we removed batch effects using ComBat-seq"就完事，还要写明输入矩阵是count还是log、是否传入分组协变量、用什么方法验证校正效果。审稿人或者合作方看到这些细节，对你的分析可信度会明显加分。

最后再分享一个我个人的工作习惯。每接一份新的转录组数据，我会先花半小时自己做一张"批次效应体检表"，把元数据里所有跟技术流程相关的字段全部列出来——采样日期、RNA提取批次、建库批次、测序批次、上机日期——然后用PCA逐一把每个字段标色，看看哪个维度能把样本分开。这个体检表花的时间不长，但能精准告诉你应该用哪个变量作为batch，以及是否存在隐藏批次。我靠这个方法处理过好几个看起来"莫名其妙"的数据集，最后都找到了问题所在。做去批次这件事，方法本身只是最后一步，搞清楚数据的来源和结构，才是真正避免返工的关键。