单细胞测序分析走到降维与细胞标记这一步,往往是最考验判断力的阶段,前面QC、标准化、高变基因筛选都做完了,数据干净了,但真正的“这堆细胞到底是谁”还得靠降维和聚类来解决。Seurat这套流程走到第四步,RunPCA、RunUMAP、FindClusters、FindAllMarkers这几个函数谁都会敲,参数也就那么几个,可真正拉开分析质量差距的,恰恰是PC数量怎么定、resolution怎么调、marker阈值怎么设这些看似随手的决定。我见过太多数据本身不错、结果却因为降维参数选错而把真实亚群糊成一团的例子。这篇就把降维与细胞标记这一段的实操逻辑、参数取舍、踩过的坑一次说透,适合刚跑通Seurat基础流程、准备对自己的数据下手做生物学解释的读者,也适合跑过几次但对参数选择还心里没底的同学。
1. 降维动作开始前,Validity检查才是真正的起跑线
很多教程会把NormalizeData、FindVariableFeatures、ScaleData这三步带过去,直接讲RunPCA,但从实操经验看,降维结果好不好,七成取决于这三步做得稳不稳。尤其是ScaleData这一步,看着只是标准化,实际决定了后续PCA里哪些信号被放大、哪些被压制。
1.1 高变基因选2000还是3000,取决于组织复杂度
FindVariableFeatures默认nfeatures=2000,选vst方法。这个数在血液、PBMC这种细胞类型相对清晰的样本里够用,但如果是肿瘤微环境、发育组织、脑组织这类细胞状态跨度大的样本,2000个基因容易把稀有细胞类型的特征基因筛掉。我的习惯是:常规样本用2000,复杂组织直接拉到3000,并且把selection.method设成“vst”,不要用“mean.var.plot”,后者对低表达基因不友好。
还有一个细节容易被忽略:FindVariableFeatures返回的是一张排序表,可以用VariableFeaturePlot画出来看,再用LabelPoints把前20个高变基因标出来。如果这20个基因里出现大量线粒体基因、核糖体基因、或者某种单一细胞类型的marker霸榜,说明数据里可能存在某一类细胞占比过高,或者有残留的死亡细胞,这时候应该在QC阶段回头处理,而不是硬着头皮往下走。
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 3000) top20 <- head(VariableFeatures(sce), 20) VariableFeaturePlot(sce) %>% LabelPoints(points = top20, repel = TRUE)1.2 ScaleData回归协变量,不是越多越好
ScaleData的核心是把每个基因的表达值做z-score标准化,让不同基因在PCA里有可比性。但真正决定这一步要不要加vars.to.regress的,是数据里有没有明显的技术噪音源。常见的回归对象是percent.mt(线粒体比例)和细胞周期打分(S.Score、G2M.Score)。
这里有个非常容易踩的坑:回归线粒体比例确实能削弱死亡细胞带来的干扰,但如果你的生物学问题本身就和线粒体含量相关(比如代谢重编程、某些分化状态),回归掉它等于把真实信号一起抹了。细胞周期同理,研究增殖相关问题时不能回归。我的判断原则是:先不回归跑一遍,再回归跑一遍,对比两次UMAP的cluster结构,如果回归后cluster数量和解构逻辑变化很大,就要慎重;如果只是把几个明显由细胞周期驱动的簇合并掉,那回归是对的。
sce <- ScaleData(sce, vars.to.regress = c("percent.mt", "S.Score", "G2M.Score"))1.3 基因名格式和重复名,PCA前必须清一遍
ScaleData跑完如果报错提示“duplicate row names”或者“non-unique features”,说明基因注释表里有重复基因名(常见于不同转录本或注释版本混用)。Seurat在新版本里对重复名处理比较严格,建议在创建对象之初就用make.unique或者只保留蛋白编码基因来规避。另一个常见问题是基因名格式(ENSG还是SYMBOL),如果一行代码里混用两种格式,高变基因筛选和后续marker鉴定都会对不上号。这些看起来是琐碎细节,但出错时往往表现为“PCA跑出来一片糊”或者“marker全空”,排查起来非常浪费时间。
2. PCA的PC数量,是整个降维链条上最容易被随手定死的参数
RunPCA本身没什么可调的,真正需要动脑子的是后面RunUMAP、FindNeighbors、FindClusters里的dims参数——也就是到底取前多少个PC。这个数取少了,稀有亚群的信息被压缩掉;取多了,噪音PC混进去,聚类反而散乱。
2.1 ElbowPlot只是参考,不是判决书
几乎每篇教程都会让你画ElbowPlot(sce),看“肘部”在哪就取到哪。实际用起来你会发现,很多数据的肘部根本不明显,曲线平滑下降,你说第10个是肘部也行,说第15个也行。这时候就需要第二层判断。
我的做法是三重交叉验证:ElbowPlot看整体趋势,JackStraw虽然慢但能给出每个PC的显著性(大数据集可以抽样做),再加一个PCHeatmap看前20个PC在各个细胞上的表达模式。如果某个PC的heatmap显示出明显的批次结构(比如一半细胞全红、一半全蓝,而且和样本来源高度一致),那这个PC大概率是批次噪音,就不要再往后取了。
ElbowPlot(sce, ndims = 50) sce <- JackStraw(sce, num.replicate = 100) sce <- ScoreJackStraw(sce, dims = 1:20) JackStrawPlot(sce, dims = 1:20)2.2 dims选10还是30,用实测对比说话
网上常见的推荐是取10到20,但我在多个组织样本上实测下来,这个区间偏保守。PBMC这类细胞类型少、结构清晰的样本,dims=15到20足够;而肿瘤、脑、胚胎组织,往往要取到25到30才能把稀有亚群分开。
一个实操技巧是:先用dims=1:30跑一遍UMAP和聚类,再用dims=1:15跑一遍,把两张UMAP并排看。如果30维版本里出现了15维版本看不到的小簇,而这些小簇在marker上确实有特异性表达,说明多维是对的;如果30维版本只是把一个原本完整的簇拆成几块、marker却区分不出来,说明多的那些PC是噪音。这个对比花不了几分钟,但能省下后面返工的一整天。
2.3 PC取多了之后的典型症状
PC取过量最直接的表现是UMAP图上出现大量细碎的小岛,FindClusters给出的cluster数远超你对该组织的预期,而且FindAllMarkers跑出来每个cluster的差异基因都很少、logFC很低。这种情况下第一反应不该是“我的数据里有新细胞类型”,而应该检查dims是不是定高了。我遇到过一次,dims从30降到20,cluster数从24降到11,marker一下子清晰了,这就是典型的噪音PC污染。
3. UMAP和tSNE,选哪个以及参数背后到底在调什么
降维可视化里最常见的争论就是UMAP还是tSNE。这两个都是非线性降维,都能把高维数据压到二维,但底层机制不同,画出来的图解决的问题也不一样。
3.1 n.neighbors和min.dist,UMAP的两个手感旋钮
RunUMAP的默认参数是n.neighbors=30、min.dist=0.3、spread=1。n.neighbors控制的是每个细胞在构建邻接图时考虑多少个“邻居”,值越大,越强调全局结构,局部细节被平滑;值越小,越强调局部结构,但全局关系容易失真。min.dist控制的是低维空间里点与点之间的最小距离,值越小,同一簇内的点挤得越紧,簇与簇之间的边界越清晰。
我的调参经验是:常规分析用默认值;如果发现某个已知应该存在的稀有亚群在UMAP上被大类细胞吞掉了,把n.neighbors降到15到20试试,往往能让它冒出来;如果发现图上一堆细碎小簇、视觉上很乱,把n.neighbors升到40到50,能得到更干净的整体结构。但要注意,这些参数只影响可视化,不改变聚类结果,所以调参是为了看清结构,不要拿它当聚类依据。
sce <- RunUMAP(sce, dims = 1:20, n.neighbors = 30, min.dist = 0.3, spread = 1)3.2 tSNE的perplexity,和它那点“伪结构”问题
tSNE的perplexity大致对应“每个点期望的邻居数”,默认30。perplexity偏小,会强调局部,容易把一个连续状态拆成多个假簇;偏大,会把真实亚群揉在一起。更需要注意的是,tSNE在簇间距离上没有定量意义——两个簇在图上看得很近,不代表它们生物学上相似。我见过有人拿着tSNE图量簇间距,这就是典型的误用。
相比之下,UMAP在全局结构的保留上略好一些,簇间相对位置有一定参考价值,但也仅限参考。两者共同的铁律是:所有基于距离的定量结论,都要回到PCA空间或者表达矩阵本身去算,不能从二维图上看。
3.3 二维图不能承载的几件事
- 不能用来判断细胞数量差异:UMAP上点密度和真实细胞比例不是线性关系。
- 不能用来做组间比较:不同批次、不同参数跑出来的UMAP坐标系不可比,如果要做组间对比,应该用相同参考做整合后再降维,或者直接比较各cluster的细胞比例表。
- 不能作为聚类依据:聚类是在高维PCA空间里做的,UMAP只是画出来给你看的。
这三点听起来像废话,但真正在组会上被追问“为什么A组这个簇比B组大”时,如果拿UMAP图去回答,基本就露馅了。
4. resolution决定了你最终交出去的是几群细胞
FindClusters是这一步最关键的函数,而resolution是它唯一真正需要你深思的参数。这个参数控制的是Louvain/Leiden算法在社区发现时的粒度,值越大,切出来的cluster越多。
4.1 resolution从0.2到2.0,cluster数是怎么变的
以一个两万细胞左右的肿瘤样本为例,resolution=0.2大概出5到6个簇,0.5出8到10个,0.8出12到15个,1.2可能出20个以上,2.0能出30个以上。这不是线性关系,通常在某几个值附近会出现跳变。
实操建议是不要一上来就0.8。先跑0.5,把marker看一眼;如果发现某个大类里明显能分出亚群(比如T细胞里既有CD4又有CD8,还有调节性T的迹象),再把resolution提到0.8或1.0重跑。Seurat的FindClusters支持直接给一个序列,它会一次性生成多组结果存放在不同列里,方便对比:
sce <- FindClusters(sce, resolution = c(0.3, 0.5, 0.8, 1.0, 1.2)) # 结果存在 sce@meta.data 的 RNA_snn_res.X 列里然后用clustree画一张聚类树,能直观看出不同resolution下cluster是怎么分裂合并的,这个图对判断“该停在哪一档”特别有帮助。
4.2 Louvain还是Leiden,algorithm参数别忽略
Seurat的algorithm参数默认是1(Louvain),可选2或4(Leiden的不同实现)。两者的差别在于Leiden对Louvain的一个已知问题做了修正——Louvain可能产生内部连通性不好的“坏簇”,Leiden通过迭代 refinement 保证每个簇内部是连通的。在细胞数较多、结构复杂的数据里,我倾向于直接用Leiden(algorithm=4),跑出来的簇更稳定、可重复性更好。代价是跑得稍慢一点,但对两万细胞以内完全不是问题。
4.3 过度聚类的识别信号
过度聚类最典型的信号有三个:一是相邻cluster的top marker高度重叠,翻来覆去就是那几个基因;二是某个cluster只有几十个细胞,且marker散乱、没有明确生物学身份;三是clustree上看出某一档resolution之后簇数暴增但并没有新的marker支撑。
遇到这种情况,正确的做法是往回降resolution,而不是硬着头皮去给每个小簇编细胞类型名。我见过有人把一个20细胞的簇命名为“未知亚群X”,结果换一批数据就再也复现不出来——这种命名除了给自己添麻烦没有别的作用。
5. marker基因鉴定,FindAllMarkers的默认参数藏着不少陷阱
聚类给出编号之后,下一步就是FindAllMarkers找出每个cluster的特征基因。这个函数的默认参数对新手非常不友好,直接跑出来的结果经常是几百个基因一大片,真正有用的没几个。
5.1 默认阈值的两个问题
FindAllMarkers的默认设置是logfc.threshold=0.25、min.pct=0.1、test.use=“wilcox”。logfc阈值0.25意味着只要一个基因在目标簇里比在其他簇里高一点点就算差异基因,在细胞数多、表达噪音大的数据里,这会放进大量假阳性。min.pct=0.1又太宽松,一个基因只要在10%的细胞里表达就算数,对于那些只在少数细胞里高表达的基因会漏掉,同时又放进一堆低表达噪音。
我的常用配置是:
markers <- FindAllMarkers(sce, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.5, test.use = "wilcox")only.pos=TRUE只保留上调基因,省掉一半结果;min.pct提到0.25,logfc提到0.5,能过滤掉大部分噪音。如果数据细胞数特别多(十万以上),wilcox检验会给出极小的p值,这时候p_val_adj几乎没有区分度,更该看的是avg_log2FC和pct.1/pct.2的对比。
5.2 二次过滤,比调函数参数更重要
FindAllMarkers跑完拿到的是一张表,真正的筛选在表上做。我的过滤规则是:p_val_adj < 0.05、avg_log2FC > 0.5、pct.1 > 0.5、pct.2 < 0.3。这四个条件组合下来,每个cluster通常能留下10到30个基因,数量刚好够判断细胞类型,又不至于淹没重点。
| 过滤条件 | 作用 | 常用阈值 |
|---|---|---|
| p_val_adj | 控制多重检验假阳性 | < 0.05 |
| avg_log2FC | 保证差异幅度 | > 0.5(严格可到1.0) |
| pct.1 | 目标簇内表达比例 | > 0.5 |
| pct.2 | 其他簇表达比例 | < 0.3 |
需要提醒的是,pct.1和pct.2的阈值要结合cluster大小调整。一个只有50个细胞的小簇,pct.1很难做到0.5以上,这时候应该看pct.1 - pct.2的差值,而不是死守绝对阈值。
5.3 DotPlot和FeaturePlot必须交叉看
光看marker表格容易误判。一个基因在某簇logFC很高,可能是因为它在其他所有簇里几乎不表达(分母小导致比值大),并不代表它在该簇里表达量高。这时候要把top marker拿去画DotPlot和FeaturePlot。
DotPlot能同时看出表达比例(点大小)和平均表达量(颜色深浅),比表格直观得多。FeaturePlot则是把表达量投射到UMAP上,特别适合判断一个基因是不是只在某个簇里特异表达。我通常会对每个候选marker做这两张图,确认它在UMAP上的表达分布和它在表格里的统计结果一致,再下结论。
top_markers <- markers %>% group_by(cluster) %>% top_n(5, avg_log2FC) DotPlot(sce, features = unique(top_markers$gene)) + RotatedAxis() FeaturePlot(sce, features = c("CD3D", "CD79A", "LYZ", "NKG7"), ncol = 2)6. 从cluster编号到细胞类型,注释这一步的三种打法
cluster只是数字编号,1、2、3并不代表任何生物学意义,把它们翻译成T细胞、B细胞、单核细胞才是这一步的终点。注释方式大致分三类,各有适用场景。
6.1 手动注释,靠经典marker表打底
手动注释的核心是维护一张自己常用的marker表。不同组织、不同物种的marker不一样,但常见大类基本稳定:T细胞看CD3D/CD3E/CD8A/CD4,B细胞看CD79A/MS4A1/CD19,NK看NKG7/GNLY/KLRD1,单核和巨噬看LYZ/CD14/CD68/CD163,树突状细胞看FCER1A/CST3,成纤维看COL1A1/DCN/LUM,内皮看PECAM1/VWF,上皮看EPCAM/KRT18,浆细胞看MZB1/JCHAIN。
手动注释的流程是:先看每个cluster的top 10 marker,对照marker表初步命名;再用DotPlot把这些候选marker一次性画出来,横向对比所有cluster的表达模式;对于拿不准的cluster,回到FeaturePlot看它的分布是否符合该类细胞的预期。这种方式灵活、可控,缺点是对少见细胞类型依赖经验,而且不同人对同一簇可能给出不同命名。
6.2 自动注释工具,方便但必须人工复核
SingleR、CellTypist、scType这类工具用参考数据集做相关性打分,能给出候选细胞类型和置信度。它们在大类识别上表现不错,但有两个局限:一是参考数据集覆盖不到的细胞类型它认不出来,二是相近类型之间容易混淆(比如单核和巨噬、CD4记忆和naive)。
我的用法是把自动注释结果作为“第二意见”,和自己的手动判断对照。如果两者一致,信心就高;如果不一致,就去看那个cluster的marker到底更支持哪一方。绝不能直接拿自动注释结果贴到图里就完事,审稿人或者同行一眼就能看出那种没有经过人工核对的注释。
library(SingleR) ref <- celldex::HumanPrimaryCellAtlasData() pred <- SingleR(test = GetAssayData(sce, slot = "data"), ref = ref, labels = ref$label.main) sce$singleR_label <- pred$labels table(sce$seurat_clusters, sce$singleR_label)6.3 混合注释,效率和质量之间最务实的折中
实际项目里我基本都走混合路线:先用SingleR跑一遍拿到初步标签,再用FindAllMarkers的top基因去验证;对于自动注释置信度低的cluster,手动补充;最后用一张完整的DotPlot把所有cluster的marker模式画在一起,做整体一致性检查。
这张总DotPlot是最终交付前必须要看的图。有些cluster在单独看marker时像是某种细胞,放到全局图上一对比,发现它其实是另一类细胞的过渡态或者双细胞(doublet),这时候就需要回到上游做doublet剔除,或者把它单独标注为过渡状态。这一步不做,后面的差异分析和通路富集都可能建立在错误的细胞分组上。
7. 这一步最容易翻车的几个场景
降维和注释环节的问题,往往不是代码报错,而是结果“看起来对但其实不对”,这种问题最难自查。下面这几个是我和身边同行都遇到过的高频翻车场景。
7.1 聚类按样本来源分开,而不是按细胞类型
这是最经典的问题。UMAP上一眼看去cluster分得很清楚,但把样本来源(orig.ident)投射上去一看,每个cluster几乎对应一个样本——这说明聚类捕捉的是批次效应,不是生物学差异。多批次数据在做降维前必须做整合(比如Harmony、CCA、RPCA),整合后再跑PCA和UMAP。判断整合是否成功的标志是:UMAP上不同样本的同类细胞混在一起,而不是各占一块。
7.2 线粒体高表达细胞主导了某个cluster
如果某个cluster的top marker里出现大量MT-开头的基因,说明这个簇大概率是死亡或濒死细胞。这种情况不该在注释阶段硬给它命名,而应该回到QC阶段收紧percent.mt阈值,或者在降维前把这类细胞剔除。同理,如果一个cluster的marker全是核糖体基因(RPS/RPL系列),可能是细胞状态异常或者建库偏差,需要结合实验记录判断。
7.3 找不到marker,可能不是注释方法的问题
有些cluster跑完FindAllMarkers几乎没剩下几个显著基因,这时候不要急着换检验方法或者降阈值。常见原因有三个:一是resolution定得太高导致簇内异质性不足,任何基因都区分不开;二是这个cluster本身就是其他几个簇的混合体,内部没有统一特征;三是这个cluster的细胞数太少,统计效力不够。对应的处理分别是降resolution、检查是否该合并、以及考虑把过小的簇并入邻近簇。找到是哪一种,比盲目调参数有用得多。
7.4 命名过度自信,是注释阶段最大的隐患
最后说一个不是技术层面的坑。单细胞注释的结论天然带有不确定性,一个cluster命名为“CD8效应记忆T细胞”还是“CD8耗竭T细胞”,往往取决于几个marker的表达强弱,而这几个marker的表达又受测序深度、细胞状态、样本处理影响。我的习惯是在命名时保留一定的余地,比如标注为“T细胞(CD8为主)”,而不是一口咬定某个精细亚型;同时在文章或者报告的方法部分写清楚判断依据是哪几个marker、用的什么阈值。这样即使后续有争议,也有据可查、可复现。
我在实际项目里最深的一点体会是:降维和注释这两步,参数没有标准答案,只有“对这批数据合不合理”。同一套参数在PBMC上堪称完美,换到肿瘤样本可能就一塌糊涂。所以与其背参数,不如把每一个参数的调整都当成一次小实验——改一个值、跑一遍、对比结果、记录理由。跑的次数多了,对自己手里的数据会形成一种直觉,看一眼UMAP就知道哪里可能不对劲。这种直觉没法从教程里抄来,只能靠自己一遍遍试出来。