☰
单细胞Seurat降维与细胞标记:PCA、UMAP参数及marker判读
2026/10/27 3:42:15 网站建设 项目流程

单细胞测序的流程走到降维这一步,往往是最容易出现“数据看着没问题、结论却站不住”的阶段。前面质控、标准化、高变基因筛选、聚类该做的都做了,跑出来一个 Seurat 对象,聚类结果也给了十几个 cluster,但打开 UMAP 图一看,问题才真正开始:这些 cluster 到底对应什么细胞?PC 数取多少才算合理?UMAP 上挨着的两群能不能说它们关系近?FindAllMarkers 跑出来几千行结果,阈值怎么设才不至于把真信号过滤掉、把噪音留下来?这篇就围绕单细胞 Seurat 降维与细胞标记这条主线,把 PCA 主成分数量的确定、UMAP 与 tSNE 的参数取舍、差异基因筛选的阈值逻辑、marker 基因的人工判读这几个环节拆开讲清楚。适合已经能把 Seurat 对象跑到聚类、但对“cluster 0 到底是什么细胞”还没底的人;也适合那些结果跑出来了、却在写结果描述时不知道该怎么下笔的人。阅读前建议手边有一个跑到FindNeighbors/FindClusters之后的 Seurat 对象,边看边跑,效果会更好。

1. 降维与注释的定位:为什么这一步最容易积累隐性问题

1.1 从表达矩阵到“细胞类型”之间隔着两层抽象

很多新手会把“降维”和“聚类”混在一起理解,实际上这是两个方向完全不同的操作。聚类是在高维空间(通常是 PCA 空间)里按距离把细胞分组,输出的是一个离散的标签;降维是把高维数据压到二维或三维,输出的是坐标。前者用于“分组”,后者用于“看图”。真正决定 cluster 归属的是聚类那一步,UMAP 只是把聚类结果画出来给人看。这个认知如果不建立起来,后面所有的参数调整都会失去方向:你会去调 UMAP 的min.dist试图把两群分开,但真正应该调整的是聚类分辨率或者 PCA 维数。

我见过不少分析卡在这一步,本质原因是从业者把可视化工具当成了分析工具。UMAP 有一个很要命的特性:它保留局部邻域结构的能力很强,但对全局距离的还原非常不可靠。两个在 UMAP 上隔得很远的 cluster,在真实表达空间里可能只差几十个基因;反过来,两群在 UMAP 上贴在一起,也可能只是共享了一部分高变基因带来的投影假象。所以看图的时候心里要有一条线:UMAP 回答的是“谁和谁像”,不回答“谁和谁差多少”。

这也是为什么这一篇要先把定位讲清楚再谈参数。降维与细胞标记这两件事,本质上是把数值结果翻译成生物学语言的过程。翻译得准不准,取决于你对每一层抽象的理解深度,而不是取决于用了哪个包、哪个函数。Seurat 的参数只是工具,判断依据还得自己建立。

1.2 线性降维与非线性降维的分工不能混

PCA 是线性降维,UMAP 和 tSNE 是非线性降维,这个分类不是学术上的文字游戏,而是直接决定它们各自该在流程的哪个位置出现。PCA 的作用是把上万个基因压缩成几十个正交主成分,同时尽可能保留原始方差。它的核心价值在于“去冗余”:单细胞数据里几千个高变基因之间存在大量共线性,直接拿去算距离会被重复信息拖垮,所以先用 PCA 把维度降到 10 到 50 之间,再在 PCA 空间里做聚类和邻居图构建,这是标准做法。

UMAP 和 tSNE 的位置在 PCA 之后,用途是可视化。它们接收的是 PCA 的坐标矩阵,而不是原始表达矩阵。这个细节很多人会忽略:RunUMAP默认用的就是dims = 1:30这样的主成分范围,而不是原始基因。所以如果你的 PCA 维数选得离谱,UMAP 的图一定跟着失真,这两者是串联关系,不是并联关系。

打个比方,PCA 像是把一本厚书按主题压缩成几十页摘要,UMAP 是把这几十页摘要画成一张思维导图。摘要压缩得不好,导图画得再漂亮也没用;但导图画得好不好看,也不改变摘要本身的内容准确性。理解了这个层级关系,后面调参时就不会乱:想改变分群结果,去动 PCA 维数和聚类分辨率;只想让图更好看,才去动 UMAP 的邻域参数。

1.3 一个被反复踩的坑:把 UMAP 当成聚类结果的验证

有个很典型的做法:跑完聚类得到 cluster 编号,画 UMAP,发现 cluster 3 和 cluster 7 在图上混在一起,于是判断“聚类有问题”,回头调分辨率重跑。这个反应本身不算错,但判断依据是错的。UMAP 上混在一起,可能只是因为这批细胞的差异集中在少数几个主成分上,而 UMAP 的邻域参数把这点差异抹平了;也可能是因为聚类分辨率偏高,本来就该拆成两群的细胞被强行分开,然后又因为投影关系在二维图上重叠。

正确的做法是回到 PCA 空间去看。可以用DimPlot(reduction = "pca", dims = 1:2)和dims = 3:4分别看一下,也可以直接算 cluster 之间的平均距离矩阵,或者用FindAllMarkers的结果判断两群之间有没有稳定的差异基因。如果差异基因明确且成规模,那分群是成立的,UMAP 上混着只是显示问题;如果差异基因只有零星一两个、且 log2FC 很低,那才是聚类过分割,需要降低分辨率或者减少 PCA 维数。

注意:UMAP 的坐标数值没有任何生物学意义,横纵轴不标单位是对的,任何基于“两个 cluster 在 UMAP 上的欧氏距离”做定量结论的说法都需要打问号。

2. PCA 主成分数量的确定:三种方法和一个经验区间

2.1 为什么主成分数是个必须认真对待的参数

RunPCA跑完之后,dims这个参数会一路传到FindNeighbors、FindClusters、RunUMAP三个函数里。它决定了后续所有步骤在多少维的空间里计算距离。取得太少,会把真实的生物学差异当成噪音丢掉,表现是聚类结果过于粗糙、多个明显不同的细胞群体被合并成一群;取得太多,会把技术噪音、批次效应、测序深度差异这些不想要的信息也纳入计算,表现是聚类过度分割、出现一堆只有几十个细胞的小 cluster,且这些 cluster 的 marker 基因往往是线粒体基因或核糖体基因。

这个参数的敏感度比很多人想象的高。同样一份数据,dims = 1:10和dims = 1:40跑出来的 cluster 数量可能相差一倍以上。所以“随便取 30”这种做法在小数据量上还行,数据量一大就必须看图表来定。下面三种方法是我实际用下来最稳的组合。

2.2 ElbowPlot:最快的第一轮筛选

ElbowPlot(obj, ndims = 50)画的是每个主成分的标准差,也就是这个主成分承载了多少方差。曲线通常会在某个位置明显变平,那个拐点就是建议的主成分数。这个方法的好处是零成本,RunPCA跑完后直接就能画,不需要额外计算。

实操上有几个细节。第一,不要只看图的形状就下结论,把ElbowPlot的返回值(其实是个 ggplot 对象)对应的obj[["pca"]]@stdev拿出来看数值会更准。比如前 10 个主成分的标准差从 12 掉到 8,10 到 20 之间从 8 缓降到 6,20 之后基本在 5.5 附近持平,那拐点就在 18 到 20 之间,取dims = 1:20是合理区间。第二,标准差下降的绝对幅度比形状更重要,有的数据曲线很平缓,看不出明显拐点,这时候就要靠累计方差贡献来兜底。

我一般会把ElbowPlot当成第一道粗筛,它给出的区间上下浮动 5 都是可接受的。真正定下来还得看后面的稳定性测试。这里有个小技巧:把ElbowPlot和VizDimLoadings配合用,看前几个主成分在哪些基因上载荷高。如果 PC1 和 PC2 的 top 基因是明显的细胞类型 marker(比如免疫球蛋白基因、血红蛋白基因),说明主成分抓到了生物学信号;如果 PC1 的 top 基因全是MALAT1、RPL系列、MT-开头的线粒体基因,说明这个主成分主要反映技术差异,这时候即便方差很高,也不该把它当成主要的判别维度。

2.3 累计方差贡献与稳定性测试

累计方差贡献的做法是把每个主成分的方差占比累加,看累计到多少时接近 80% 到 90%。Seurat 里没有直接给这个数的函数,需要自己算:

pca_stdev <- obj[["pca"]]@stdev var_ratio <- pca_stdev^2 / sum(pca_stdev^2) cum_ratio <- cumsum(var_ratio) plot(cum_ratio[1:50], type = "b", xlab = "PC", ylab = "Cumulative variance ratio") which(cum_ratio >= 0.9)[1]

这里有个容易被忽略的点:单细胞数据的累计方差曲线通常爬升得非常慢,因为基因数动辄两万,前 50 个主成分的累计贡献可能只有 30% 到 40%。所以用“累计达到 90%”这个标准在很多数据集上会得出一个荒谬的大数字,比如 200 个主成分,这显然不能用。我的经验是把标准放宽:只要累计贡献曲线的斜率明显变缓(二阶差分接近零),就可以停。对多数组织样本来说,这个位置落在 15 到 30 之间;数据异质性特别强的样本(比如包含大量不同分化阶段的细胞)可能要到 30 到 40。

稳定性测试是我最推荐的收尾方法,虽然要多跑几次。做法是取三个候选值,比如dims = 1:20、1:30、1:40,分别跑FindClusters,比较 cluster 数量和 marker 基因的一致性。如果 20 和 30 跑出来的 cluster 数量只差一两个、且多出来的那个 cluster 有明确的 marker,那取 30;如果 30 和 40 跑出来的结果里多出一堆只有几十个细胞、marker 是线粒体基因的碎 cluster,说明 40 已经过拟合,退到 30。这套流程跑下来大概多花十几分钟,但能省掉后面反复返工的时间。

2.4 PCA 空间里的几个必查项

定完维数之后,别急着往下走,先在 PCA 空间里做几项体检。第一项是看DimPlot(obj, reduction = "pca"),如果细胞在 PC1-PC2 平面上呈现出与某个技术变量强烈相关的分层(比如按样本来源、按测序批次分开),说明存在未校正的批次效应,这时候应该回到前面做整合,而不是硬着头皮往下跑。

第二项是看 PCA 的前几个主成分与已知协变量的相关性。比如把每个细胞的nCount_RNA、nFeature_RNA、percent.mt拿出来,算它们与 PC1 到 PC5 的相关系数。如果 PC1 与percent.mt的相关系数超过 0.5,说明第一大主成分主要反映线粒体含量差异,这通常是细胞状态或质量差异,不是细胞类型差异。这种情况下要么回头加强质控,要么在ScaleData时把这些协变量回归掉(vars.to.regress),但要谨慎,回归过猛会把真实的生物学信号一起抹掉。

第三项是检查RunPCA时有没有设features参数。默认用的是前面FindVariableFeatures选出的高变基因,这个默认值一般没问题。但如果你在找高变基因时没有正确设置nfeatures,比如取了 2000 个高变基因但其中一半是核糖体基因,那 PCA 的质量会明显下降。这种情况下的补救办法是重新找高变基因,用FindVariableFeatures(selection.method = "vst", nfeatures = 2000)之后,手动把MT-、RP[SL]、HB[AB]开头的基因从高变基因列表里剔除,再跑 PCA。

3. UMAP 与 tSNE:参数选择的实际影响

3.1 UMAP 的三个关键参数

RunUMAP里真正影响图形形态的参数主要是三个:n.neighbors、min.dist、spread。n.neighbors控制局部与全局的平衡,值越大,UMAP 越倾向于保留全局结构,图形会更“团”,各群体之间的相对位置更接近真实关系;值越小,越强调局部邻域,图形会碎成更多小岛,同一群体内部的结构细节更丰富。默认值 30 对大多数数据集是合适的,但如果细胞数超过十万,可以适当调大到 50;如果细胞数只有几千,调到 15 到 20 会让细节更清楚。

min.dist控制点与点之间的最小距离,默认 0.3。这个参数直接决定图的“紧凑程度”:调到 0.1 以下,同一 cluster 内的点会挤成一团,看起来更干净,但会掩盖亚群结构;调到 0.5 以上,点会散开,适合观察过渡态细胞,但整体会显得杂乱。spread与min.dist联动,默认 1,一般不需要动,除非你要调整整个图的视觉密度。

实际操作里我建议这样调:先跑默认参数出一版图作为基准,如果发现某个你确信存在的亚群在图上被压扁了、看不出来,就把min.dist降到 0.1 到 0.2 再跑一次;如果发现整体图太散、群体边界模糊,就把n.neighbors调大。每次只动一个参数,跑完对比,这样能建立起对参数敏感度的直观感受。切忌一次改三个参数,改完之后根本不知道是哪个起的作用。

3.2 tSNE 在什么场景下仍然值得用

现在多数人默认用 UMAP,tSNE 用得少了,但它并没有被淘汰。tSNE 的优势在局部结构:同一细胞类型内部的亚结构、连续分化过程中的过渡态,tSNE 往往展示得比 UMAP 更清晰。它的核心参数是perplexity,默认 30,含义是每个点考虑的有效邻居数量。这个值对图形影响极大,经验范围是 5 到 50,一般认为取细胞数的平方根量级比较合适,比如 5000 个细胞取 30 左右,20000 个细胞可以取 50。取太小(比如 5),图会碎成无数小簇,每个簇只有几个点;取太大(比如 100 以上),图会糊成一团,所有细胞挤在一起看不出结构。

tSNE 的代价是慢和不可重现的全局结构。RunTSNE在大数据集上跑十几分钟很正常,而且每次跑出来的全局布局都不一样,只有局部邻域关系是稳定的。所以如果你的目的是给论文配一张主图,UMAP 更合适;如果是要论证某个亚群确实独立存在于主群之外,tSNE 在较高 perplexity 下给出的局部紧凑性反而更有说服力。

我个人的习惯是两个都跑,写在方法和结果里时以 UMAP 为主图,tSNE 作为补充放在附件。这样做的好处是审稿人或同行提出“你这个分群是不是 UMAP 参数的产物”时,可以直接拿出第二种降维方法的一致性证据。

3.3 降维图上的距离该怎么解读

前面提过 UMAP 的全局距离不可靠,这里展开说一下具体到什么程度。有研究做过系统性比较,发现在 UMAP 上两个群体的中心距离与它们在真实高维空间中的距离之间的相关性只有中等水平,相关系数在 0.5 到 0.7 之间浮动,具体取决于数据本身的结构。也就是说,你不能看着图说“A 群离 B 群比离 C 群近,所以 A 和 B 更相似”,这个推断有一半概率是错的。

能做的解读是:同一个 cluster 内部的点聚集程度可以粗略反映该群体的异质性,聚得越散说明内部差异越大;两个 cluster 之间存在大量中间过渡点时,提示可能存在分化连续体或者双细胞。不能做的解读是:cluster 之间的相对位置、图的整体形状、轴的走向。这些都属于“看起来有信息但其实没有”的范畴。

有个具体的判断技巧:把DimPlot的group.by换成某个特定 marker 基因的表达量(用FeaturePlot),看这个基因的表达是不是沿着 UMAP 上的一条连续带分布。如果是,说明存在分化轨迹,这时候可以考虑做拟时序分析;如果表达是离散地集中在几个孤立的 cluster 里,说明是明确的细胞类型差异。这个判断比看 cluster 位置靠谱得多。

4. 差异基因筛选与 marker 基因的判定

4.1 FindAllMarkers 的参数逻辑

FindAllMarkers是这套流程里计算量最大的函数,参数设不好要么跑得慢,要么结果没法用。核心参数有四个:test.use、logfc.threshold、min.pct、only.pos。

test.use默认是 Wilcoxon 秩和检验,这个是单细胞领域的常规选择,因为它不假设数据服从正态分布,而单细胞表达数据是典型的零膨胀分布,大量基因在多数细胞里表达量为零。其他选项里,"bimod"是似然比检验,对低表达基因更敏感但假阳性也更多;"t"是 t 检验,速度快但对分布假设敏感;"negbinom"是负二项模型,理论上更贴合测序数据的特性,但计算量大很多,几万个细胞跑下来可能要几个小时。我的建议是默认用 Wilcoxon,只有在明确需要和已发表结果对齐时才换。

logfc.threshold默认 0.25,含义是要求基因在目标 cluster 与其他 cluster 之间的平均 log2 倍数变化超过这个值才进入检验。这个默认值偏宽松,跑出来的结果动辄几千行。实际使用中我一般设到 0.5 甚至 1.0。为什么?因为 log2FC = 0.5 意味着表达量差 1.41 倍,这个差异在单细胞数据里往往落在技术噪音范围内;设到 1.0 意味着差 2 倍,才算是比较明确的差异。当然,设太高会漏掉一些真实的低幅度差异基因,所以这一步要结合后面的结果过滤一起考虑,不要指望一个参数解决所有问题。

min.pct默认 0.1,含义是基因至少在 10% 的细胞里有检出才参与检验。这个参数的作用是过滤掉那些只在极少数细胞里表达的基因,避免它们因为偶然的高表达而产生虚假的倍数变化。对于稀有细胞类型,这个值可能要降到 0.05,否则稀有群体的 marker 会被全部滤掉。反过来,如果你的数据里混了大量低质量细胞,把min.pct提到 0.25 能有效减少噪音。

only.pos = TRUE是我强烈建议打开的。默认是 FALSE,会同时返回上调基因和下调基因,导致结果表里同一个基因可能在 cluster A 里是上调、在 cluster B 里是下调,读起来非常混乱。打开之后只返回上调基因,每个 cluster 的 marker 列表干净很多。

markers <- FindAllMarkers( obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.5, test.use = "wilcox" )

4.2 结果表格怎么过滤才算合理

跑出来的markers是个数据框,列包括p_val、avg_log2FC、pct.1、pct.2、p_val_adj、cluster、gene。这里最容易犯的错是只看p_val不看p_val_adj。单细胞数据里每个 cluster 要检验上万个基因,多重检验校正之后p_val_adj会比p_val大好几个数量级。用p_val < 0.05作为过滤标准,几乎等于不过滤。

合理的过滤组合是这样的:p_val_adj < 0.05保证统计显著性;avg_log2FC > 0.5保证效应量;pct.1 - pct.2 > 0.2保证这个基因在目标 cluster 里的检出率明显高于其他 cluster。第三条经常被忽略,但很关键。有些基因的pct.1是 0.3、pct.2是 0.28,虽然统计上显著(因为细胞数多),但实际检出率只差 2 个百分点,这种基因作为 marker 是没有说服力的。

提示:过滤阈值不要一刀切。对于细胞数少的 cluster(比如只有 50 个细胞),p_val_adj很难达到显著,这时候可以放宽到 0.1,同时把效应量阈值提高。对于细胞数上千的大 cluster,标准可以严格一些。

另外一个实用技巧是按 cluster 分组取 top N。用markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC)取每个 cluster 的前 20 个基因,做成热图或者 dotplot,这是判断 cluster 身份最直观的方式。取 20 这个数不是随便定的:太少可能漏掉关键 marker,太多会混入一堆功能不明确的基因,反而干扰判断。

4.3 常见细胞类型的 marker 判读思路

判读 marker 这件事,靠的不是死记硬背基因列表,而是理解“这个基因在什么细胞里该高表达”。下面这张表是我常用的外周血和常见组织样本的参考,实际项目里还要结合组织来源调整。

细胞类型常用 marker判读要点
T 细胞CD3D, CD3E, CD3G, TRAC, IL7RCD3 系列是最稳的,但有 dropout,建议三个一起看
CD4 TIL7R, CD4, CCR7CD4 本身表达量低,单看容易误判
CD8 TCD8A, CD8B, GZMK效应态会高表达颗粒酶基因
NK 细胞NKG7, GNLY, KLRD1, PRF1与 CD8 T 有重叠,关键区别是 CD3 阴性
B 细胞MS4A1, CD79A, CD79B, CD19MS4A1(CD20)最特异
浆细胞MZB1, JCHAIN, DERL3同时高表达免疫球蛋白基因
单核细胞CD14, LYZ, VCAN, S100A8, S100A9S100A8/A9 在中性粒细胞里也高,需结合 CD14
巨噬细胞CD68, CD163, C1QA, C1QB组织样本里常见,与单核细胞靠 C1Q 系列区分
树突状细胞FCER1A, CST3, CD1C细胞数通常很少,容易在聚类时被并入单核
血小板PPBP, PF4, ITGA2B常作为污染指标出现
红细胞HBB, HBA1, HBA2高表达说明红细胞污染未清除干净
增殖细胞MKI67, TOP2A, STMN1可以跨细胞类型出现,不是独立细胞类型

判读时的核心逻辑是“不看单个基因,看基因组合”。比如一个 cluster 同时高表达CD3D、CD3E、IL7R但CD8A阴性,那基本可以定为 CD4 T 细胞;如果CD3D阴性而NKG7、GNLY高表达,那就是 NK。单个基因的表达受 dropout 影响很大,尤其在小 cluster 里,某个 marker 在FeaturePlot上看起来一片灰,不代表这个群体不表达它,可能只是这个基因在这批细胞里检出率低。

还有一个经常被忽略的点:阴性 marker 和阳性 marker 同样重要。判断一个 cluster 是不是单核细胞,除了看CD14、LYZ高不高,还要确认CD3D、MS4A1是阴性的。如果两个都阳性,那要么是双细胞,要么是这个 cluster 混了两群细胞,得回去重新审视聚类分辨率。

4.4 可视化组合:FeaturePlot、VlnPlot、DotPlot 各管什么

FeaturePlot看空间分布,适合快速判断一个基因是不是集中在某个 cluster 上。它的坑在于颜色标尺的默认设置会把低表达区域压成灰色,导致看起来像“完全不表达”。调min.cutoff和max.cutoff能改善,比如设成min.cutoff = "q10"、max.cutoff = "q90",用分位数代替绝对值,视觉效果会好很多。

VlnPlot看表达量分布,适合确认某个基因在目标 cluster 里是不是整体高表达、还是只有少数细胞高表达。后一种情况往往意味着这个基因不是这个群体的特征,而是某个亚群或者细胞状态的特征。看VlnPlot时要注意纵轴,不同基因的表达量范围差别很大,放在一起对比时最好加pt.size = 0避免点把图形糊住。

DotPlot是三者里信息密度最高的,一个点同时编码了两个信息:点的大小是检出率(pct.exp),点的颜色是平均表达量(avg.exp.scaled)。判读的时候先看点的大小,检出率低于 20% 的基本可以不看;再看颜色深浅,最后横向对比不同 cluster。DotPlot还有个实用功能是可以分组显示,用group.by把 cluster 按细胞类型归并之后再看,能更清楚地看出细胞类型层面的 marker 格局。

DotPlot(obj, features = c("CD3D","CD3E","MS4A1","CD79A","CD14","LYZ", "NKG7","GNLY","FCER1A","PPBP","HBB"), group.by = "seurat_clusters") + RotatedAxis()

5. 细胞类型注释的实施路径

5.1 人工注释的标准流程

注释这一步没有标准答案,但有标准流程。我的做法分四步走。第一步,把FindAllMarkers的结果按 cluster 分组,每个 cluster 取 top 20 基因,逐个 cluster 过一遍,对每个 cluster 给出一个初步猜测。这一步不需要精确,只需要判断“大概是淋巴系还是髓系”。第二步,挑出每个猜测对应的经典 marker 基因,用DotPlot一次性画出所有 cluster 在所有候选 marker 上的表达格局。这一步的目的是验证第一步的猜测,同时发现猜错的地方。第三步,对存疑的 cluster 单独深挖,用FeaturePlot看几个关键基因的空间分布,判断是一个群体还是两个混在一起。第四步,把所有 cluster 的注释结果写成一个命名向量,用RenameIdents批量改名,同时保留原始的 cluster 编号作为一列 metadata。

这个流程里最关键的是第二步。DotPlot一次能看十几个 marker 在十几个 cluster 上的格局,比一个个画FeaturePlot高效得多,而且能直接看出“这个 marker 在两个 cluster 上都高”这种交叉情况。我习惯把候选 marker 按细胞类型分成几组,画三到四张DotPlot,每张覆盖一组相关的细胞类型,这样判读的时候思路更清晰。

第四步的命名向量写法有个小技巧:不要用中文名直接做 ident,用英文或缩写,中文名放在 metadata 的另一列里。原因是后续很多函数对 ident 的字符处理会有编码问题,用英文能避免莫名其妙的报错。

new_cluster_ids <- c( "0" = "CD4 T", "1" = "CD14 Mono", "2" = "CD8 T", "3" = "B", "4" = "NK", "5" = "FCGR3A Mono", "6" = "DC", "7" = "Platelet" ) obj <- RenameIdents(obj, new_cluster_ids) obj$cell_type <- Idents(obj)

5.2 打分法:AddModuleScore 的正确用法

AddModuleScore的思路是给每个细胞算一个“模块得分”:先算这组基因的平均表达量,再减去随机背景基因集的平均表达量,得到相对富集程度。这个函数的价值在于处理那些 marker 表达不集中、但整体基因集有协同变化的情况,比如细胞周期、代谢状态、通路活性。

用法上有个容易踩的坑:AddModuleScore的ctrl参数默认是 100,含义是随机抽 100 组背景基因,这个值对结果影响不大但会显著影响运行时间。数据量大时可以把ctrl降到 20 到 50,结果差别很小。另一个坑是基因列表里如果有重复基因,函数会静默去重,但去重后的基因数变了,得分就不可比了,所以用之前先unique()一下。

AddModuleScore输出的列名是Cluster1、Cluster2这种自动编号,容易和 cluster 编号混淆。建议加完分之后立刻改个有意义的名字,或者直接在features参数里给列表命名,Seurat 会把这个名字带到输出列名里(需要较新版本)。

需要说明的是,打分法适合做辅助判断,不适合单独作为注释依据。因为打分结果受基因集选择、背景基因选择、数据标准化方式的多重影响,同一个细胞类型用不同基因集打分,排名可能不一致。我的用法是把它作为DotPlot的补充:当DotPlot上某个 cluster 的 marker 表现模棱两可时,用打分法看整体趋势,两者一致才下结论。

5.3 三类疑难 cluster 的处理

第一类是双细胞 cluster。特征是细胞数少、marker 基因同时涵盖两种不相关的细胞类型、nFeature_RNA和nCount_RNA明显偏高。比如一个 cluster 同时高表达CD3D和CD14,那基本可以判定是 T 细胞和单核细胞的双细胞。处理方式是标记出来,在后续分析里剔除,或者回到前面用专门的双细胞检测工具重新过滤。不要试图把它注释成某种“新型细胞类型”,那是自欺欺人。

第二类是低质量 cluster。特征是percent.mt高、nFeature_RNA低、marker 基因散乱无明显富集。这种情况通常是前面质控标准偏松,或者样本本身存在大量濒死细胞。处理方式是在 metadata 里标记,后续分析排除。如果这类细胞占比超过 10%,建议回到质控那一步重新设定阈值重跑,而不是在这里打补丁。

第三类是状态相关 cluster。特征是 marker 基因不属于任何特定细胞类型,而是细胞周期基因(MKI67、TOP2A)、应激基因(HSPA1A、HSPB1)、干扰素响应基因(ISG15、IFI6)。这类 cluster 反映的是细胞状态而非细胞类型,注释时应该标注成“CD4 T(增殖态)”这种复合形式,而不是单独列成一种细胞类型。这一点在写结果时特别重要,因为它会影响后续细胞类型比例统计的解读。

6. 常见问题与排查速查表

6.1 报错与异常现象对照

现象可能原因排查与处理
RunUMAP报错找不到 reduction没跑 PCA 或 reduction 名字写错检查obj@reductions里有没有pca,reduction = "pca"大小写要一致
FindAllMarkers跑几小时不动细胞数多、min.pct太低、没开并行提高min.pct到 0.25,用future并行,或改用加速方案
所有 cluster 的 top marker 都是线粒体基因质控过松或 PCA 维数过高回到质控检查percent.mt阈值,降低dims
某个 cluster 没有任何显著 markercluster 太小或与邻近 cluster 差异微弱降低logfc.threshold到 0.25,提高聚类分辨率重新分群
UMAP 图所有细胞糊成一团n.neighbors过大或 PCA 维数过少降到 15 到 20,或增加dims范围
同一细胞类型散落在多个 cluster聚类分辨率过高或存在批次效应降resolution到 0.3 到 0.5,检查是否需整合
DotPlot报基因不存在基因名大小写或物种不匹配用grep在rownames(obj)里模糊搜索确认
注释后细胞比例明显异常双细胞或低质量细胞未剔除检查nFeature_RNA分布和双细胞比例

6.2 参数速查与推荐区间

参数函数默认值推荐区间说明
dimsRunUMAP/FindNeighbors1:30(需指定)15–40结合 ElbowPlot 和稳定性测试确定
resolutionFindClusters0.80.3–1.2细胞数越多可适当提高
n.neighborsRunUMAP3015–50小数据偏小,大数据偏大
min.distRunUMAP0.30.1–0.5想突出亚群就调小
perplexityRunTSNE305–50约为细胞数的平方根量级
logfc.thresholdFindAllMarkers0.250.5–1.0配合效应量过滤
min.pctFindAllMarkers0.10.1–0.25稀有细胞群降到 0.05
nfeaturesFindVariableFeatures20002000–3000超过 3000 容易引入噪音

6.3 我从反复返工里总结的几条经验

第一条,每次只改一个参数。这是老生常谈,但真正做到的人不多。因为降维和聚类的参数之间有交互作用,一次改两个,结果变好了你不知道是谁的功劳,结果变差了你也不知道该退哪个。养成习惯:跑一版、记录参数、画图、对比、改一个、再跑。

第二条,把每一版的图和参数都存下来。用带日期和参数组合的文件名,比如umap_dims30_res0.8_nn30.png,配合一个简单的文本记录本写清楚这版和上一版的差异。这个习惯在项目做到第三周、你已经跑了二十版图的时候会救你的命。我吃过这个亏:客户问“为什么最后选了 0.6 的分辨率”,我翻遍文件夹找不到当时的对比图,只能重跑。

第三条,注释结果要留两列。一列是seurat_clusters原始编号,一列是注释后的cell_type。永远保留原始编号,因为后续如果发现某个 cluster 注释错了,或者数据量增加需要重新聚类,原始编号是你回溯的唯一线索。只保留注释名、丢掉原始编号这种做法,会在返工时让你多花一整天。

第四条,不要迷信自动注释工具。现在有不少基于参考数据集的自动注释方案,跑起来很方便,但它们的效果高度依赖参考数据集的匹配度。跨组织、跨物种、跨疾病状态的数据,自动注释的结果经常和实际 marker 表达对不上。正确用法是把它当成一个“提出假设”的工具,用人工 marker 判读去验证。两者一致的可以直接用,不一致的以人工判读为准,并把分歧记录下来,这本身就是有价值的信息。

第五条,留意负向证据。判读细胞类型时,看到“这个 cluster 高表达LYZ”就下结论是单核细胞,这是最容易出错的地方,因为LYZ在很多髓系细胞里都高。真正让结论站得住的是“高表达LYZ且CD3D、MS4A1、NKG7都是阴性”这种组合证据。养成习惯,每确定一个阳性 marker,顺手确认两三个该细胞类型不该有的阴性 marker,误判率会明显下降。

最后分享一个判断自己这套流程是否可靠的土办法:把最终注释好的对象存一份,隔两周再打开,不看之前的注释,重新按流程走一遍,看两次结果的一致率有多少。一致率在 90% 以上,说明你的判读标准是稳定的;低于 80%,说明有些 cluster 的归属其实是在赌,需要回头把那些 cluster 的 marker 重新看一遍。这个方法听起来笨,但对建立自己的判读标准非常有效,尤其是刚开始独立做分析的时候,做上三五次,心里的那杆秤就立起来了。

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

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

立即咨询