我做了四五年单细胞转录组分析,见过太多人在SCTransform和Harmony这两个环节翻车。明明数据质量挺好,跑完聚类却乱七八糟,最后发现不是数据问题,是参数设置埋了雷。这篇文章专门聊SCTransform和Harmony的参数设置,把每个参数背后到底在干什么、为什么这么设、设错了会出什么问题,一条条掰开讲清楚。
如果你正在处理多样本单细胞数据,被批次效应搞得头疼,或者发现SCTransform跑完后的结果和预想差很多,这篇文章就是给你写的。新手可以把里面的参数组合直接抄走,已经跑过流程的老人也能对照检查一下自己有没有踩进某些隐蔽的坑里。
1. 先搞清楚SCTransform到底帮你做了什么
很多教程会告诉你“用SCTransform替代LogNormalize”,但没讲清楚它和传统标准化到底差在哪里。不理解原理就去调参数,基本等于盲人摸象。这一节先把SCTransform的核心逻辑讲明白,后面聊参数才有意义。
1.1 为什么大家开始放弃log-normalize
传统Seurat流程里,NormalizeData用的是log1p(CPM/10000)这类操作,本质是先把UMI数转成相对丰度,再取对数压缩动态范围。这个方案的问题在于:它假设所有基因的方差只和表达量均值有关,但单细胞数据里还存在技术噪声(比如测序深度差异)和生物学噪声(比如细胞周期、线粒体应激),这些因素会让高表达基因的方差被系统性放大。
你去做ScaleData的时候,通常会顺手做一步z-score标准化,但z-score的均值和方差估计都会被技术因素污染。结果就是,你挑出来的高变基因可能不是因为生物学差异大,而是因为某些细胞测序深度特别深。传统流程里有人用“回归掉UMI数和线粒体比例”来补救,但本质上还是在用线性回归处理非线性关系,效果有限。
SCTransform换了个思路:它直接用负二项回归建模UMI计数,把测序深度、线粒体比例这些技术因素作为协变量放进模型里,然后取残差作为校正后的表达量。这一步同时完成了标准化、方差稳定化和技术因素校正,比“NormalizeData + ScaleData + regress out”这条老路干净得多。
1.2 SCTransform的“回归”和“标准化”是一次完成的
很多人不知道,SCTransform跑完以后,SCT这个assay里其实包含了三类数据:
- counts:和原始RNA assay的counts一样,但只保留了参与建模的基因。
- data:对数化后的表达量(log1p),用于差异表达分析和可视化。
- scale.data:模型残差标准化后的结果,主要用于PCA和聚类。
关键理解是:后续PCA用的是scale.data,也就是“去掉技术因素后的残差”。所以如果你在SCTransform之后又画了FeaturePlot去看某个基因的表达量,默认用的其实是data里的log1p值,而不是残差。这会导致一个经典困惑:为什么基因表达图看起来没问题,但聚类结果里这个基因却没有贡献?因为聚类看的是scale.data,FeaturePlot看的是data,两者不是一回事。
这个信息在处理“某个marker基因在umap上明明表达,但找差异基因时居然不显著”这类问题时特别有用。真实场景中,我见过有人因为这个差异去反复调整FindMarkers参数,最后才发现是自己对数据的理解偏了。
1.3 SCTransform核心参数逐个拆解
先给出一份我常用的SCTransform参数清单,后面逐个讲为什么这么设。
obj <- SCTransform( obj, assay = "RNA", new.assay.name = "SCT", vars.to.regress = c("percent.mt"), variable.features.n = 3000, return.only.var.genes = TRUE, conserve.memory = TRUE, seed.use = 42, verbose = FALSE )- assay:指定从哪个assay取原始counts,默认是RNA。如果你之前跑过cellranger、STARsolo等流程,需要确认导入Seurat后counts确实在RNA assay里,否则会报错或者拿到空值。
- new.assay.name:默认叫SCT。如果你想保留多套标准化结果做对比,可以改成不同的名字。这是我比较推荐的做法,尤其是想比较SCT和LogNormalize结果差异的时候。
- vars.to.regress:需要回归掉的协变量。最常见的做法是回归percent.mt(线粒体比例),因为这个指标通常代表细胞应激或破损程度。但这里有一个很大的坑:并不是说“能回归的都回归掉”就是好事。细胞周期分数、样本ID这类变量别随便放进去,理由稍后单独细说。
- variable.features.n:要保留的高变基因数量,默认3000。如果样本类型复杂度高(比如肿瘤组织包含多种微环境细胞),我会调到4000-5000。如果数据比较简单(纯细胞系、纯T细胞亚群),2000-3000就够了。基因太少会丢失稀有亚群的信号,太多则把噪声也当成高变特征,PCA前几个主成分可能被无关基因主导。
- return.only.var.genes:默认为TRUE,意思是SCTransform后SCT assay里的scale.data只保留高变基因。这能显著省内存,也不影响后续PCA,因为PCA默认就是在高变基因上做的。但如果你之后想用某个不在高变基因列表里的基因去做热图或者打分,就需要小心了——这时候你可能需要临时把return.only.var.genes设成FALSE,或者用RNA assay的数据去提取表达量。
- conserve.memory:强烈建议设为TRUE。SCTransform对内存的消耗非常大,尤其是在10万细胞以上的数据集里。这个参数开启后会让SCTransform分块处理,内存占用能下降不少。缺点是速度会慢一点,但比起内存溢出来说完全值得。
- seed.use:随机数种子,保证结果可复现。写论文或者做项目时,这一项必须设置,否则每次运行结果可能在局部细节上有差异。
- verbose:设为FALSE可以省去大量中间日志,不然控制台刷屏刷得很难受,而且SCTransform本身日志量巨大,容易刷掉前面的报错信息。
1.4 vars.to.regress最容易被误用的坑
vars.to.regress这个参数坑特别多,我单独拎出来讲。
首先说percent.mt。常规认知里,线粒体比例高代表细胞状态差或者细胞破碎,所以把它回归掉是合理的。但我遇到过一个案例:处理某个带有线粒体基因高表达特征的细胞类型(比如某些代谢活跃的细胞)时,强行回归percent.mt后,这个细胞类型的特征基因信号被大幅削弱,聚类时这个群体直接消失了。
所以一个稳健的做法是:先用默认设置(回归percent.mt)跑一遍,看看主要细胞类型能不能分出来。如果某个你预期的细胞类型消失了,再尝试不回归percent.mt,对比两个结果。千万不要机械执行“别人都回归percent.mt所以我也回归”。
再说细胞周期。很多人喜欢用CellCycleScoring算出S.score和G2M.score,然后放进vars.to.regress。这个操作对增殖能力差异很大的组织(比如肿瘤)有时是必要的,因为细胞周期信号太强会把其他生物学信号盖掉。但如果你研究的就是干/祖细胞的增殖状态,或者你的目标细胞类型恰好和细胞周期高度相关,那这一回归就把你最关心的信号给干掉了。
我的经验是:先别急着回归细胞周期。跑完聚类看结果,如果发现聚类主要被细胞周期相关基因分开,而不是按细胞类型分开,再回归不迟。这个“先跑后决策”的思路,在单细胞分析里适用性极广。
另外,千万不要把样本ID或者批次信息放进vars.to.regress。SCTransform是每个细胞独立建模的,不是用来做批次整合的。如果把样本ID放进去,SCTransform会尝试把样本间的表达差异全部当作技术因素去掉,这可能把真实的生物学差异也一并抹掉。批次整合请交给Harmony、CCA或者scVI去处理,各司其职,别混用。
2. Harmony不是无脑跑的:参数背后的统计学含义
Harmony是目前整合多样本/多批次数据最常用的工具之一。它的效果确实好,但也正因为效果好,很多人不管三七二十一直接拿默认参数一顿跑,结果跑出来的整合结果是“看起来都混在一起了”,但生物学差异也被磨平了。要避免这个问题,必须理解Harmony在数学上做了什么。
2.1 Harmony在单细胞流程里的位置
Harmony不是替代SCTransform的,它是在PCA之后运行的。标准流程是:
- SCTransform:完成标准化和技术因素校正。
- RunPCA:在SCT assay的scale.data上跑PCA,得到降维后的细胞嵌入。
- RunHarmony:以PCA的嵌入结果作为输入,根据group.by.vars指定的批次变量做整合。
- 后续的FindNeighbors、FindClusters、RunUMAP全部基于Harmony校正后的嵌入。
这里有个很多人搞不清楚的点:Harmony到底改了什么东西?它不是改表达矩阵,而是在PCA降维后的低维空间里,对细胞的坐标做迭代校正。它先在每个聚类簇里估计“这个簇里不同批次的代表细胞”,然后通过混合多个批次的代表细胞来重新放置每个细胞的位置,最终让来自不同批次的同类细胞在低维空间里重叠。
这意味着Harmony不会改变你的基因表达值,只改变细胞在低维空间里的相对位置。如果你后续要做差异表达分析,用的还是SCTransform校正后的数据,而不是Harmony的输出。搞清楚这一点,就不会纠结“Harmony跑完后表达矩阵在哪里”这种问题了。
2.2 group.by.vars怎么选,别把生物学变量当批次
group.by.vars是一个字符串向量,表示“你认为哪些列代表了技术批次”。最常见的选择是sample ID或者library batch。
这个参数是Harmony整合效果的核心。设对了,整合效果好;设错了,后果往往很严重。举几个反面例子:
第一,把生物学分组(比如疾病组vs对照组)放进group.by.vars。这样Harmony会强行把疾病组和对照组的同类细胞往一起拉,最终导致疾病和对照之间的差异被磨平,后续找差异基因基本找不到,或者找到的都是假阳性。这个错误出现的频率高得惊人,尤其是有代码习惯的人复制粘贴上一轮脚本时最容易犯。
第二,把太多变量塞进group.by.vars,比如同时放sample ID、library batch、date、plate号。Harmony会尝试在所有维度上整合,结果就是每个维度都被弱化,而且计算量剧增。正确的做法是只放“你确信代表技术批次”的1到2个变量。如果你不确定某个变量是不是技术批次,可以在整合前先画个不需要整合的UMAP,看细胞是否按这个变量分堆。如果确实分堆了,再考虑放进去。
第三,样本数量太少时用Harmony要谨慎。我见过只有两个样本的数据集,其中一个样本细胞数特别少,Harmony跑完以后基本等于把所有细胞硬拉到一起,聚类结果完全丧失生物学意义。这种情况下,先用BBKNN、merge后直接聚类等方式做预实验,或考虑直接放弃整合,改用严格的样本间差异分析策略。
2.3 theta、lambda、sigma这些参数影响什么
Harmony的默认参数在多数数据集上表现都还行,但你如果想追求更好的效果,或者遇到了某些特殊场景,就需要调整下面几个核心参数。
先看默认参数是什么:
obj <- RunHarmony( obj, group.by.vars = "sample", theta = 2, lambda = 1, sigma = 0.1, max.iter.cluster = 20, epsilon.cluster = 1e-4, epsilon.harmony = 1e-4, plot_convergence = FALSE, verbose = FALSE )- theta:多样性惩罚参数,默认2。theta越大,Harmony就越强烈地推动不同批次在同一个聚类里混合;theta越小,混合力度越弱。如果你发现整合不足(不同批次还是明显分离),可以适当调大theta,比如3或4。如果你发现整合过度(不同细胞类型被强行拉到一起),可以调小theta,比如1或0.5。
- lambda:岭回归正则化参数,默认1。它控制批次代表细胞估计的平滑程度。lambda越大,估计的批次代表细胞越平滑,对异常值的容忍度越高;lambda越小,代表细胞估计越容易受到少数细胞的影响。多数情况下默认值1就可以,不用动。只有在样本量极不平衡时才需要考虑调大,比如一个样本有5万细胞,另一个只有1000,适当增大lambda可以让小样本的代表细胞估计更稳定一些。
- sigma:高斯核的带宽参数,默认0.1。它影响细胞距离相似度的衰减速度。sigma越小,只有非常近的细胞才被认为是相似的;sigma越大,较远的细胞也会相互影响。大多数情况下0.1是一个稳健的默认值。如果细胞类型很多且离散度高,可以试0.2;如果细胞类型之间边界很模糊,可以试0.05。
- max.iter.cluster:聚类迭代的最大次数,默认20。如果你发现Harmony没有收敛,日志里一直显示还有迭代在跑,可以适当提高这个值。但如果数据量很大,这个值也不要设太高,否则运行时间会显著拉长。实际处理20万细胞时,我会把max.iter.cluster降到10,因为20次迭代的运行成本太高,而且通常在10次左右就已经收敛了。
- epsilon.cluster和epsilon.harmony:收敛阈值,默认都是1e-4。一般不用调。如果你跑大数据集时觉得太慢,可以放宽到1e-3,能显著减少迭代次数,但整合精度会略微下降。我的建议是先默认跑一遍,观察日志里的收敛情况,再决定要不要放宽。
- plot_convergence:设为TRUE可以画出收敛曲线。第一次跑某个数据集时,建议打开看一下,确认Harmony确实在第几轮迭代后收敛了,免得白白浪费算力。
2.4 维数选择的蝴蝶效应
RunHarmony的输入来自RunPCA得到的PC嵌入。默认情况下,RunPCA会计算50个PC,但RunHarmony默认只取前30个PC作为输入。这个“dims”的选择对整合效果影响很大,经常被忽略。
如果dims设得太小,比如只取5-10个PC,你等于只保留了表达变化最剧烈的前几个方向,这通常是批次效应主导的方向,而真正的生物学信号可能排在后面,直接就被截掉了。结果是Harmony整合出来的空间里,细胞类型分不开。
如果dims设得太大,比如取到50个,那你把大量噪声方向也放进来了,Harmony会花很多精力去整合噪声,整合效果反而变差,还增加运行时间。
我的经验是:先用ElbowPlot画出每个PC的方差占比曲线,找到曲线变平的拐点。一般取拐点前后浮动5个PC。如果拐点不太清楚,保守起见取20-30个PC。对于大多数10x平台产生的单细胞数据,20-30是比较舒服的范围。你也可以用两种维度设置跑一遍UMAP做对比,看哪个结果里细胞类型分得更清楚且批次混合更均匀,拿实际结果说话。
另外注意,RunPCA和RunHarmony时的dims必须保持一致。千万别在RunPCA时用了npcs=50,然后到RunHarmony时只取dims=10,这样浪费了PCA后半段的信息,还有可能在后续FindNeighbors时维度对不上。建议把这些都写在同一段代码里,用同一个变量控制,减少出错概率。
3. 完整实操流程与参数组合建议
理论部分讲完了,这里给出一套可以直接复制到RStudio跑的完整代码。我以10x Genomics的PBMC数据为例,数据格式是cellranger输出的三个文件:barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz。实际项目中,你只需要把数据路径和metadata列名替换成自己的即可。
3.1 环境准备和必要包检查
需要确保R版本不低于4.1,Seurat版本不低于4.0。seurat-disk和harmony这两个包缺一不可。Harmony的安装比较简单,但从GitHub安装时偶尔会遇到依赖包编译问题,建议优先用CRAN版本。
library(Seurat) library(harmony) library(tidyverse) library(RColorBrewer) # 检查版本 packageVersion("Seurat") packageVersion("harmony")如果你的数据是cellranger输出,建议用Read10X读取,然后创建Seurat对象时加上min.cells和min.features过滤条件。这里min.cells=3表示基因至少在3个细胞中表达才保留,min.features=200表示细胞至少检测到200个基因才保留。这两个过滤条件能去掉绝大部分空液滴和破碎细胞,但是别把过滤阈值设得太狠,否则稀有细胞类型可能被直接过滤掉。
data_dir <- "path/to/your/cellranger/outs/filtered_feature_bc_matrix" counts <- Read10X(data.dir = data_dir) obj <- CreateSeuratObject( counts = counts, project = "my_project", min.cells = 3, min.features = 500 )3.2 SCTransform+Harmony完整流程
下面这一段是我日常项目里用得最多的一套参数组合。注意看,我在读取数据后先计算了percent.mt,然后传给SCTransform的vars.to.regress。
# 添加质量指标 obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-") # 可视化QC,确认过滤阈值 VlnPlot(obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) plot1 <- FeatureScatter(obj, feature1 = "nCount_RNA", feature2 = "percent.mt") plot2 <- FeatureScatter(obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") plot1 + plot2 # 过滤细胞 obj <- subset(obj, subset = nFeature_RNA > 500 & nFeature_RNA < 5000 & percent.mt < 20)QC这块有个细节:percent.mt的阈值不要一上来就拍脑袋定成5%或者10%。先画VlnPlot看分布,再决定阈值。如果大多数细胞都集中在2%-5%,那阈值定在10%就很宽松;如果大多数细胞都有10%左右的线粒体比例,贸然过滤到5%会把一大半细胞删掉。要基于数据本身做判断,不要套模板。
# SCTransform obj <- SCTransform( obj, vars.to.regress = "percent.mt", variable.features.n = 3000, conserve.memory = TRUE, seed.use = 42, verbose = FALSE ) # PCA obj <- RunPCA(obj, npcs = 30, verbose = FALSE) # 确认主成分 ElbowPlot(obj, ndims = 30)如果你之前没有跑过标准化,SCTransform会自动完成NormalizeData和ScaleData的操作,所以不需要再用NormalizeData。但如果你的数据里还保留了RNA assay,后续做FeaturePlot时默认绘制的可能是RNA assay的log1p表达量,所以没关系,不用额外处理。
# Harmony整合 obj <- RunHarmony( obj, group.by.vars = "sample", theta = 2, lambda = 1, sigma = 0.1, max.iter.cluster = 20, epsilon.cluster = 1e-4, epsilon.harmony = 1e-4, plot_convergence = TRUE, verbose = TRUE )这里group.by.vars=“sample”是假设你的metadata里有一列叫sample,表示每个细胞的样本来源。如果你的批次信息在别的列,比如batch、library,就改成对应的列名。
# 聚类与UMAP obj <- FindNeighbors(obj, reduction = "harmony", dims = 1:30) obj <- FindClusters(obj, resolution = 0.8) obj <- RunUMAP(obj, reduction = "harmony", dims = 1:30, seed.use = 42) # 可视化 DimPlot(obj, reduction = "umap", group.by = "sample") + ggtitle("By Sample") DimPlot(obj, reduction = "umap", group.by = "seurat_clusters") + ggtitle("By Cluster")我特别提醒一下:FindNeighbors里的reduction必须指定成“harmony”,很多人漏掉这一步,导致后续聚类用的还是未经整合的PCA结果。如果你在日志里看到类似“Using PCA as reduction”的提示,那就是没写上reduction参数。这种情况尤其常见在只需要跑单样本数据的脚本上,因为单样本根本不需要Harmony,所以没问题,但一换成多样本整合,漏掉reduction参数,结果就完全不同了。
3.3 参数组合速查表
下面是我在不同场景下的参数组合建议,可以直接当作参考配置抄走。
| 使用场景 | SCTransform核心设置 | Harmony核心设置 | 注意事项 |
|---|---|---|---|
| 常规PBMC/血液样本 | vars.to.regress = "percent.mt", variable.features.n = 3000 | theta=2, lambda=1, dims=1:20 | 通用首选,稳定但缺乏针对性 |
| 肿瘤组织(高复杂度) | variable.features.n = 5000 | theta=3, dims=1:30 | 高变基因多一点,否则稀有亚群信号容易丢 |
| 脑组织/复杂组织 | vars.to.regress = c("percent.mt"), variable.features.n = 4000 | theta=2, lambda=1, dims=1:30 | 注意细胞类型注释难度 |
| 极小样本量(2-3个样本) | return.only.var.genes = TRUE | theta=1 | theta调低,过度整合风险更大 |
| 大样本量(>20万细胞) | conserve.memory = TRUE, variable.features.n = 3000 | max.iter.cluster = 10 | 保证服务器内存,注意运行时间 |
| 细胞周期影响明显 | vars.to.regress = c("percent.mt", "S.score" , "G2M.score") | theta=2 | 先试试不回归的区别,确认细胞周期确实是混淆因素再回归 |
3.4 可视化检查批次整合效果
跑完Harmony之后,第一件事不是急着去看聚类注释,而是检查整合效果到底好不好。两个图必看:一个按样本上色的UMAP,一个按聚类上色的UMAP。
按样本上色的UMAP用来检查批次混合程度。理想情况是不同样本的细胞在同一个细胞类型区域里均匀混合,而不是各自聚成独立的一堆。但这里有个容易误判的点:如果某种细胞类型只出现在某一个样本中,那这个区域的细胞就只能来自这个样本,看起来“没有混合”其实是正常的生物学现象。所以检查整合效果时,一定要结合细胞类型注释来看,不要只看“是否完全混合”。
按聚类上色的UMAP用来检查聚类分辨率是否合适。如果用了0.8的resolution但聚类数量还是很少(比如只有4-5个),说明数据复杂度低,或者之前的基因筛选/整合过度了。如果聚类数量特别多(比如30+个),可以考虑降低resolution到0.4-0.5,或者先用SingleR做初步注释再看是否需要细分。
我见过一个更直观的检查方法:对某一个已知marker基因做FeaturePlot,看这个marker是否只在预期细胞类型里表达。如果marker表达模式良好,说明整合没有破坏生物学信号;如果marker表达被打散了,说明Harmony的强度过大,需要调低theta。
4. 常见问题与排查技巧实录
这一节把我在实际项目中遇到的高频问题整理出来,每一个都是真实场景,排查思路和解决方案可以直接套用。
4.1 问题一:SCTransform之后找不到counts数据
这种情况最常见于使用旧版Seurat(4.x之前)或者没有正确设置new.assay.name。如果你在SCTransform后运行FeaturePlot(obj, features = "CD3D"),有可能会报错说找不到这个基因。这是因为SCT assay的scale.data只包含了高变基因(return.only.var.genes=TRUE)。非高变基因的CD3D并没有进入scale.data。但FeaturePlot默认读取的是data槽,而不是scale.data,所以如果报错,多半是你把DefaultAssay切换得不对。
排查思路:
DefaultAssay(obj) <- "SCT" FeaturePlot(obj, features = "CD3D")如果还报错,看下这个基因是否被过滤掉了。在小细胞量的数据集中,很多基因会因为min.cells过滤条件被删除。这时你可以临时把DefaultAssay切换成RNA再试,因为RNA assay保留了全部检测到的基因。
4.2 问题二:Harmony跑了但聚类还是按样本堆
我在一个六样本的PBMC数据集上遇到过这种情况:跑了Harmony,但UMAP上样本还是各占一块区域,没有混合。排查后发现两个原因:一个是在FindNeighbors时忘了指定reduction="harmony",导致聚类用的还是PCA嵌入;另一个是theta设成了默认的2,但样本间差异实在太大,混合力度不够。
解决办法是两步走:
- 检查FindNeighbors/FindClusters/RunUMAP是否都用了reduction="harmony”。
- 把theta调大,从2调到4,同时适当增大max.iter.cluster(比如从20调到50),让Harmony有更多迭代次数去混合。这里要注意,调大theta可能带来过度整合的风险,所以每次只调一个参数,观察一下UMAP变化再决定下一步。
有时候问题也不在Harmony本身,而在上游的PCA维度选择。如果dims=1:10,而主要批次效应恰好集中在前10个PC之外,那么Harmony根本没有把那个PC纳入进来,自然整合不到。这种情况的解决办法是增加dims,比如改成1:30,看看是否有改善。
4.3 问题三:回归掉线粒体比例后细胞类型消失了
之前我处理一个肿瘤浸润免疫细胞的数据集,在SCTransform的vars.to.regress里回归了percent.mt,结果发现一群高表达线粒体基因的巨噬细胞亚群完全消失了。后来仔细查了一下,这群细胞高表达mitochondrial相关的基因是它们的生物学特征,而不纯粹是技术噪声。强行回归掉percent.mt,等于把这群细胞的身份特征给抹掉了。
解决方案是:换一套参数重新跑SCTransform,不回归percent.mt,对比两套结果。如果主要细胞类型在两种设置下都能分出来,说明不回归也没关系。如果少了某个细胞类型,再用回归percent.mt的结果来做后续分析,但需要在论文里说明这个处理对特定细胞类型可能的影响。
还有一种更稳妥的做法:把percent.mt作为“已知混杂因素”在Harmony整合时处理,而不是在SCTransform里回归。Real world数据里,批次效应、样本间差异、线粒体比例经常交织在一起,你需要在不同环节用不同工具处理不同因素,不能指望一个参数解决所有问题。
4.4 问题四:数据量大后内存爆掉
SCTransform在大数据集上非常吃内存。我之前处理一个20万细胞的数据集,16G内存的服务器直接OOM。后来发现是return.only.var.genes和conserve.memory这两个参数的组合出了问题。
如果你把return.only.var.genes设为FALSE(比如后续想用非高变基因做分析),SCT assay里会保存全部基因的scale.data,内存占用会爆炸式增长。更优的策略是:SCTransform时保持return.only.var.genes=TRUE,等需要用非高变基因时再随时提取RNA assay的原始counts做分析。这样可以大幅降低内存占用。
另外,如果你不需要保留原始RNA assay的全部信息,可以在一开始读入数据后先用NormalizeData,或者直接删除掉RNA assay(但一般不建议删,因为后续找marker基因时经常要用到counts数据)。最稳妥的方案是在创建Seurat对象时,就只保留必要的细胞和基因,减少数据量。
4.5 避坑清单速查
把上面提到的问题和关键点整理成一个速查表,方便你保存对照。
| 检查项 | 检查内容 | 推荐做法 |
|---|---|---|
| SCTransform回归变量 | 是否把细胞周期、样本ID、批次信息放入了vars.to.regress | 只用percent.mt,必要时再用细胞周期,样本ID和批次留给Harmony |
| 高变基因数量 | 是否设置了合理的高变基因数量 | 标准数据3000,复杂组织4000-5000 |
| return.only.var.genes | 是否与你的后续分析冲突 | 保持默认TRUE,需要时再临时修改 |
| conserve.memory | 大数据集是否开启 | 数据集超过5万细胞时强烈建议开启 |
| dims选择 | PCA和Harmony的维度是否一致 | 用ElbowPlot选拐点,常用1:20或1:30 |
| reduction指定 | FindNeighbors/RunUMAP是否用了harmony | 检查日志确认没有用PCA |
| theta值 | 是否根据整合效果动态调整 | 混合不足调大,过度整合调小 |
| group.by.vars | 是否把生物学分组误当批次变量 | 只放技术批次变量,不放疾病组/对照组 |
| 收敛性检查 | 是否确认Harmony收敛 | 用plot_convergence=TRUE查看曲线 |
4.6 一个容易被忽略的细节:细胞注释和Harmony的先后顺序
最后说一个很多教程都不提前讲清楚的事情:如果你打算用SingleR或CellTypist做自动注释,注释结果是基于原始表达矩阵的,和Harmony整合后的低维嵌入没有直接关系。所以先用注释工具给细胞打上类型标签,再基于Harmony整合后的嵌入去看不同类型细胞的分布,是一种非常高效的工作流。这样不但可以检查整合效果,还能快速发现“某些细胞类型是否在不同样本间有偏移”。
反过来,如果你先跑聚类再注释,聚类结果完全取决于你设置的resolution和Harmony参数,如果这些参数设得不好,注释结果也会被带偏。先注释再聚类,至少能保证你对细胞类型的判定是独立于整合参数的,后续调参时也有一个稳定的参照系。我在实际项目中几乎总是先做一轮初步注释,再回头调整Harmony参数,这样能少走很多弯路。
5. 写在最后的个人体会
这套SCTransform+Harmony的组合,我前前后后跑了不下几百个数据集,最大的感受就是:没有一套参数能吃遍天下,但理解了参数背后的原理之后,你会知道在什么场景下该动哪个旋钮。调参不是玄学,它是在“保留真实生物学信号”和“去除技术批次效应”之间找平衡,而找平衡的依据,永远来自你对自己数据的理解和一次次比较实验。
我自己的习惯是:每一个新数据集,都会先让SCTransform和Harmony用默认参数跑一遍,把QC图、elbow图、整合前后的UMAP图全部存下来留档,然后再根据这些图决定要不要调整参数。很多时候,问题不是参数不够好,而是你根本没来得及看默认结果就急着往下走了。多花十分钟看看那些图,后面能省下好几个小时的Debug时间。
最后分享一个小技巧:给每次分析建立独立的输出目录,把不同参数组合跑出来的UMAP图和聚类统计保存成带标签的文件,比如umap_theta2_dims30.png、umap_theta4_dims30.png,这样你可以随时翻回去对比哪个参数组合在你的数据上表现最好,也能避免复现实验时搞混了版本。这一点听起来很平常,但我见过太多人栽在这上面了,宁可多做一步存档,也不要事后抓狂。