Seurat中sctransform标准化实战:原理、参数与避坑指南
2026/9/19 7:47:52 网站建设 项目流程

单细胞测序做到标准化这一步,基本就告别“跑通流程”的阶段了。我见过太多人把Seurat流程从头到尾走一遍,聚类图出来花花绿绿挺好看,结果一深究marker基因,发现全是核糖体基因和线粒体基因在作妖。问题十有八九出在标准化环节。sctransform这个方案在Seurat生态里已经不算新东西了,但真正把它用对、用透的人并不多。这篇文章就围绕Seuratsctransform标准化的实战展开,把原理、参数、踩坑点、和下游分析的衔接一次讲清楚。适合已经跑过基础Seurat流程、想进一步提升数据质量的单细胞分析人员,也适合正在纠结选LogNormalize还是SCTransform的同行参考。

1. 为什么标准化这一步值得单独拿出来讲

1.1 单细胞数据标准化的核心矛盾

单细胞转录组数据和传统bulk RNA-seq最大的区别在于稀疏性和技术噪音的耦合。一个细胞里某个基因表达量是0,可能是真的不表达,也可能是测序深度不够没捕获到。标准化要解决的核心问题是:让不同细胞之间的表达量可比,同时尽量不引入新的偏差。

LogNormalize的思路很直接——每个细胞的counts除以该细胞的总counts,再乘以一个缩放因子(默认10000),最后log1p转换。这个操作本质上假设所有细胞的总RNA含量是相似的,差异只来自测序深度。但实际情况是,不同细胞类型的总RNA含量差异可以很大,比如神经元细胞和淋巴细胞的总RNA含量就不在一个量级。这个假设不成立的时候,LogNormalize就会把生物学差异错误地归因到技术因素上。

sctransform的出发点完全不同。它用负二项分布广义线性模型来建模每个基因的counts,把测序深度作为模型的一个协变量。换句话说,它不假设所有细胞总RNA含量相同,而是通过模型来估计每个基因在给定测序深度下的期望表达量,然后对残差进行方差稳定化变换。这个思路的数学基础来自2019年Hafemeister和Satija发表在Genome Biology上的那篇论文,核心思想是让标准化后的数据在方差上更稳定,减少高表达基因对下游分析的过度影响。

1.2 sctransform解决了LogNormalize的哪些痛点

我整理了一个对比表格,把两种方法在实际使用中的差异列出来:

对比维度LogNormalizeSCTransform
核心假设所有细胞总RNA含量相似不依赖此假设,用模型估计
技术噪音处理简单缩放,不建模负二项GLM建模测序深度
方差稳定性高表达基因方差大方差稳定化变换
下游分析兼容性直接用于PCA、聚类需注意与旧流程的兼容
计算资源极低较高,尤其大样本
对稀有细胞类型的影响可能被淹没相对更敏感

从实际项目经验来看,sctransform在以下几种场景下优势明显:一是样本间测序深度差异大的时候,比如有的样本测了5000个细胞,有的测了20000个细胞;二是细胞类型异质性高的数据,比如脑组织、肿瘤微环境;三是需要做跨样本整合的时候,标准化后的残差数据在整合时表现更稳定。

sctransform也不是万能药。它的计算开销明显更大,尤其是细胞数超过5万的时候,内存和时间的消耗需要提前规划。另外,sctransform输出的Pearson残差有正有负,和LogNormalize输出的非负表达矩阵在数值范围上完全不同,下游做差异分析或者可视化的时候需要调整参数。

1.3 什么情况下该选SCTransform

我的判断标准比较直接:如果数据来自多个样本、多个批次,且样本间测序深度或细胞组成差异明显,优先考虑SCTransform。如果是单样本、测序深度均匀、细胞类型简单的数据,LogNormalize完全够用,没必要为了用而用。

还有一个容易被忽略的点:SCTransform线粒体基因和核糖体基因的处理。默认情况下,SCTransform会把所有基因都纳入模型,包括那些技术噪音很大的基因。如果不提前过滤,这些基因的残差可能会主导下游的PCA。所以我在实际流程里,通常会在SCTransform之前先做一轮基因过滤,把线粒体基因、核糖体基因、以及表达细胞数过少的基因去掉。

2. SCTransform标准化的核心参数与实操细节

2.1 核心参数逐个拆解

SCTransform函数本身参数不算多,但每个都值得说清楚:

seurat_obj <- SCTransform( seurat_obj, assay = "RNA", new.assay.name = "SCT", reference.SCT.model = NULL, do.correct.umi = TRUE, ncells = 5000, residual.features = NULL, variable.features.n = 3000, variable.features.rv.th = 1.3, vars.to.regress = NULL, do.scale = FALSE, do.center = TRUE, clip.range = sqrt(ncells), conserve.memory = FALSE, return.only.var.genes = TRUE, seed.use = 1448145, verbose = TRUE )

ncells:这个参数控制用于估计模型参数的细胞数。默认5000,意思是随机抽取5000个细胞来拟合负二项模型。如果你的数据细胞数少于5000,就用全部细胞。如果细胞数超过10万,可以适当调大这个值,但计算时间会线性增加。我一般会在细胞数超过2万的时候把ncells设成10000,实测下来模型更稳定。

variable.features.n:标准化后保留的高变基因数量,默认3000。这个数字不是越大越好。3000个高变基因对于大多数组织类型已经足够,如果做到5000以上,可能会引入一些噪音基因。我试过在PBMC数据上用5000个高变基因,结果发现前几个PC的贡献率反而下降了,说明多出来的基因没有提供有效信息。

vars.to.regress:这个参数在SCTransform里的行为和LogNormalize流程不同。在LogNormalize里,ScaleDatavars.to.regress是对标准化后的数据做线性回归取残差。而在SCTransform里,vars.to.regress是把这些变量作为额外的协变量加入负二项GLM模型。这意味着回归的效果更温和,不会像ScaleData那样把数据“洗”得太干净。我一般会回归percent.mt(线粒体比例),但不会回归nCount_RNA,因为SCTransform本身已经建模了测序深度。

clip.range:默认是sqrt(ncells),对于5000个细胞就是约70.7。这个参数控制Pearson残差的截断范围,防止极端值影响下游分析。如果你的数据里有特别极端的离群细胞,可以适当调小这个值,比如设成sqrt(ncells)/2

return.only.var.genes:默认TRUE,只返回高变基因的标准化数据。如果你后续要做基因集打分或者通路分析,需要用到全部基因的标准化数据,就要把这个参数设成FALSE。但要注意,返回全部基因会显著增加内存占用。

2.2 实操流程:从原始数据到标准化完成

我以PBMC数据为例,走一遍完整的流程。假设你已经有了一个初步过滤后的Seurat对象:

library(Seurat) library(ggplot2) # 假设seurat_obj已经完成了基础过滤 # 包括:nFeature_RNA > 200, nFeature_RNA < 2500, percent.mt < 5 # 第一步:先做一轮基因过滤 # 去掉线粒体基因、核糖体基因、以及表达细胞数少于3的基因 mt_genes <- grep(pattern = "^MT-", x = rownames(seurat_obj), value = TRUE) rb_genes <- grep(pattern = "^RP[SL]", x = rownames(seurat_obj), value = TRUE) genes_to_remove <- c(mt_genes, rb_genes) # 计算每个基因的表达细胞数 gene_cell_counts <- rowSums(seurat_obj[["RNA"]]@counts > 0) low_expr_genes <- names(gene_cell_counts[gene_cell_counts < 3]) # 合并要移除的基因 all_remove <- unique(c(genes_to_remove, low_expr_genes)) # 保留剩余基因 keep_genes <- setdiff(rownames(seurat_obj), all_remove) seurat_obj <- subset(seurat_obj, features = keep_genes) # 第二步:运行SCTransform seurat_obj <- SCTransform( seurat_obj, assay = "RNA", new.assay.name = "SCT", ncells = 5000, variable.features.n = 3000, vars.to.regress = "percent.mt", do.scale = FALSE, do.center = TRUE, clip.range = sqrt(5000), return.only.var.genes = TRUE, seed.use = 1448145, verbose = TRUE )

这里有几个细节值得展开说。基因过滤放在SCTransform之前,而不是之后。原因是SCTransform的模型拟合会受到极端基因的影响,如果线粒体基因占比很高,模型可能会把大量方差分配给这些基因。提前过滤可以让模型更专注于有生物学意义的基因。

do.scale = FALSEdo.center = TRUESCTransform默认不做scale,只做center。这是因为Pearson残差本身已经经过了方差稳定化变换,再做scale可能会过度校正。我试过把do.scale设成TRUE,结果发现聚类结果反而变差了,一些稀有细胞类型被合并到了大类里。

seed.use:这个参数控制随机数种子,因为ncells抽样是随机的。固定种子可以保证结果可重复。我一般用默认值,但在需要严格复现的项目里会显式指定。

2.3 标准化后的数据提取与使用

SCTransform跑完之后,标准化数据存在seurat_obj[["SCT"]]@scale.data里。注意,这里的scale.data存的是Pearson残差,不是传统意义上的scale后的数据。提取方式:

# 提取标准化后的数据矩阵 sct_data <- GetAssayData(seurat_obj, assay = "SCT", slot = "scale.data") # 查看维度 dim(sct_data) # 应该是 3000 x 细胞数(如果return.only.var.genes = TRUE) # 提取高变基因 variable_genes <- VariableFeatures(seurat_obj, assay = "SCT")

如果你需要全部基因的标准化数据,比如做基因集打分,有两个方案:一是SCTransform时设return.only.var.genes = FALSE,二是用GetResidual函数对特定基因集计算残差:

# 对特定基因集计算Pearson残差 genes_of_interest <- c("CD3D", "CD3E", "CD4", "CD8A", "CD8B") residuals <- GetResidual( seurat_obj, features = genes_of_interest, assay = "SCT", umi.assay = "RNA" )

GetResidual这个函数在实际项目中非常实用。比如你想看某个通路的所有基因在细胞间的变化,但不想重新跑一遍SCTransform,就可以用这个函数按需计算。

3. 标准化与下游分析的衔接:PCA、聚类、整合

3.1 PCA降维的注意事项

SCTransform跑完之后,下一步通常是RunPCA。这里有一个容易踩的坑:默认的RunPCA会使用scale.dataslot,而SCTransformscale.data是Pearson残差,有正有负。这本身没问题,但如果你之前习惯了LogNormalize流程,可能会对PC的解读产生困惑。

# 运行PCA seurat_obj <- RunPCA( seurat_obj, assay = "SCT", npcs = 50, features = VariableFeatures(seurat_obj, assay = "SCT"), verbose = TRUE ) # 查看PC贡献率 ElbowPlot(seurat_obj, ndims = 50)

我一般会看前30个PC的累积贡献率,如果前20个PC的累积贡献率不到80%,说明高变基因的选择可能有问题,或者数据本身的异质性不够。在PBMC数据上,前20个PC通常能解释85%以上的方差。

另一个细节是**npcs的设置**。默认是50,但对于大多数数据集,前30个PC已经足够。我试过在10万细胞的数据集上跑50个PC,发现从第25个PC开始贡献率就低于1%了,后面的PC基本是噪音。所以我的建议是:先跑50个PC看ElbowPlot,然后根据拐点选择后续聚类用的PC数。

3.2 聚类分辨率的选择

FindNeighborsFindClusters是聚类的核心步骤。SCTransform标准化后的数据在聚类时,分辨率参数resolution的敏感度和LogNormalize流程有所不同。因为Pearson残差的方差更稳定,聚类结果对resolution的变化相对不那么敏感。

# 构建邻居图 seurat_obj <- FindNeighbors( seurat_obj, reduction = "pca", dims = 1:30, k.param = 20, verbose = TRUE ) # 聚类 seurat_obj <- FindClusters( seurat_obj, resolution = 0.8, algorithm = 1, verbose = TRUE )

我一般会试多个分辨率,比如0.4、0.6、0.8、1.0、1.2,然后看聚类数的变化。在PBMC数据上,resolution = 0.8通常能得到8-12个cluster,对应主要的免疫细胞类型。如果resolution = 0.4resolution = 1.2的聚类数差异很大,说明数据的层次结构比较明显,可能需要考虑用clustree包来可视化不同分辨率下的聚类关系。

**k.param**这个参数控制KNN图的邻居数,默认20。对于细胞数超过5万的数据集,可以适当调大到30。但要注意,k.param太大会导致聚类过于平滑,稀有细胞类型可能被合并。

3.3 跨样本整合时的SCTransform策略

如果你有多个样本需要整合,SCTransform的使用方式有两种:先整合后标准化先标准化后整合。这两种策略的差异很大。

先标准化后整合:对每个样本单独跑SCTransform,然后用SelectIntegrationFeatures选择跨样本的高变基因,再用FindIntegrationAnchorsIntegrateData整合。这种策略的优点是每个样本的标准化是独立的,不会受到其他样本的影响。缺点是计算量大,尤其是样本数多的时候。

先整合后标准化:先把所有样本合并成一个Seurat对象,然后跑一次SCTransform。这种策略的优点是计算量小,缺点是如果样本间批次效应很强,标准化时可能会把批次效应错误地建模成生物学差异。

我一般推荐先标准化后整合,尤其是样本间差异明显的时候。具体流程:

# 假设有sample1和sample2两个Seurat对象 # 分别跑SCTransform sample1 <- SCTransform(sample1, vars.to.regress = "percent.mt", verbose = FALSE) sample2 <- SCTransform(sample2, vars.to.regress = "percent.mt", verbose = FALSE) # 选择整合特征 features <- SelectIntegrationFeatures( object.list = list(sample1, sample2), nfeatures = 3000, assay = c("SCT", "SCT") ) # 准备整合 sample1 <- RunPCA(sample1, assay = "SCT", features = features, verbose = FALSE) sample2 <- RunPCA(sample2, assay = "SCT", features = features, verbose = FALSE) # 找锚点 anchors <- FindIntegrationAnchors( object.list = list(sample1, sample2), anchor.features = features, assay = c("SCT", "SCT"), reduction = "rpca", dims = 1:30 ) # 整合 combined <- IntegrateData(anchorset = anchors, dims = 1:30)

这里有一个关键细节FindIntegrationAnchorsreduction参数设成"rpca"(Reciprocal PCA),这是Seurat v4之后推荐的方案,比默认的"cca"更快,而且对细胞数大的数据集更友好。我实测下来,rpca在10万细胞的数据集上比cca快3-5倍,整合效果差异不大。

整合完成后,combined[["integrated"]]@scale.data里存的是整合后的数据。后续的PCA、聚类、UMAP都基于这个assay。

4. 常见问题与排查技巧实录

4.1 SCTransform跑得特别慢怎么办

这是被问得最多的问题。SCTransform的计算复杂度大致是O(基因数 × 细胞数),但实际耗时受多个因素影响。我整理了一个排查表:

现象可能原因解决方案
跑了几小时没结束细胞数太多,ncells默认5000但基因数多减少基因数,先过滤低表达基因
内存爆了return.only.var.genes = FALSE设成TRUE,或分块处理
单样本快,多样本慢每个样本单独跑,总计算量大考虑先合并再标准化
特定基因导致模型不收敛极端高表达基因提前过滤或设clip.range

我的经验是:在跑SCTransform之前,先把基因数控制在20000以内。人类基因组大约有20000-25000个蛋白编码基因,加上非编码基因可能到30000以上。但单细胞数据里真正有表达的基因通常只有10000-15000个。提前过滤掉那些在少于3个细胞中表达的基因,可以显著减少计算量。

另一个技巧是future包做并行SCTransform本身不支持并行,但你可以用future框架把多个样本的标准化并行化:

library(future) plan("multisession", workers = 4) # 对每个样本并行跑SCTransform sample_list <- list(sample1, sample2, sample3, sample4) sample_list <- lapply(sample_list, function(x) { SCTransform(x, vars.to.regress = "percent.mt", verbose = FALSE) })

注意,并行化会显著增加内存占用,每个worker都需要独立的内存空间。如果内存不够,反而会拖慢速度。

4.2 标准化后聚类结果和预期不符

这个问题比较棘手,因为原因可能有很多。我一般按以下顺序排查:

第一步:检查高变基因。用VariableFeaturePlot看看高变基因的分布。如果前几个高变基因是线粒体基因或者核糖体基因,说明基因过滤没做好。

# 查看高变基因 top_genes <- head(VariableFeatures(seurat_obj, assay = "SCT"), 20) plot1 <- VariableFeaturePlot(seurat_obj, assay = "SCT") plot2 <- LabelPoints(plot = plot1, points = top_genes, repel = TRUE) plot2

第二步:检查PC的贡献率。如果前几个PC的贡献率异常高(比如PC1贡献率超过50%),说明数据里有强烈的批次效应或者技术噪音。

第三步:检查聚类分辨率。如果resolution = 0.8得到的cluster数远多于预期,可能是分辨率太高了。试着降到0.4或0.6。

第四步:检查vars.to.regress。如果回归了太多变量,可能会把生物学信号也回归掉。我一般只回归percent.mt,不回归nCount_RNAnFeature_RNA,因为SCTransform本身已经处理了测序深度。

4.3 Pearson残差有负值,下游分析怎么处理

这是SCTransform和LogNormalize流程最大的差异之一。Pearson残差有正有负,均值接近0。这带来两个问题:

问题一:差异分析。传统的差异分析工具(如FindMarkerswilcox检验)假设数据是非负的。用Pearson残差做FindMarkers时,需要把slot设成"data"而不是"scale.data"SCTransformdataslot存的是log1p校正后的counts,是非负的。

# 用data slot做差异分析 markers <- FindMarkers( seurat_obj, ident.1 = "CD4 T", ident.2 = "CD8 T", assay = "SCT", slot = "data", test.use = "wilcox" )

问题二:可视化FeaturePlotVlnPlot默认用dataslot,所以可视化不受影响。但如果你用DoHeatmap,它默认用scale.data,也就是Pearson残差。热图上的颜色范围会以0为中心,正负分明。这本身没问题,但解读时要注意,负值不代表不表达,只代表低于模型期望值。

我个人的习惯是:差异分析和可视化用dataslot,降维和聚类用scale.dataslot。这样既保证了统计检验的合理性,又利用了Pearson残差的方差稳定性。

4.4 和Harmony、LIGER等整合方法的兼容性

SCTransform标准化后的数据可以和其他整合方法配合使用。比如Harmony,它直接在PCA空间做整合,所以只要PCA是基于SCTransformscale.data跑的,Harmony就能正常工作:

library(harmony) # 先跑PCA seurat_obj <- RunPCA(seurat_obj, assay = "SCT", npcs = 50) # 用Harmony整合 seurat_obj <- RunHarmony( seurat_obj, group.by.vars = "sample", reduction = "pca", assay.use = "SCT", dims.use = 1:30 ) # 后续用harmony reduction做UMAP和聚类 seurat_obj <- RunUMAP(seurat_obj, reduction = "harmony", dims = 1:30) seurat_obj <- FindNeighbors(seurat_obj, reduction = "harmony", dims = 1:30) seurat_obj <- FindClusters(seurat_obj, resolution = 0.8)

Harmony的优势是速度快、内存占用低,尤其适合样本数多但每个样本细胞数不多的情况。我试过在20个样本、总共5万细胞的数据集上,Harmony整合只需要几分钟,而Seurat自带的FindIntegrationAnchors要跑半小时以上。

但Harmony也有局限:它假设批次效应是线性的,对于复杂的非线性批次效应,效果可能不如Seurat的anchor-based方法。所以我的建议是:先试Harmony,如果整合效果不理想,再换Seurat的整合流程

5. 一些实战中的经验与避坑建议

5.1 关于基因过滤的时机

我见过不少人在SCTransform之后才做基因过滤,比如去掉线粒体基因再跑PCA。这个顺序其实不太对。SCTransform的模型拟合会受到所有输入基因的影响,如果线粒体基因占比很高,模型会把大量方差分配给这些基因,导致高变基因的选择偏向技术噪音。

正确的顺序是:先过滤基因,再跑SCTransform。过滤的标准我一般用:表达细胞数少于3的基因去掉,线粒体基因去掉,核糖体基因去掉。但有一个例外:如果你研究的生物学问题本身就和线粒体功能或核糖体功能相关,那就不能去掉这些基因。比如研究线粒体疾病或者核糖体病,这些基因反而是核心。

5.2 关于vars.to.regress的选择

SCTransformvars.to.regressScaleDatavars.to.regress在行为上有本质区别。前者是把变量作为协变量加入GLM模型,后者是对标准化后的数据做线性回归。这意味着SCTransform的回归更温和,不会把数据“洗”得太干净。

我一般只回归percent.mt,偶尔会回归nCount_RNA(如果测序深度差异特别大)。但不建议回归太多变量,因为每回归一个变量,模型就多一个参数,过拟合的风险就增加一分。而且回归太多变量可能会把生物学信号也回归掉,比如细胞周期基因和细胞类型是相关的,回归细胞周期可能会削弱细胞类型的差异。

5.3 关于return.only.var.genes的取舍

默认return.only.var.genes = TRUE,只返回3000个高变基因的标准化数据。这对于大多数下游分析(PCA、聚类、UMAP)已经足够。但如果你需要做基因集打分(比如用AddModuleScore)或者通路分析,就需要全部基因的标准化数据。

我的做法是:先用默认参数跑一遍SCTransform做降维聚类,如果需要全部基因的标准化数据,再用GetResidual按需计算。这样既节省了内存,又保证了灵活性。GetResidual的计算速度很快,对几百个基因的基因集,几秒钟就能算完。

5.4 关于随机种子的设置

SCTransformncells抽样是随机的,所以每次跑的结果会有微小差异。在需要严格复现的项目里,一定要固定seed.use。我一般用默认值1448145,这个数字没什么特殊含义,就是Seurat作者随手设的。但如果你在跑多个样本,建议每个样本用不同的种子,避免抽样偏差。

还有一个细节:RunPCARunUMAP也有随机性,虽然影响比SCTransform小,但在严格复现的场景下也需要固定种子。我一般会在流程开始的时候设一个全局种子:

set.seed(12345)

然后在每个关键步骤(SCTransformRunPCARunUMAPFindClusters)都显式指定seed.use参数。

5.5 关于内存管理

SCTransform是内存消耗大户。一个5万细胞、3000高变基因的Seurat对象,scale.data矩阵大约是50000 × 3000 × 8字节 ≈ 1.2GB。如果return.only.var.genes = FALSE,基因数变成20000,内存占用就变成8GB。再加上原始counts矩阵、data矩阵、PCA结果、UMAP结果,总内存占用可能超过20GB。

我的建议是:如果内存有限,一定要设return.only.var.genes = TRUE。另外,在跑完SCTransform之后,可以用DietSeurat函数精简Seurat对象,去掉不必要的assay和slot:

# 精简Seurat对象 seurat_obj <- DietSeurat( seurat_obj, assays = "SCT", dimreducs = c("pca", "umap"), graphs = "SCT_nn", misc = FALSE )

DietSeurat会保留指定的assay和降维结果,去掉其他内容。我一般在SCTransformRunPCA之后跑一次DietSeurat,可以节省30%-50%的内存。

5.6 关于和BCR单细胞分析的衔接

如果你做的是BCR单细胞分析(B细胞受体分析),SCTransform标准化后的数据可以和BCR克隆型信息结合。具体做法是:先用SCTransform标准化表达数据,做降维聚类,然后把BCR克隆型信息映射到UMAP上。

# 假设bcr_data是一个包含克隆型信息的数据框 # 包含cell_id和clone_id两列 # 把克隆型信息加到Seurat对象的meta.data里 seurat_obj$clone_id <- bcr_data$clone_id[match(colnames(seurat_obj), bcr_data$cell_id)] # 在UMAP上可视化克隆型 DimPlot(seurat_obj, reduction = "umap", group.by = "clone_id")

这里有一个关键点:BCR分析通常只关注B细胞,所以你可能需要先从所有细胞中提取B细胞,再单独做SCTransform。因为B细胞的总RNA含量和其他细胞类型差异较大,如果混在一起标准化,B细胞的表达特征可能会被其他细胞“稀释”。

我一般会先用所有细胞做一轮SCTransform和聚类,鉴定出B细胞cluster,然后提取B细胞子集,再单独跑一轮SCTransform。这样B细胞内部的异质性才能被充分捕捉。

5.7 关于SICAR标准化程序的参考

SICAR标准化程序是一个在单细胞分析社区里被讨论较多的标准化流程参考。它的核心思路和SCTransform有相似之处,都强调对技术噪音的建模,但在具体实现上有差异。SICAR更侧重于跨平台数据的标准化,比如把Smart-seq2和10x Genomics的数据放在一起分析。

如果你需要做跨平台整合,可以参考SICAR的思路:先对每个平台的数据单独做标准化,然后用锚点方法整合。SCTransform本身不区分平台,所以跨平台整合时,建议在每个平台的数据上分别跑SCTransform,再用FindIntegrationAnchors整合。

我在实际项目里试过把10x数据和Smart-seq2数据放在一起,用SCTransform分别标准化后再整合,效果比直接合并后标准化要好。因为Smart-seq2的测序深度远高于10x,直接合并会让模型偏向Smart-seq2的数据分布。

6. 完整流程回顾与参数速查

把整个流程串起来,从原始数据到标准化完成,再到下游分析,核心步骤和参数如下:

# 1. 基础过滤 seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5) # 2. 基因过滤 mt_genes <- grep("^MT-", rownames(seurat_obj), value = TRUE) rb_genes <- grep("^RP[SL]", rownames(seurat_obj), value = TRUE) gene_counts <- rowSums(seurat_obj[["RNA"]]@counts > 0) low_expr <- names(gene_counts[gene_counts < 3]) keep <- setdiff(rownames(seurat_obj), c(mt_genes, rb_genes, low_expr)) seurat_obj <- subset(seurat_obj, features = keep) # 3. SCTransform标准化 seurat_obj <- SCTransform( seurat_obj, assay = "RNA", new.assay.name = "SCT", ncells = 5000, variable.features.n = 3000, vars.to.regress = "percent.mt", do.scale = FALSE, do.center = TRUE, clip.range = sqrt(5000), return.only.var.genes = TRUE, seed.use = 1448145, verbose = TRUE ) # 4. PCA降维 seurat_obj <- RunPCA( seurat_obj, assay = "SCT", npcs = 50, features = VariableFeatures(seurat_obj, assay = "SCT"), seed.use = 42 ) # 5. 聚类 seurat_obj <- FindNeighbors(seurat_obj, reduction = "pca", dims = 1:30) seurat_obj <- FindClusters(seurat_obj, resolution = 0.8, random.seed = 42) # 6. UMAP可视化 seurat_obj <- RunUMAP(seurat_obj, reduction = "pca", dims = 1:30, seed.use = 42) # 7. 差异分析(用data slot) markers <- FindMarkers( seurat_obj, ident.1 = "0", ident.2 = "1", assay = "SCT", slot = "data", test.use = "wilcox" )

参数速查表:

参数推荐值说明
ncells5000-10000细胞数多时调大
variable.features.n3000一般不超过5000
vars.to.regresspercent.mt不建议回归太多
do.scaleFALSEPearson残差已方差稳定
clip.rangesqrt(ncells)极端值多时调小
return.only.var.genesTRUE内存有限时必选
npcs30-50看ElbowPlot拐点
resolution0.4-1.2多试几个值

这套流程我在多个项目里跑过,包括PBMC、脑组织、肿瘤样本,整体稳定性不错。但每个数据集都有其特殊性,参数需要根据实际情况微调。比如脑组织的细胞类型异质性高,variable.features.n可以适当调大到4000;肿瘤样本的线粒体基因占比高,percent.mt的过滤阈值可以放宽到10%。

最后分享一个我踩过的坑:有一次跑一个10万细胞的数据集,SCTransform跑了整整一个下午没结束,后来发现是return.only.var.genes设成了FALSE,导致模型要计算所有20000多个基因的残差。改成TRUE之后,时间缩短到40分钟。所以在跑大样本之前,一定要检查这个参数

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询