单细胞测序这几年火得一塌糊涂,但凡做肿瘤、发育、免疫相关的课题组,手里没几张UMAP图都不好意思跟人打招呼。但很多人跑完单样本分析后,卡在最关键的一步——多样本整合。我在实际项目中接过不少这样的活儿,发现大家对"整合"的理解普遍停留在"把数据合并在一起跑个harmony"的层面,对整合的本质逻辑、参数背后的意义、以及整合后如何验证效果,其实并不清楚。这篇文章把我这两年处理单细胞数据整合的经验完整梳理一遍,从为什么需要整合、数据准备阶段的潜在陷阱,到Harmony和Seurat CCA的底层逻辑差异,再到一套可以复现的完整流程和验证方法,希望帮你少走几个弯路。
1. 为什么单细胞数据不能直接合并跑分析
1.1 批次效应是什么,从一次失败的合并说起
先讲一个我早期踩过的坑。当时手里有两批PBMC数据,一批是自己实验室抽血后立刻上机的,另一批是合作的医院采集后冻存了三个月才做测序的。我天真地把两个矩阵直接cbind,然后RunPCA、RunUMAP,结果画出来的图完美分成了两团——我还以为发现了什么新的细胞亚群,高兴了好几天。后来用已知的marker基因一标,发现两团都是CD4 T细胞,只是分别来自两个批次。这就是最典型的批次效应:技术差异被误当成了生物学差异。
批次效应的来源很多,最常见的包括样本采集和处理时间差、建库试剂盒批号不同、测序深度不一致、冻存vs新鲜样本、甚至是同一个人不同时间操作带来的细微差别。单细胞测序对实验条件极度敏感,同样的细胞悬液分成两份,用不同批次的10x芯片跑出来,基因表达谱都可能产生系统性偏移。这种偏移在bulk RNA-seq里可以用limma或者ComBat之类的方法处理,但在单细胞水平上问题更复杂——因为每个细胞是一个独立的样本,我们不能像bulk那样对"同一个样本"做配对矫正。
如果直接合并跑分析,常见的后果有三个:聚类结果按照样本来源而不是细胞类型分开;同一细胞类型在不同批次中的marker基因表达被平均后失真;后续的差异表达分析会出现大量假阳性或假阴性。我见过有人因为不整合直接跑差异分析,最后找出来的"差异基因"其实全是批次相关的基因,比如线粒体基因、核糖体基因这类对细胞活性敏感的基因。
注意:线粒体基因占比高通常意味着细胞状态差,如果不同批次的细胞活性差异大,这个比例本身就会成为最强的"批次特征",直接合并时它甚至会主导PCA的前几个主成分。
1.2 什么时候必须整合,什么时候不需要
虽说整合是多样本分析的标配,但说实话,不是所有情况都适合无脑整合。我自己总结了一套判断标准,分享给大家参考。
第一种情况,单个样本内部有多个生物学重复(比如同一个处理组的三个重复),且主要目的是找组间差异基因,这种情况我建议保守一点:可以先用整合做细胞类型注释,但差异表达分析回到原始count矩阵上做,不要用整合后的校正值。原因后面会详细说,整合的矫正过程会压缩真实的生物学差异。
第二种情况,样本跨越了不同组织、不同发育时间点、或者不同物种,这种"大跨度的生物学差异"本身就足够强,通常不需要用整合算法去矫正,直接合并也能分得很开。强行整合反而可能把真实的组织特异性或发育阶段特异性信号抹掉。比如你拿肝脏和肺的数据去整合,harmony会试图让两批数据"看起来一样",结果就是某些组织特异性的marker被拉平了。
第三种情况,同一组织来源、不同个体、不同批次的样本,要合并做细胞类型图谱,这是最需要整合的场景。比如构建肿瘤微环境的细胞图谱,样本来自十几个病人,每个病人一个批次,如果不去除批次效应,聚类结果大概率是按病人分的,根本没法画出一张能代表整体的细胞类型图。
另外还有一点容易被忽视:如果你的数据只有两个样本,一个是处理组一个是对照组,且样本量很小,我强烈建议不要用复杂的整合算法。这时候批次效应和生物学差异几乎无法区分,任何整合方法都可能过度矫正。我碰到过不少这种情况,最后给的建议都是分别分析,然后做比较,而不是强行整合。
2. 整合前的数据准备:这一步决定整合质量的上限
2.1 数据质控和过滤的统一标准
很多人以为整合就是从读入数据开始的,其实整合的质量早在质控阶段就已经被决定了一大半。这里最关键的一点是:多样本质控标准必须统一,且这个标准要充分考虑不同样本的实际情况。
通常的质控指标有三个:每个细胞的UMI总数(nCount_RNA)、检测到的基因数(nFeature_RNA)、线粒体基因比例(percent.mt)。问题在于,不同样本的测序深度不一样,如果直接用一个全局阈值去卡,比如所有样本都卡nFeature大于1000,那么测序浅的样本会被卡掉大量细胞,导致样本间的细胞组成发生偏移。这种偏移在整合时会被算法当成批次信号,然后被矫正掉——但这不是技术批次,是人为质控造成的假象。
我的做法是分两步。第一步,用比较宽松的全局阈值做初步过滤,主要去掉明显的空液滴和破碎细胞,比如nFeature低于300的、percent.mt高于20%的。第二步,分别查看每个样本的分布曲线,针对个别样本适当调整阈值。比如某个样本整体测序深度偏低,nFeature的中位数只有800,那我就把它自己的下限设在500,而不是硬性卡1000。这样做的目的是保留每个样本的有效细胞,又不引入过多噪音。
这里再提一个很多教程不会讲但很实用的点:doublet过滤(双细胞检测)。整合分析对doublet的存在特别敏感,因为doublet本质上是一个"混合表达谱",它在PCA降维后往往落在两个真实细胞类型之间的过渡区域。整合算法会把这部分信号当成真实的中间态或者新的细胞亚群来处理,干扰锚点识别。建议在质控阶段用DoubletFinder或scDblFinder把每个样本的doublet预测出来并过滤掉,尤其在细胞上样量高的样本里,doublet比例可能达到5%到10%。
2.2 归一化和高变基因选择的细节
质控完成后,归一化是下一个关键节点。这一块常见的坑是表达式不对:整合分析需要在LogNormalize和SCTransform之间做选择,而且这个选择会影响后续的整合策略。
如果你打算用Harmony整合,那么标准的流程是LogNormalize(Seurat默认的NormalizeData),因为Harmony本身是在PCA嵌入空间里做矫正的,它对输入数据的分布要求没有SCTransform那么强。如果你用Seurat自带的方法(找锚点那套,后面细讲),两个都可以,但SCTransform一般效果更好、速度也更快,尤其适合高深度数据。
我一直用LogNormalize配Harmony,主要原因是我经常需要在整合后直接跑差异表达分析,而LogNormalize的结果更容易解释,下游工具的兼容性也更好。SCTransform的优点是对dropout和测序深度造成的噪音有更强的校正能力,但它计算量更大,而且在某些情况下会过度矫正,导致真实的生物学差异被压缩。
归一化之后是高变基因(HVG)选择。这一步很多人只是简单用FindVariableFeatures默认的2000个基因,这其实会埋下隐患。默认的vst方法选择的高变基因倾向于高表达基因,但如果你的样本里有比较大比例的特定细胞类型(比如肿瘤组织里有大量T细胞),高变基因就会被这些细胞类型的marker主导,那些在稀有细胞类型中起区分作用的低表达基因会被忽略。整合算法是基于高变基因做PCA的,如果高变基因没有代表性,整合的效果必然大打折扣。
我的习惯是把nfeatures从2000提到3000或4000,并且在确认每个样本的基因数足够的前提下使用SelectIntegrationFeatures来合并各样本的高变基因列表,而不是简单用某个单样本的结果。还有一个小技巧:用VariableFeatures检查一下选出来的基因里有没有包含免疫球蛋白基因(IGHG1、IGLC2这类)、线粒体核糖体基因(MRPS、MRPL系列),如果有,最好手动剔除。这些基因的表达水平和细胞类型、批次都有强关联,会强烈干扰降维结构,但它们本身不是我们关心的生物学信号。
注意:整合前的PCA也需要留个心眼。默认RunPCA会使用所有高变基因,但如果你发现PC1主要反映的是测序深度或者线粒体比例,这说明技术噪音没有被清洗干净,不建议直接用这个结果跑Harmony,而是应该回去检查质控和归一化参数。
3. 主流整合方法的底层逻辑:Harmony、CCA、scVI怎么选
3.1 Harmony:迭代聚类矫正的思路
Harmony是目前最流行的单细胞整合工具,没有之一。它的核心思想可以这样理解:先把所有细胞放到PCA降维后的低维空间里,然后反复做两件事——聚类和矫正。
具体来说,Harmony先对细胞做软聚类(soft clustering),每个细胞会被分配到若干个聚类中心,同时它会统计每个聚类中心里各批次样本的占比。如果一个聚类中心里的细胞几乎全来自批次A,说明这个聚类很可能是批次特异的,Harmony会计算一个矫正向量,把来自不同批次的细胞往一起拉;然后基于矫正后的坐标重新聚类,再重新计算差异,反复迭代直到收敛。这里面有两个重要参数:theta控制矫正力度,lambda控制每个聚类中心做矫正时允许的多样性程度。
Harmony的优势是速度快、内存占用低、对大规模数据友好,而且它的矫正是在PCA嵌入空间完成的,所以对下游的聚类和UMAP几乎是无缝衔接。在多数常规场景下,Harmony是性价比最高的选择。我跑过十几万细胞、二十多个样本的数据,Harmony计算只用了五分钟不到,这是其他方法很难做到的。
3.2 Seurat CCA:锚点对齐的思路
Seurat自带的整合方法(FindIntegrationAnchors + IntegrateData)思路完全不同,它用的是典型相关分析(CCA)加锚点对齐。CCA做的事情是找两组或多组数据共同的低维表示,使得不同数据集在同一维度上的相关性最大化。也就是说,它关注的是批次间共有的"结构",而不是批次内各自的特征。
在这个共享空间里,算法会寻找"锚点对"——即不同样本中彼此最相似的细胞对。这些锚点对就像两块布料需要缝合的位置标记,IntegrateData会根据这些锚点的信号,把不同样本的细胞往同一个"参考坐标系"里变换。
CCA方法的优点是保留了更多生物学信号,因为它锚定的是跨样本的相似细胞,对细胞类型组成差异很大的样本效果反而好。缺点是计算量大,几十万细胞的时候速度明显变慢,而且对内存要求很高。另外,如果样本间的细胞类型组成极度不平衡——比如一个样本全是T细胞,另一个全是上皮细胞——CCA可能找不到足够的锚点,整合效果就不好。
3.3 深度学习方法scVI
scVI是另一种思路的代表,它用变分自编码器(VAE)对单细胞表达数据建模。它的做法是把每个细胞的表达向量、批次的标签输入到一个神经网络里,让网络学习一个低维隐变量空间,在这个空间里批次效应被显式模型化为一个可分离的变量,从而得到批次矫正后的表示。
scVI在数据量足够大、测序深度差异很大的场景下效果非常出色,尤其对于UMI深度的建模和缺失值的处理做得比传统方法精细。但它的缺点也很明显:训练时间长,对GPU有依赖,模型参数较多,新手调参容易翻车。它对小数据集、几百个细胞的场景反而可能过拟合。
3.4 选型建议
我直接给一个经验性的选型表,纯个人经验,供参考:
| 场景 | 推荐方法 | 理由 |
|---|---|---|
| 10万+细胞、20+样本常规图谱 | Harmony | 速度快,稳定性好,默认首选 |
| 少样本(2-4个)、细胞组成差异大 | Seurat CCA | 锚点逻辑更适合小规模精细整合 |
| 测序深度差异极大(10x vs 其他平台) | scVI | 深度模型对测序深度的建模更精确 |
| 跨物种、跨组织整合 | Seurat CCA或scVI | Harmony在差异过大时可能过度矫正 |
| 只想要快速预览结果 | Harmony | 几行代码就能出图 |
不过我想强调一句,选型不是一道单向选择题。我在实际项目里经常是先跑Harmony看初步结构,如果发现某个区域整合不好,再针对那个亚群做二次整合,有时候会用Seurat的CCA方法再精调一次。方法之间不是互斥的,组合使用很常见。
4. 手把手完整流程:Harmony整合实战
4.1 环境准备与数据加载
这里我用Seurat + harmony包演示一套完整的整合流程。假设你已经有了多个10x样本的矩阵文件(标准格式是barcodes.tsv、features.tsv、matrix.mtx,或者你已经用Cell Ranger跑完得到了每个样本的filtered_feature_bc_matrix目录)。
library(Seurat) library(harmony) library(tidyverse) # 假设样本信息在samples向量中 samples <- c("sample1", "sample2", "sample3", "sample4") seurat_list <- list() for (s in samples) { seurat_list[[s]] <- Read10X(data.dir = paste0("data/", s, "/filtered_feature_bc_matrix/")) } # 创建Seurat对象,并给每个细胞加上样本标签 seurat_obj <- lapply(names(seurat_list), function(s) { CreateSeuratObject(counts = seurat_list[[s]], project = s, min.cells = 3, min.features = 200) }) %>% merge(y = .) # 在meta.data里加一列sample信息 seurat_obj$sample <- seurat_obj$orig.ident这里有个容易踩的坑:min.cells = 3这个参数的意思是基因至少在3个细胞中出现才保留。在整合场景中,多样本合并后细胞总数很大,min.cells可以适当提高,比如设为5或10,这样可以滤掉大量只在单细胞中出现的噪音基因,减少后续计算的负担。但注意不要设得太高,否则稀有细胞类型的marker会被删掉。
4.2 质控、归一化与PCA
# 计算线粒体基因比例 seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-") # 质控过滤 seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 500 & nFeature_RNA < 6000 & percent.mt < 20) # 归一化 seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000) # 高变基因 seurat_obj <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 3000) # 标准化 seurat_obj <- ScaleData(seurat_obj, vars.to.regress = "percent.mt") # PCA seurat_obj <- RunPCA(seurat_obj, npcs = 30, verbose = FALSE)标准化这一步,我要多说两句。ScaleData默认是对所有基因做z-score标准化,很多教程直接跑默认参数,但如果你不打算用SCTransform,我建议把vars.to.regress设置成percent.mt甚至nCount_RNA,把线粒体比例、测序深度这些技术变量回归掉。有些人对这种做法有争议,担心会误伤生物学信号,但在整合这个场景里,我的经验是利大于弊,尤其是样本间测序深度差异大的情况下。
关于PCA的维度选择,npcs = 30是经验值。理论上PC数量应该根据数据的复杂程度来定,复杂组织(比如大脑)可能需要更多PC才能捕获足够的信息。你可以用ElbowPlot辅助判断,但不要太迷信那个拐点。在整合场景中,我倾向于保守一点,把PC数量设得偏高(30到50),因为Harmony虽然是在PCA空间做矫正,但如果PC数量太少,有些稀有细胞类型的信息就没有进入矫正过程,后续聚类会丢失亚群。
4.3 Harmony运行与参数调整
PCA跑完以后就是Harmony的主角时间:
# 运行harmony,按sample变量矫正 seurat_obj <- RunHarmony(seurat_obj, group.by.vars = "sample", theta = 2, lambda = 1, reduction.save = "harmony")这里的group.by.vars指定了你要矫正的批次变量,最简单的场景就是样本ID。如果你的数据里还包含不同的实验批次、不同的测序平台,可以把它们都加进去,比如group.by.vars = c("sample", "platform"),Harmony会依次对每个变量做矫正。不过要注意,每多加一个矫正变量,过度矫正的风险就增加一分。我见过有人把无关紧要的变量也放进去矫正,结果整合完细胞类型都糊了。变量选择的标准只有一个——这个变量是否造成技术性系统差异,而不是生物学差异。
theta参数的默认值是2,它控制矫正力度:值越大矫正越激进,值越小越保守。如果你发现整合后还是有一些样本形成明显的独立簇,可以尝试把theta调高到3或4;如果发现过度矫正——比如不同细胞类型被强行拉到一起,把theta降到1试试。lambda参数一般保持默认1,它在聚类中心层面控制多样性的容忍度,我用过几次调整lambda的场景,但概率很低,大部分数据默认就好。
然后基于Harmony之后的嵌入表示做聚类和UMAP:
# 用harmony embedding做聚类 seurat_obj <- FindNeighbors(seurat_obj, reduction = "harmony", dims = 1:30) seurat_obj <- FindClusters(seurat_obj, resolution = 0.5) # UMAP seurat_obj <- RunUMAP(seurat_obj, reduction = "harmony", dims = 1:30)这里的关键是reduction = "harmony",一定不能漏。我见过不少人在跑完Harmony后,FindNeighbors和RunUMAP还是默认用pcareduction,等于整合了个寂寞。这条真的值得标红。
4.4 Seurat CCA流程对照
有些情况我会选择CCA,把流程也放出来供对照:
# 拆分对象 seurat_list <- SplitObject(seurat_obj, split.by = "sample") # 各自归一化 seurat_list <- lapply(seurat_list, NormalizeData, normalization.method = "LogNormalize") seurat_list <- lapply(seurat_list, FindVariableFeatures, nfeatures = 3000) # 找整合锚点 anchors <- FindIntegrationAnchors(object.list = seurat_list, dims = 1:30) # 执行整合 integrated <- IntegrateData(anchorset = anchors, dims = 1:30) # 整合后的数据需要重新标准化和PCA DefaultAssay(integrated) <- "integrated" integrated <- ScaleData(integrated) integrated <- RunPCA(integrated, npcs = 30) integrated <- RunUMAP(integrated, reduction = "pca", dims = 1:30)注意一个关键差异:CCA整合后的分析是在integratedassay上做的,而不是RNAassay。这意味着后续找marker基因、做差异表达、看基因表达量的时候,要用DefaultAssay切换回RNA,因为integratedassay里的数值是经过矫正的,不代表真实的表达量。这是个经典大坑,我身边因为这个问题分析出错误结果的人不在少数。
5. 整合效果怎么看:别只盯着UMAP
5.1 生物学保守性 vs 批次混合度
整合效果评估有两个维度:批次混合度和生物学保守性。这两个维度有时候是矛盾的,整合好的数据应该是在两者之间取得平衡。
批次混合度好理解——不同样本的细胞应该相互混合在一起,而不是各自抱团。UMAP上如果同一个样本的细胞聚成一大块孤立区域,基本说明整合没做好。但纯看UMAP容易误判,因为UMAP本身可以调整参数来改变视觉效果。我一般用两个定量指标来辅助判断:一个是计算每个细胞周围邻居中来自其他样本的比例,另一个是后面要讲的LISI分数。
生物学保守性是另一个容易被忽略的维度——整合不能把不同细胞类型之间的差异也抹掉。判断方法是:看已知的marker基因在整合后的UMAP上是否仍然清晰地区分对应的细胞类型。比如PBMC数据,CD3D应该只在T细胞区域表达,MS4A1只在B细胞区域。如果整合后这些marker变得模糊,甚至完全消失,多半是矫正过头了。
我经常用的一个快速验证方法:画出整合前后(PCA embedding vs Harmony embedding)的UMAP并排对比,同一组marker分别标两次。如果在PCA embedding中样本批次的主导地位明显,而在Harmony embedding中细胞类型的主导地位明显,同时marker没有被抹平,那这次的整合效果就基本达标了。
5.2 量化指标:LISI、kBET
定量评估有几个成熟工具,其中LISI(Local Inverse Simpson Index)是我最常用的。它衡量的是每个细胞局部邻域内的样本多样性:如果某个细胞的邻域里全是同一批次的细胞,LISI分数低;如果邻域里的细胞来自各个批次且比例均匀,LISI分数高。整合好的数据应该是高iLISI(批次混合度好)同时高cLISI(细胞类型保守度好)。Seurat有一个配套的lisi包可以算,代码不复杂,这里给个示意:
library(lisi) # 从Seurat对象提取embedding和meta信息 emb <- Embeddings(seurat_obj, reduction = "harmony") meta <- seurat_obj@meta.data # 计算LISI lisi_res <- computeLISI(emb, meta, batch_col = "sample") summary(lisi_res$sample)kBET(k-nearest neighbor batch effect test)是另一种经常被引用的方法,它的思路是检验每个细胞邻域内的批次构成是否与全局构成一致。简单说,它会对每个细胞做卡方检验,看看局部批次比例是否符合预期。结果会给出一个"接受率",接受率越高说明批次混合得越好。这两个指标相互补充,如果时间和算力允许,建议都算一下。
但我要泼一盆冷水:这些指标只是参考,不能替代生物学判断。我见过一些数据在LISI上分数很高,整合得很"完美",但仔细一看,某些稀有细胞类型被彻底抹掉了——因为把它们往混合方向矫正的时候,它们的特异性信号被当成了批次噪音。所以指标只是工具,最后还是要靠marker基因和生物学注释来把关。
5.3 标记基因验证
说到标记基因,这是整合效果评估里最硬核的一环。我的流程是:整合完成后,先不急着画画画写写写,先挑一组标志性的marker基因,用FeaturePlot标到UMAP上,逐一检查它们的表达模式是否符合预期。
这里有个容易被忽略的细节:如果用的是Harmony整合,FeaturePlot默认用的是RNAassay的数据来展示基因表达量,这个没问题;但如果用的是CCA整合,必须确保DefaultAssay(seurat_obj) <- "RNA",否则你画的是整合后的矫正表达值,可能跟实际情况对不上。另外,检查marker表达分布的时候,建议把pt.size调小一点(比如0.5),把order参数设为TRUE,让高表达细胞显示在上面,不然大量低表达细胞会把高表达信号遮住,看不出区分度。
除了看marker图,我还会顺便检查几个"应该没有"的基因。比如在纯T细胞数据里看看是否有B细胞的marker基因被异常激活,这种跨谱系的信号往往是整合过度矫正的警告信号。
6. 实战中踩过的坑与应对
6.1 细胞类型比例失衡导致的过度矫正
这是我在真实项目中遇到最多的问题。假设有10个样本,其中9个样本含有单核细胞,只有1个样本因为某次实验操作的原因,单核细胞几乎没有被捕获到。Harmony在矫正的时候发现单核细胞区域的聚类中心里有一个样本缺失,它不会认为这是生物学差异,而是倾向于把它当成"这批样本单核细胞没有被测到"的技术事故,于是它会把其他9个样本的单核细胞信号往这个缺失样本的方向拉,试图"填补"这个空白。结果就是单核细胞的真实信号被平均化,很多标记基因的表达被拉低,这个细胞类型在图谱上变得模糊。
这种情况的处理办法有两个方向:一是检查缺失的原因,如果确实是技术问题(比如样本在准备过程中丢失了某类细胞),可以在整合前就把这个样本排除或者单独处理;二是调整Harmony的lambda参数,让聚类中心允许更大的多样性,降低过度矫正的风险。但遇到这种情况我建议先把样本的细胞组成做一遍完整分析,确认缺失的原因再决定策略。
6.2 整合后marker基因"消失"
有次我在一个数据集上跑整合,发现CD8A这个基因在整个数据集中检测到的阳性细胞数量急剧下降,几乎"消失"了。这个问题的本质是:在整合矫正的过程中,低表达但细胞类型特异的基因信号被压缩了。尤其在使用PCA做降维时,那些方差贡献小的基因在保留的PC维度中占比很小,Harmony在矫正PC坐标时对这些基因的直接保护就弱。
应对方案有两个。一是在整合前确认高变基因选择是否充分,CD8A这类基因单独看表达量可能不高,但它的变异程度确实很高,应该会被选入高变基因列表,如果你发现没选进去,就手动加进去。二是整合后不要直接从矫正后的低维坐标去推断marker表达,特征基因表达应该回到原始count矩阵用DotPlot或VlnPlot来确认。我一般要求自己在写结论之前,所有关键的marker都要回到RNA assay里去验证一遍。
6.3 整合对差异表达分析的影响
这条建议适用于所有人:整合后的embedding可以用于聚类和可视化,但差异表达分析一定要回到原始count数据上做。原因在于,Harmony实际上是对细胞的低维坐标做了数学变换,它并没有生成一个新的"表达矩阵",任何基于矫正后embedding去做差异分析的做法,都是在拿降维后的信息丢失结果当输入,不靠谱。
CCA整合的情况风险更大,因为integratedassay里的值是经过变换的,已经偏离了原始计数分布。如果你拿integratedassay去跑FindMarkers,得到的差异基因可能只是矫正过程中引入的伪差异。正确做法是用DefaultAssay切回RNA,然后FindMarkers跑原始count数据。这时候你可能会发现一些整合后聚类定义的细胞群体之间差异不显著,这是正常的,因为整合时已经压缩了部分组间差异,你只要明确说明差异表达的输入是原始count就行。
6.4 内存和计算时间的控制
整合分析的计算资源消耗经常超出预期。我跑过最大的一个数据集是32万细胞、48个样本,Harmony用时不到十分钟,但之前的Standardize、PCA以及最后的UMAP才是真正的性能瓶颈。几个实用经验:
ScaleData是内存杀手,如果基因数太多(超过2万个),建议先用高变基因跑,或者开启SplitObject分块处理;- Seurat对象会同时保存RNA assay和integrated assay,内存开销翻倍,可以在确认整合效果后删掉不用的assay;
- UMAP的
n.neighbors参数默认30,对大数据集来说耗时会明显增加,可以适当降到15,虽然细节会有损失,但宏观结构几乎不变; - 纯CPU跑Harmony也是可行的,不需要GPU,不必为了它专门配一台带显卡的机器。
7. 写在最后的个人经验
单细胞数据整合这件事,表面上看是个技术活,本质上其实是"如何在区分技术和生物学信号之间找平衡"的判断活。我经常提醒自己,整合算法不是万能的魔法,它只是给了我们一个可以同时观察所有样本的坐标系,至于这个坐标系是否有生物学意义,最终还是需要用自己的领域知识去验证。
工具选型上,Harmony目前是我在80%场景下的首选,它确实稳定、快速、好调参,而且与该流程兼容性极好。但我不建议把Harmony当作唯一会用的方法,理解CCA锚点逻辑和scVI的建模思路,在遇到复杂场景时你会多几条退路。
如果只让我给一条最重要的建议,那就是——整合完成后,一定要做系统的、量化的效果验证,而不是只看一张好看的UMAP。把LISI分数、marker表达、样本混合结构、以及生物学合理性这四关都过了,你的整合结果才算真正可靠。另外多说一句,做整合之前,先把最终要回答的科学问题想清楚,你是想描述整体细胞图谱、还是研究某个亚群的状态变化?目的不同,整合策略和矫正力度都要跟着调整,这是比任何参数都重要的事。