单细胞转录组测序这几年从新鲜技术变成了很多课题组的常规武器但真正让做分析的人头疼的往往不是测序本身而是手上攒了七八个样本、横跨三四个批次之后怎么把它们揉成一个能用的对象。我见过太多人在这步翻车要么整合完细胞类型全混在一起要么过度整合把真实的生物学差异抹平了最后图是画出来了结论却站不住脚。这篇就围绕单细胞数据整合这件事把我在实际项目里踩过的坑、试过的方案、以及每个参数背后的逻辑尽量讲透。不管你是刚拿到一批数据的新手还是已经跑过几轮整合但结果总不太对劲的老手下面这些内容应该都能对上你的某个具体问题。1. 先搞清楚整合到底在解决什么问题1.1 批次效应的本质不是错误而是技术噪声很多人一上来就把批次效应当成bug恨不得除之而后快。这个理解是偏的。批次效应的来源非常具体不同时间点建库、不同操作人员、不同试剂批次、不同测序仪通道甚至同一个样本分两次上机都会引入系统性的表达偏移。它的本质是技术因素导致的、与生物学无关的表达差异。关键在于与生物学无关这个限定。如果两个批次之间本来就存在真实的生物学差异比如一个批次全是处理组、另一个批次全是对照组那你在整合时就必须格外小心因为算法分不清哪些差异该保留、哪些该消除。这就是为什么实验设计阶段要尽量让每个批次里都包含各种实验条件做到批次与生物学条件不共线。我见过一个项目四个批次里前两个全是对照组、后两个全是处理组整合完之后处理效应几乎被抹平这种设计上的硬伤后期任何算法都救不回来。所以整合的目标从来不是让所有细胞混在一起而是在保留真实生物学差异的前提下消除技术性的批次偏移。这个目标听起来简单但它是后面所有参数选择的判断标准。1.2 不整合会怎样三个典型症状不整合直接合并比如简单地把多个对象的count矩阵cbind起来会出现什么我总结下来有三个典型症状你可以对照自己的数据看看。第一个是聚类按样本分而不是按细胞类型分。UMAP图上每个样本自成一坨同一类细胞在不同样本里被拆成好几个cluster。这时候你去做差异分析得到的全是样本间差异不是细胞类型间差异。第二个是标记基因表达出现断层。比如某个T细胞标记基因在A样本的T细胞里高表达在B样本的T细胞里却几乎不表达但两个样本的T细胞明明是同一类。这种断层就是批次造成的。第三个是细胞类型比例失真。如果不同批次对某些细胞类型的捕获效率不同合并后比例会严重偏离真实情况下游的组成分析全部不可信。这三个症状里第一个最容易观察到也最容易被误判成我的样本异质性就是这么大。判断方法很简单看cluster和样本的对应关系如果高度重合基本就是批次问题。1.3 什么情况下可以不整合不是所有多样本数据都需要整合。如果你的所有样本是同一时间、同一操作者、同一批次建库测序的批次效应可能非常微弱这时候强行整合反而可能引入假象。判断标准是先做一次不整合的合并看UMAP上样本是否明显分离。如果样本间混合得很好细胞类型也能对上那就没必要整合。另一个例外是样本间生物学差异本身就是研究目标。比如你就是要比较疾病组和对照组且两组是不同批次测的这时候整合要非常保守甚至可以考虑不整合、改用其他策略比如分别分析再比较。整合是一把双刃剑用之前先想清楚你要保留什么。2. 主流整合方案的技术路线对比2.1 从CCA到RPCA再到Harmony三代方法的演进逻辑单细胞整合方法这几年迭代很快但主线其实很清晰。第一代以Seurat的CCA典型相关分析为代表思路是找到两个数据集之间表达模式最相关的维度然后在这些共享维度上做对齐。它的优点是能捕捉比较复杂的非线性关系缺点是计算量大、对细胞数敏感而且容易过度整合。第二代是RPCA鲁棒主成分分析Seurat v3之后主推。它不再找共享维度而是把每个数据集投影到一个公共空间同时保留各自的特异结构。相比CCARPCA在保留生物学差异上更稳计算也更快是目前Seurat工作流里的默认推荐。第三代是Harmony思路完全不同。它不重构表达矩阵而是在PCA降维之后迭代地调整每个细胞的低维坐标让批次间的同一类细胞靠拢。Harmony的优势是快、内存友好、对大数据集友好而且它直接在降维空间操作不改变原始表达值所以下游做差异分析时不会因为整合而失真。这也是为什么Harmony在很多场景下成了首选。三代方法没有绝对的优劣关键看你的数据规模、批次复杂度和研究目标。下面这张表是我根据实际使用经验整理的对比参数和表现都是常见实践下的典型值。方法核心思路适用规模计算速度过度整合风险典型场景CCA找共享相关维度对齐中小规模5万细胞慢较高批次差异大、样本少RPCA投影到公共空间保留特异结构中等规模中等中等Seurat默认工作流Harmony降维空间迭代校正坐标大规模10万细胞快较低多样本、多批次常规整合2.2 为什么Harmony在多数场景下更省心我自己的项目里Harmony用得最多原因很实际。第一它不碰原始表达矩阵这意味着你整合完之后做差异表达分析用的还是真实的count值不会因为整合算法编造了表达量而得到假阳性。第二它的输入就是PCA坐标你前面怎么降维、选多少主成分它都能接灵活。第三它对批次数量不敏感五六个批次和十几个批次都能处理不像CCA那样批次一多就崩。但Harmony也有它的边界。它假设批次效应在降维空间里是可以用线性方式校正的如果批次效应非常复杂、非线性很强Harmony可能校正不干净。另外Harmony的校正强度由参数控制调不好同样会过度整合。所以省心是相对的参数还是得调。2.3 选型决策三个问题帮你定方案面对一个具体项目怎么选我一般问自己三个问题。第一个问题细胞总数多少超过10万优先HarmonyCCA基本不用考虑RPCA也要看机器内存。5万以下三个都能试。第二个问题批次之间差异有多大如果UMAP上样本分离得非常厉害说明批次效应强可以考虑CCA或RPCA这种力度大的方法如果只是轻微分离Harmony足够。第三个问题下游要不要做严格的差异分析如果要优先Harmony因为它不改表达值。如果只是做细胞类型注释和比例分析三个都行。这三个问题问完方案基本就定了。实际项目里我经常是先用Harmony快速跑一版看整体结构如果发现某些批次校正不理想再回头试RPCA做对比。3. Harmony整合的完整实操链路3.1 整合前的数据准备质控和标准化不能省Harmony虽然是在降维空间操作但前面的质控、标准化、降维一步都不能马虎因为这些步骤直接决定了PCA坐标的质量而PCA坐标就是Harmony的输入。输入是垃圾输出必然是垃圾。质控阶段我通常关注三个指标每个细胞检测到的基因数nFeature_RNA、每个细胞的总count数nCount_RNA、以及线粒体基因占比percent.mt。阈值没有统一标准得看组织类型。比如PBMC这种细胞类型单一的样本nFeature在200到2500之间比较合理percent.mt控制在10%以下但如果是肿瘤组织或心肌这类线粒体本来就丰富的组织percent.mt的阈值可能要放宽到20%甚至更高。我一般会先画小提琴图看分布再定阈值而不是套用网上的固定数字。标准化用LogNormalizescale因子默认10000。这一步是把每个细胞的count除以自己的总count再乘10000然后取log目的是消除测序深度差异。做完标准化后找高变基因HVG通常选2000个。HVG的选择会影响后续PCA进而影响Harmony所以这一步也别偷懒可以用FindVariableFeatures的默认vst方法也可以根据数据特点调整selection.method。提示如果你的数据里包含不同物种或不同组织HVG最好分别找再取交集直接合并找HVG可能被某个占主导的样本带偏。3.2 多样本合并的正确姿势merge还是split这一步是很多人出错的地方。Seurat里合并多个样本有两种方式merge和split。merge是把多个对象合成一个大对象split是在一个大对象里按样本分层。做Harmony整合我推荐用merge因为Harmony需要所有细胞的PCA坐标在同一个矩阵里merge之后跑一次PCA就行。合并的时候有个细节样本信息一定要保留在metadata里。通常用add.cell.ids参数给每个样本的细胞加前缀避免细胞条码重复。合并后检查一下table(obj$orig.ident)确认每个样本的细胞数对得上。# 假设有四个样本对象 s1 s2 s3 s4 obj - merge(s1, y c(s2, s3, s4), add.cell.ids c(S1, S2, S3, S4), project integrated_project) table(obj$orig.ident)合并后重新做标准化、找HVG、scale、PCA。注意PCA要在合并后的对象上做不是在单个样本上做。因为Harmony需要所有细胞共享同一套主成分。3.3 跑Harmony参数怎么设、结果怎么看Harmony的核心函数是RunHarmony输入是PCA降维后的对象指定批次变量。批次变量就是你的样本来源通常是orig.ident。library(harmony) obj - RunHarmony(obj, group.by.vars orig.ident, reduction pca, dims.use 1:30, assay.use RNA, max.iter.harmony 20)几个关键参数。dims.use选多少主成分这个和PCA阶段保持一致一般20到30。选太少信息不够选太多噪声进来。max.iter.harmony是迭代次数默认10我一般设20让它充分收敛但也不是越大越好太大可能过拟合。theta参数控制校正强度默认2批次差异大的时候可以调到3到4但调太高会过度整合。跑完之后用Harmony的降维结果做UMAP和聚类。这时候看UMAP如果样本混合得比较均匀、细胞类型能聚成合理的cluster就说明整合有效。如果还是按样本分离可能是dims.use太小或theta太低如果所有细胞糊成一团、细胞类型分不开那就是过度整合要降theta或减dims。obj - RunUMAP(obj, reduction harmony, dims 1:30) obj - FindNeighbors(obj, reduction harmony, dims 1:30) obj - FindClusters(obj, resolution 0.5) DimPlot(obj, reduction umap, group.by orig.ident) DimPlot(obj, reduction umap, label TRUE)3.4 整合效果的量化评估别只看UMAPUMAP好看不代表整合成功这是我最想强调的一点。UMAP是高度压缩的二维投影它可能掩盖很多问题。我一般会配合几个量化指标一起看。第一个是批次混合熵batch mixing entropy衡量每个cluster里不同批次细胞的混合程度越高越好。第二个是细胞类型一致性用已知标记基因检查每个cluster的细胞类型是否在批次间一致。第三个是kBET或LISI这类专门的批次评估指标LISI分iLISI批次混合和cLISI细胞类型分离理想情况是iLISI高、cLISI低。实际操作里如果不想装太多包至少做两件事一是看每个cluster的批次组成用table(obj$seurat_clusters, obj$orig.ident)如果某个cluster几乎全来自一个批次要警惕二是用几个经典标记基因画feature plot看表达是否在批次间连续。4. 过度整合比不整合更隐蔽的坑4.1 过度整合的三种表现过度整合是整合分析里最危险的问题因为它看起来很成功——UMAP上细胞混得特别均匀cluster边界清晰一切都很美好。但实际上真实的生物学差异被抹掉了。表现一本该分开的细胞类型被合并。比如初始T细胞和记忆T细胞它们本来有明确的转录差异过度整合后可能被强行拉到一起变成一个cluster。表现二处理组和对照组的差异消失。你明明做了处理整合后两组细胞完全混合差异基因一个都找不出来。这时候不是处理没效果是整合把效果消掉了。表现三标记基因表达被平均化。某个基因在A群细胞高表达、B群低表达过度整合后所有细胞表达都差不多feature plot上一片均匀。这三种表现里第二种最容易被误判。很多人看到处理组对照组混合第一反应是我的处理没效果其实可能是整合参数太激进。4.2 怎么判断是不是整合过头了判断过度整合我有一套自己的检查流程。首先在整合前的PCA或UMAP上看一遍记住哪些细胞类型是明显分开的。整合后再看如果原本分开的现在合并了就要怀疑。其次用已知的、可信的标记基因做验证。选几个你百分百确定应该差异表达的基因看整合后它们的表达模式是否还保留。如果全被抹平就是过度整合。第三对比不同参数下的结果。把theta从2调到1或者dims从30降到20重跑一遍看细胞类型是否重新分开。如果参数一降就分开说明原来的设置确实过头了。第四看iLISI和cLISI的平衡。如果iLISI很高但cLISI也很高意味着细胞类型也被混了就是过度整合。4.3 回调参数的具体操作确认过度整合后怎么调回来按影响程度从大到小依次调这几个。先降theta。这是校正强度从默认2降到1甚至0.5让校正力度减弱。这是最直接的手段。再减dims.use。主成分少了保留的生物学结构反而更粗不容易被过度校正。从30降到20试试。然后可以调max.iter.harmony。迭代次数少了校正不充分反而保留了更多原始结构。但这个手段比较粗暴一般不首选。最后如果怎么调都过度整合可能要考虑换方法比如改用RPCA或者干脆不整合、用其他策略处理批次。注意调参不是越低越好。theta太低、dims太少又会回到不整合的老问题。理想状态是在批次混合和生物学保留之间找平衡点这个点因数据而异没有万能值。5. 整合后的下游分析注意事项5.1 差异表达分析用哪个assay很关键整合完之后做差异表达assay的选择直接影响结果可靠性。如果你用的是Harmony它没有改变原始表达值所以差异分析直接用RNA assay的count或data就行不需要特殊处理。但如果你用的是Seurat的CCA或RPCA整合整合结果存在integratedassay里这个assay里的值是校正过的不建议直接拿来做差异分析因为校正过程可能引入假信号。正确做法是用整合结果做聚类和UMAP但差异分析回到原始RNA assay。Seurat官方也推荐这个流程。具体操作是在FindMarkers里指定assay RNA并且用slot data标准化后的值。markers - FindMarkers(obj, ident.1 0, ident.2 1, assay RNA, slot data, only.pos FALSE, min.pct 0.25, logfc.threshold 0.25)5.2 细胞类型注释整合后注释更稳还是更乱整合的一个好处是细胞类型注释可以跨样本共享。因为整合后同一类细胞聚在一起你只需要注释一次所有样本都适用。这比每个样本单独注释再手动对齐要高效得多。但整合也可能让注释变难。如果过度整合某些细胞类型被合并你可能会把一个cluster注释成单一类型实际上它包含两种。所以注释阶段要结合多个标记基因不要只看一个。我一般会准备一套marker panel每个cluster至少用三到五个基因交叉验证必要时做亚聚类。另外整合后的cluster里如果批次组成极不均匀比如某个cluster 90%来自一个样本要警惕这个cluster可能是批次特异的假cluster而不是真实的细胞类型。5.3 细胞比例分析整合会不会扭曲比例细胞比例分析是很多课题的重点但整合对比例的影响经常被忽略。整合本身不改变细胞数量所以理论上比例不变。但整合会改变聚类结果如果整合把某些细胞归到了不同的cluster比例就会变。更隐蔽的问题是如果不同批次对某类细胞的捕获效率不同整合后这类细胞的比例反映的是捕获到的比例不是真实比例。这个问题在整合前后都存在不是整合引入的但整合容易让人忽略它。做比例分析时我建议同时报告整合前后的比例并说明差异来源。如果差异很大要谨慎解读。6. 几个真实项目里的踩坑记录6.1 批次与条件共线一个救不回来的设计前面提过这个坑这里展开说。有个项目六个样本分两批测第一批三个全是对照第二批三个全是处理。整合的时候不管用什么方法处理效应都被大幅削弱。因为算法看到的是批次间差异它无法区分这是技术批次还是生物学处理。最后的补救是放弃整合改成在每个批次内部做配对分析再比较两批的结果。虽然统计效力下降但至少结论方向是对的。这个教训是实验设计阶段就要避免批次和条件共线每个批次里都要有各种条件。6.2 细胞数极度不平衡导致的淹没另一个项目一个样本有3万细胞另一个只有2千。整合后小样本的细胞几乎被大样本淹没在UMAP上完全贴在大样本的cluster里看不出自己的结构。这是因为Harmony的校正是在全局做的细胞数少的批次权重低。解决办法有两个一是下采样大样本让各样本细胞数接近二是用group.by.vars时给不同批次加权部分工具支持。我一般选下采样简单直接虽然损失一些细胞但整合结果更均衡。6.3 整合后marker基因消失的排查过程有次整合完一个原本很明确的marker基因在feature plot上几乎看不到了。排查过程是这样的先确认基因名没写错再看不整合的版本里这个基因表达是否正常——正常。然后怀疑过度整合把theta从2降到1重跑基因表达恢复了。最后定位到是theta设太高。这个排查链路的关键是保留一个不整合的版本作为对照。任何时候整合结果不对劲先和不整合的版本比能快速定位是整合的问题还是数据本身的问题。6.4 内存爆掉的几种情况和应对Harmony虽然省内存但数据量特别大的时候还是会爆。我遇到过几次总结下来几个原因一是细胞数超过20万PCA矩阵太大二是dims.use设太高比如设了50三是同时开了多个R会话。应对方法细胞数大就分批整合或下采样dims.use控制在30以内跑之前清一下内存gc()一下。如果还是不够考虑用Harmony的Python版本或者换用更省内存的工具。7. 写在最后的一点个人体会整合这件事工具和方法固然重要但更关键的是你对自己数据的理解。你得知道哪些差异是真实的、哪些是技术造成的才能判断整合结果对不对。我见过太多人把整合当成一个跑完就行的步骤参数全用默认结果出来也不验证最后结论出问题都不知道错在哪。我的习惯是每次整合都至少跑两个参数版本一个保守一个激进对比着看。保守版本保留更多生物学差异激进版本批次混合更好然后根据研究目标选。如果两个版本结论一致那结论就很稳如果不一致就要深挖原因。另外整合不是终点是起点。整合完的object要拿去做注释、差异分析、轨迹分析等等每一步都可能引入新问题。整合阶段埋下的雷后面会一个个爆出来。所以这一步多花点时间绝对值。最后分享一个实用小技巧把每次整合的关键参数和对应的UMAP图存下来建一个自己的参数-效果对照表。跑的项目多了你会慢慢摸出适合自己常用组织类型的参数范围下次直接从这个范围起步能省很多试错时间。这个表我现在攒了十几行每次新项目先查表效率高不少。
