做转录组项目只要样本一多、来源一杂批次效应就是绕不开的坎。明明是同一个处理第一轮建库的样本和第二轮建库的样本在PCA图上硬是分成两团明明是同一批细胞换了个测序机器之后差异基因列表看着就像换了物种。每次有师弟师妹拿这种图来问我我第一句话基本都是同一个先别急着跑差异分析你先把批次效应处理掉。去批次的方法很多但大家日常问得最多、也最容易搞混的就是三个ComBat、ComBat-seq、removeBatchEffect。这三个名字看着像一家子原理和适用场景其实差别很大。我这段时间刚好在整理一套多批次转录组数据的分析流程把这三个方法从原理到R代码再到使用场景完整过了一遍也踩了几个值得记录的坑。这篇就把我的整理结果写出来给正在跟批次效应搏斗的朋友做个参考。1. 批次效应来自哪里动手之前先做判断去批次不是拿到数据就直接跑函数第一步反而是在屏幕前把元数据看明白。批次效应的本质是技术层面的系统性偏差它和生物学差异混在一起让样本表达谱里既有真实信号又有技术噪音。如果连批次是怎么产生的、和数据里的哪些因素混在一起都说不清楚后面用再高级的方法也是一笔糊涂账。1.1 批次效应有哪些常见的源头先说最常见的几个来源。样本采集时间和地点不同是最容易出问题的地方我今天采的和上周采的光RNA降解程度就可能不一样。接着是建库环节试剂盒批次、反转录酶批次、加样操作人不同都会引入系统偏差。测序环节也一样同一批样本分了两次上机甚至同一个flow cell的不同lane都会带来强度差异。RNA-seq里还有个容易被忽略的因素文库构建的总量差异导致不同样本测序深度不同这个虽然不是严格意义上的批次但同样会造成类似批次的系统效应。我在实际项目里见过最典型的场景一个课题收了三个批次的临床样本第一批送了外送公司后两批换了一家最后拿回来的count矩阵合并在一起跑PCAPC1完美地把两个公司的样本分开。这种时候你不处理批次效应后续差异表达结果里一半的差异基因可能只是两个测序公司之间的技术噪音。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里正确使用ComBatComBat接收的输入是基因×样本的矩阵要求数值已经经过归一化和log变换比如log2(CPM1)或者log2(TPM1)。它的关键参数是batch和modbatch用来指定每个样本属于哪个批次mod用来放你想保留的生物学分组信息。library(sva) # 输入矩阵基因×样本log2(CPM1)转换后 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(datexpr_log, batchmeta$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数据不能直接套用ComBatRNA-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. removeBatchEffectlimma里的线性模型清理神器相比前面两个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 方法对比速查表方法输入数据类型统计原理输出形式主要适合场景主要局限ComBatlog2(CPM/TPM1)等连续矩阵高斯假设经验贝叶斯收缩连续表达矩阵芯片数据、微阵列、log尺度探索分析不适合count数据可能过度校正ComBat-seq原始整数count矩阵负二项回归离散度收缩保持整数特性的count矩阵RNA-seq差异表达前的批次校正速度慢基因数多时耗时较长removeBatchEffectlogCPM等连续矩阵线性模型残差化连续表达矩阵快速画图、聚类探索、机器学习特征无收缩对离群值敏感不能替代DE模型中的batch项选型的核心逻辑是看两件事你手上是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 常见错误与现场诊断现象可能原因排查与解决ComBat-seq报错提示count矩阵有非整数值误把CPM/RPKM或归一化后的数据传进去了回到最原始的featureCounts/HTSeq输出矩阵运行前用all(countsfloor(counts))检查去批次后PCA显示所有样本挤成一团没有在mod/group参数中传入生物学分组信息算法把分组差异也抹掉了检查是否传了mod或group补上分组变量重新校正校正后PCA还是明显按批次分开批次变量选错或者存在未被标记的隐藏批次检查元数据里是否有建库时间、试剂批次、操作人等字段尝试用batch2加第二个批次变量校正后差异基因数量暴增且marker基因方向改变过度校正可能是方法选择不适配换成ComBat-seq处理count或者在差异模型中直接加batch项不预校正ComBat结果出现NaN或极端值基因在所有样本中表达量过低或批次内样本表达全为常数先过滤低表达基因再去批次数据标准化时加上伪值如log2(CPM1)批次与分组完全混杂实验设计问题不是算法问题无法通过去批次解决尽可能补充交叉设计样本表格里列的场景我基本都亲自踩过其中不传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以及是否存在隐藏批次。我靠这个方法处理过好几个看起来莫名其妙的数据集最后都找到了问题所在。做去批次这件事方法本身只是最后一步搞清楚数据的来源和结构才是真正避免返工的关键。 转录组去批次效应实战：ComBat、ComBat-seq与removeBatchEffect选型指南
做转录组项目只要样本一多、来源一杂批次效应就是绕不开的坎。明明是同一个处理第一轮建库的样本和第二轮建库的样本在PCA图上硬是分成两团明明是同一批细胞换了个测序机器之后差异基因列表看着就像换了物种。每次有师弟师妹拿这种图来问我我第一句话基本都是同一个先别急着跑差异分析你先把批次效应处理掉。去批次的方法很多但大家日常问得最多、也最容易搞混的就是三个ComBat、ComBat-seq、removeBatchEffect。这三个名字看着像一家子原理和适用场景其实差别很大。我这段时间刚好在整理一套多批次转录组数据的分析流程把这三个方法从原理到R代码再到使用场景完整过了一遍也踩了几个值得记录的坑。这篇就把我的整理结果写出来给正在跟批次效应搏斗的朋友做个参考。1. 批次效应来自哪里动手之前先做判断去批次不是拿到数据就直接跑函数第一步反而是在屏幕前把元数据看明白。批次效应的本质是技术层面的系统性偏差它和生物学差异混在一起让样本表达谱里既有真实信号又有技术噪音。如果连批次是怎么产生的、和数据里的哪些因素混在一起都说不清楚后面用再高级的方法也是一笔糊涂账。1.1 批次效应有哪些常见的源头先说最常见的几个来源。样本采集时间和地点不同是最容易出问题的地方我今天采的和上周采的光RNA降解程度就可能不一样。接着是建库环节试剂盒批次、反转录酶批次、加样操作人不同都会引入系统偏差。测序环节也一样同一批样本分了两次上机甚至同一个flow cell的不同lane都会带来强度差异。RNA-seq里还有个容易被忽略的因素文库构建的总量差异导致不同样本测序深度不同这个虽然不是严格意义上的批次但同样会造成类似批次的系统效应。我在实际项目里见过最典型的场景一个课题收了三个批次的临床样本第一批送了外送公司后两批换了一家最后拿回来的count矩阵合并在一起跑PCAPC1完美地把两个公司的样本分开。这种时候你不处理批次效应后续差异表达结果里一半的差异基因可能只是两个测序公司之间的技术噪音。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里正确使用ComBatComBat接收的输入是基因×样本的矩阵要求数值已经经过归一化和log变换比如log2(CPM1)或者log2(TPM1)。它的关键参数是batch和modbatch用来指定每个样本属于哪个批次mod用来放你想保留的生物学分组信息。library(sva) # 输入矩阵基因×样本log2(CPM1)转换后 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(datexpr_log, batchmeta$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数据不能直接套用ComBatRNA-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. removeBatchEffectlimma里的线性模型清理神器相比前面两个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 方法对比速查表方法输入数据类型统计原理输出形式主要适合场景主要局限ComBatlog2(CPM/TPM1)等连续矩阵高斯假设经验贝叶斯收缩连续表达矩阵芯片数据、微阵列、log尺度探索分析不适合count数据可能过度校正ComBat-seq原始整数count矩阵负二项回归离散度收缩保持整数特性的count矩阵RNA-seq差异表达前的批次校正速度慢基因数多时耗时较长removeBatchEffectlogCPM等连续矩阵线性模型残差化连续表达矩阵快速画图、聚类探索、机器学习特征无收缩对离群值敏感不能替代DE模型中的batch项选型的核心逻辑是看两件事你手上是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 常见错误与现场诊断现象可能原因排查与解决ComBat-seq报错提示count矩阵有非整数值误把CPM/RPKM或归一化后的数据传进去了回到最原始的featureCounts/HTSeq输出矩阵运行前用all(countsfloor(counts))检查去批次后PCA显示所有样本挤成一团没有在mod/group参数中传入生物学分组信息算法把分组差异也抹掉了检查是否传了mod或group补上分组变量重新校正校正后PCA还是明显按批次分开批次变量选错或者存在未被标记的隐藏批次检查元数据里是否有建库时间、试剂批次、操作人等字段尝试用batch2加第二个批次变量校正后差异基因数量暴增且marker基因方向改变过度校正可能是方法选择不适配换成ComBat-seq处理count或者在差异模型中直接加batch项不预校正ComBat结果出现NaN或极端值基因在所有样本中表达量过低或批次内样本表达全为常数先过滤低表达基因再去批次数据标准化时加上伪值如log2(CPM1)批次与分组完全混杂实验设计问题不是算法问题无法通过去批次解决尽可能补充交叉设计样本表格里列的场景我基本都亲自踩过其中不传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以及是否存在隐藏批次。我靠这个方法处理过好几个看起来莫名其妙的数据集最后都找到了问题所在。做去批次这件事方法本身只是最后一步搞清楚数据的来源和结构才是真正避免返工的关键。