前几天团队里一个刚入手单细胞分析的同学跑过来问我聚类好的UMAP图上怎么把一个感兴趣的基因标出来看看它到底富集在哪些细胞类群这个问题听起来基础但实际做单细胞数据的人都知道它是整个分析流程里绕不开的一步——聚类完成后第一件事就是看已知marker基因在哪些cluster里表达或者反过来拿着一个差异筛选出来的候选基因去定位它可能的细胞身份。今天我把自己平时在UMAP图上标注基因、定位类群的那套完整做法整理出来包括底层逻辑、实际代码、配色参数和掉过的坑希望能帮到正在被这个问题卡住的人。这篇内容适合两类读者一种是刚学单细胞分析、已经把数据跑到UMAP聚类但不知道怎么解释类群的新手另一种是已经能跑流程但想系统搞清楚“为什么画出来的图总是一团蓝”“为什么基因名总是报错”“到底怎么严谨地说某个基因属于某一群”的人。我会尽量把为什么要这么做的逻辑讲透而不是只给一段能跑的命令。1. 先把目标说清楚UMAP上标基因到底在解决什么问题1.1 三种最常见的应用场景很多教程会把“在UMAP上画基因表达”当成一个纯粹的画图问题但其实这一步背后至少要区分三种完全不同目的的场景。第一种是类群注释。我们无监督聚类之后得到一堆cluster每个cluster只是一堆编号并不知道它是什么细胞。这时候就需要拿一批经典的marker基因比如CD3D、CD14、MS4A1、NKG7这些逐个在UMAP上看它们的表达分布通过“这个基因只在0号群高表达”这样的信息来判断0号群是T细胞还是别的。这种情况下基因是已知的类群是未知的我们要做的是用基因给类群贴标签。第二种是候选基因定位。这种情况很常见比如做完差异分析或者功能筛选我们手里有一个“感兴趣的基因”它可能不是经典的marker甚至可能从来没被详细注释过。这时候我们想知道它富集在哪些细胞类群、是特异表达还是普通表达从而推测它的生物学功能。这种情况下类群可能已经注释清楚了基因是新的我们反过来用类群给基因“定位”。第三种是结果验证和逻辑检查。无论前两种情况做得多么顺利最后都要回到UMAP图上肉眼确认一遍。比如你注释某个cluster是CD8 T细胞那CD8A这个基因必须在这个cluster里有明显表达CD14这种单核细胞marker最好别在这里出现。这一步看起来朴素但几乎每一次我都会做因为它能快速暴露注释错位、基因名映射错误、数据归一化异常等一堆隐患。1.2 为什么偏偏是UMAP和t-SNE差在哪这一步有必要花点时间说清楚。很多人直接跟着教程用UMAP画图但并不知道为什么选它不选t-SNE。简单讲UMAP的全称是Uniform Manifold Approximation and Projection它的核心思想是在高维空间里构造一个带权重的k近邻图然后把这个图“压缩”到二维或三维空间里让相近的点靠拢、疏远的点分开。和t-SNE相比UMAP最大的优势是运行速度快、全局结构保留得更好。t-SNE非常擅长展示局部近邻关系但它对全局距离的还原很差不同簇之间的相对远近距离经常没什么意义。而UMAP在保留局部结构的同时还能大致维持簇与簇之间的全局关系这在看一群细胞里是否存在连续过渡状态时特别有用。另外UMAP算法本身更稳定相同参数下多次运行的结果差别更小虽然也不是完全确定性但配合固定随机种子基本可以复现。说到随机种子这其实是个很实际的坑。Scanpy里跑UMAP时默认就有random_state参数Seurat里则有seed.use选项我建议无论哪套流程都用固定种子否则你同一份数据今天画和明天画图的形状可能有肉眼可见的差异。要是不固定后续写文章插图时很麻烦。2. 动手之前必须搞懂的底层逻辑你画的表达值到底从哪来2.1 用哪个矩阵、哪个slot直接决定图对不对很多人第一次在UMAP上画基因表达代码是从教程里抄的比如Scanpy的sc.pl.umap(adata, colorCD3D)Seurat的FeaturePlot(pbmc, featuresCD3D)。命令跑通了就以为万事大吉但一旦看图发现“这基因怎么在所有细胞里都显色”或者“全是蓝色”往往就要回头检查数据到底存在哪一层。先说Scanpy。标准流程里adata.X在完成归一化和对数化之后存的是一个log1p(CPM/10k)的表达矩阵。但是很多流程接着会做高变基因筛选之后可能只剩两千个高变基因这时候如果你感兴趣的基因不在高变基因列表里直接color某个基因就会报错或者画不出来。因此成熟的习惯是在做任何基因筛选或降维之前先把全量数据存一份adata.raw adata很多教程从不解释这一行为什么重要。存下raw之后即使后续adata.X被重构为高变基因子集sc.pl.umap在遇到变量名不在当前var_names里时会自动尝试从adata.raw里取表达量。如果没有提前存raw后期再想画某个被过滤掉的基因就得重新读原始数据再走一遍归一化非常麻烦。再说Seurat。Seurat对象里至少有三个层counts是原始UMI计数data是归一化并log变换后的表达量scale.data是经过中心化和标准化(scale)后的残差。很多初学者会误以为应该用scale.data来画FeaturePlot因为它“看起来最标准”这是一个大坑。scale.data把每个基因的均值拉平成0方差拉平成1所以图上会出现大量负值蓝色红色对称分布肉眼很难判断这个基因到底是高表达还是低表达。FeaturePlot默认调用的是dataslot这个默认行为是对的不要去手动改成scale.data。一句话总结画UMAP上的基因表达分布应该用log归一化后的表达值不应该用scale后的值也不应该直接画原始UMI counts。用原始counts画会导致个别表达量极高的细胞把整个色标拉爆其他细胞全变成淡色。2.2 基因名大小写、物种与基因名映射的那些事单细胞数据标注基因名看起来只是一个小细节但实际踩坑的人非常多。人类基因的标准写法是全大写比如CD3D、CD14、MS4A1小鼠基因是首字母大写其余小写比如Cd3d、Cd14、Ms4a1。如果你用的是人和小鼠混合的数据或者参考基因组版本选错那很可能你的marker基因列表里有一半在var_names里根本找不到。我见过一种情况同学从一篇文献里抄了一组人源marker结果数据是小鼠的画图时全是KeyError。排查了很久才发现是物种问题。所以做这一步之前第一件事永远是确认你的数据物种可以用adata.var_names里随便挑几个已知基因做对照。另外不同版本基因注释可能会带.后缀或者-后缀比如CD3D.1这种这种通常是因为基因名去重时加了后缀要小心处理必要时统一去后缀再匹配。如果遇到基因名明明存在但sc.pl.umap报错可以用一行代码快速检查target_gene CD3D print(target_gene in adata.var_names) # 当前矩阵 print(target_gene in adata.raw.var_names) # raw矩阵如果两者都返回False那就是基因名本身不对或者物种不对。如果只存在于raw里说明你做了高变基因过滤但好在保存了raw画图时Scanpy仍能自动取到。2.3 颜色映射、色标范围与“全图蓝汪汪”的玄机每次讲到这一步我都要强调UMAP上基因表达图的颜色不等同于“有没有表达”。Scanpy和Seurat默认的色标范围都是从数据最小值到最大值如果这个基因的表达量跨度极大比如某个细胞表达量是8而绝大多数细胞都是0到1那么整个UMAP图上大部分点都会是浅色只有那个高表达细胞是深色视觉上会显得“这个基因好像只有极少数细胞表达”但其实丰度中等但分布广泛的细胞全都淹没在浅色里了。这种情况下有两种处理一种是把色标下限设为0上限设为一个合理的分位数比如95%分位数。Scanpy写法sc.pl.umap(adata, colorCD3D, vmin0, vmaxadata.X[:, adata.var_names.get_loc(CD3D)].max())或者更简单直接用vmax手动设一个经验值画完看效果再调整。Seurat里对应的是FeaturePlot(pbmc, featuresCD3D, colsc(lightgrey, red), min.cutoffq10, max.cutoffq95)。注意min.cutoff和max.cutoff用的是分位数/绝对值字符串能非常有效地避免离群值拉爆色标。还有一个容易被忽视的问题是点的大小和透明度。当细胞数量达到几万甚至几十万时UMAP上的点会叠在一起浅色表达区域容易被深色高表达点“遮盖”导致你低估一个基因的阳性细胞比例。我一般会在细胞数超过1万时适当调小size、调低alpha或者分面板放大局部区域观察。这些都是实操里很有用的细节。3. 完整实操从聚类结果到把marker基因回标到UMAP3.1 用一份公开PBMC数据把流程跑通理论讲再多不如直接拿数据跑一遍。这里我以单细胞领域最经典的一份公开PBMC数据为例这种数据在Scanpy或Seurat的官方示例里都能直接获取。先大概说一下数据预处理到聚类这一步的完整代码因为标注基因必须建立在质量过关的UMAP聚类结果之上。import scanpy as sc # 读取公开的PBMC数据集这一步得到的adata是原始counts矩阵 adata sc.read_10x_h5(pbmc_10k_v3_filtered_feature_bc_matrix.h5) # QC过滤低质量细胞和双细胞先清理掉 sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue) adata adata[adata.obs.pct_counts_mt 20, :] # 归一化和log变换 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) # 保存raw层这一步非常关键 adata.raw adata # 高变基因筛选 sc.pp.highly_variable_genes(adata, min_mean0.0125, max_mean3, min_disp0.5) adata adata[:, adata.var.highly_variable] # PCA、邻居图、UMAP、聚类 sc.tl.pca(adata, svd_solverarpack) sc.pp.neighbors(adata, n_neighbors10, n_pcs30) sc.tl.umap(adata, random_state42) sc.tl.leiden(adata, resolution0.5)这份数据里大概有九到十个主要细胞类群T细胞、B细胞、NK细胞、单核细胞、树突状细胞等。我刻意强调adata.raw adata这一步要在过滤高变基因之前执行就是为了后续能画任意基因。在拿到聚类结果后我先跑一遍差异基因看整体质量也顺便获得每个cluster的top markersc.tl.rank_genes_groups(adata, leiden, methodwilcoxon) sc.pl.rank_genes_groups(adata, n_genes10, shareyFalse)这一步输出的主要是排序后的差异基因名单它还不足以代替“用已知marker回标UMAP”的验证。差异基因和已知marker结合来看是判断cluster身份最可靠的做法。3.2 单基因检查与类群初判聚类完成后我先确认几个最经典的marker在UMAP上的表现。这里我会一次性画多个基因方便对比marker_genes [CD3D, CD14, MS4A1, NKG7, FCGR3A, IL7R, CD8A] sc.pl.umap(adata, colormarker_genes, ncols3, vmin0, vmax5, size20, alpha0.6)这样一张大图上六个面板同时展开一眼就能看到CD3D集中在哪几个clusterCD14又集中在哪几个。以PBMC数据为例通常结果是这样的CD3D在T细胞群高表达CD14在那群单核细胞里高表达MS4A1作为B细胞的经典marker只在B细胞群出现NKG7在NK细胞群亮起来。当你看到这些marker的表达区域与Leiden聚类的cluster边界高度吻合时聚类结果就基本可信了。到这里我还不会急着给每个cluster下最终结论。因为单靠眼睛看图判断“高表达”是有主观性的我通常会在看图的同一时间输出每个cluster的平均表达量和阳性细胞比例用数字验证肉眼印象。import pandas as pd # 按leiden分组计算每个cluster中该基因的平均表达量 df pd.DataFrame({ cluster: adata.obs[leiden].values, gene_expr: adata.raw[:, CD3D].X.toarray().flatten() if hasattr(adata.raw[:, CD3D].X, toarray) else adata.raw[:, CD3D].X.flatten() }) print(df.groupby(cluster)[gene_expr].mean().round(2))实际跑出来你可能会看到0号cluster的CD3D平均表达量是3.2而5号cluster是0.01这个数字比图上的颜色分布更可靠。但我还要提醒一句平均表达量低于0.1并不等于“没有表达”单细胞数据里大量dropout是正常现象所以还要结合阳性细胞比例来判断。3.3 多基因联动用共表达判断类群身份单个marker看多了之后会发现一个问题很多基因并不只在一个类群里表达。比如IL7R它既是CD4 T细胞的marker也在部分单核细胞里有低水平表达又比如FCGR3A它在CD16单核细胞和NK细胞里都有表达。如果只凭单个基因的颜色分布下结论很容易把CD16单核和NK混为一谈。这种情况下我的习惯是看两个基因的共表达情况。在Scanpy里可以直接给color传多个基因它会并排画多个UMAP方便对照边界。但我更常用的一招是给每个细胞算一个“共表达分数”也就是把两个基因的表达量相加或取平均然后画到UMAP上。比如我想判断某个cluster到底更像NK还是更像CD8 T我可以同时把CD8A和NKG7画出来观察这个cluster是哪个基因主导。要更精细地看共表达可以用sc.tl.score_genes给一组基因算模块分数sc.tl.score_genes(adata, gene_list[CD8A, NKG7], score_nameCD8_NKG7_score) sc.pl.umap(adata, colorCD8_NKG7_score, vmin0, size20)这个分数其实是对基因集表达量做的标准化综合比单独看两个基因更抗噪。在判断类似“这群细胞是不是细胞毒性T细胞”这类问题时比逐个基因翻图高效得多。3.4 从“看基因”到“标类群”给cluster起名字并输出最终图当我们通过几轮marker回标和共表达判断确认了每个cluster的身份之后最后一步是把Leiden聚类的数字编号替换成细胞类型名称再画最终版的UMAP标注图。Scanpy里可以这样做# 以PBMC示例为例替换为实际数据中对应的注释结果 cluster_to_celltype { 0: CD4 T cells, 1: CD14 Monocytes, 2: B cells, 3: CD8 T cells, 4: NK cells, 5: CD16 Monocytes, 6: Dendritic cells, 7: Megakaryocytes } adata.obs[cell_type] adata.obs[leiden].astype(str).map(cluster_to_celltype).astype(category) sc.pl.umap(adata, colorcell_type, size15)Seurat里的对应操作也很直接pbmc - RunUMAP(pbmc, reduction pca, dims 1:30) # 手动注释示例 new.cluster.ids - c(CD4 T cells, CD14 Monocytes, B cells, CD8 T cells, NK cells, CD16 Monocytes, Dendritic cells, Megakaryocytes) names(new.cluster.ids) - levels(pbmc) pbmc - RenameIdents(pbmc, new.cluster.ids) DimPlot(pbmc, reduction umap, label TRUE, pt.size 0.5, label.size 5)到这里“在UMAP图上标识感兴趣的基因所在的类群”这件事已经完成了一大半。如果要说最终交付的图通常是两个版本一个是以细胞类型为颜色标注的UMAP图备用于论文或报告另一个是把目标基因表达叠加在UMAP上、用来证明“这个基因确实富集在某类群”的证据图。两图配合使用效果最好。3.5 如果“感兴趣的基因”不是已知marker而是新发现的候选基因前面讲的都是已知marker如何注释类群但回到标题里的一个常见场景我们手里有一个通过差异分析或者GWAS结果锁定的新基因想知道它主要在哪些类群表达。这种需求的数据操作逻辑其实是反过来的——类群已经注释好了基因是新的。做法分为三步第一步先把基因画到UMAP上看整体分布第二步按注释好的细胞类型逐群计算这个基因的阳性细胞占比和平均表达量第三步用小提琴图或点图展示这个基因在不同类群的表达丰度差异。import scanpy as sc gene 某个感兴趣的基因 # 第一步看UMAP分布 sc.pl.umap(adata, colorgene, vmin0, size15) # 第二步按细胞类型统计表达量 sc.pl.violin(adata, keysgene, groupbycell_type, rotation90, stripplotFalse)如果这个基因在某个类群中的表达明显高于其他类群那就可以初步说这个基因富集在XX类群。但要注意在PBMC这类复杂组织里很多基因的表达谱并不那么“干净”特别是一些管家基因或者应激相关基因往往在多个类群里同时表达这时候不要急着下结论建议进一步分析这类群的亚群结构或者做RNA速度、功能富集来佐证。4. 高频翻车现场我在这步踩过的坑和排查思路4.1 基因名大小写与物种错配这个问题值得再说一次因为我在团队里几乎每个月都会看到一次。人类和小鼠基因名写法的差异在前面已经提过但还有一个更隐蔽的坑如果你用的数据来自某些公共数据库基因名可能已经被转换过比如Ensembl ID、Symbol混在一起或者基因名带了版本号。我在跑一个数据集的时候遇到过var_names里一半是symbol、一半是Ensembl ID的情况其实是早期下载数据时没有统一注释版本。遇到这种情况最稳妥的办法是在做表达矩阵时就统一好基因注释比如用scanpy的sc.pp.annotate或者直接查一下参考注释文件。如果数据已经读进来了那就写一个简单的映射字典把Ensembl ID批量换成symbol# 假设gene_map是从注释文件里读出的字典 adata.var_names [gene_map.get(i, i) for i in adata.var_names]改完之后记得检查一遍确认没有重复的基因名否则后面画图会报“duplicate var_names”的错误。4.2 全图蓝汪汪一片是不是做错了这类问题几乎每周都有人问。UMAP上某个基因显示全蓝最直接的一个可能是这个基因在几乎所有细胞里都不表达或者表达量极低比如你画一个B细胞特有的marker但这份数据里恰好没有B细胞或B细胞极少那图当然全蓝。这种情况并不一定是代码问题而是生物学事实。另一种可能是色标范围被离群值拉爆了。前面提到的vmin/vmax调整就是为这个准备的。实际操作时我会先大致看下基因的表达量分布import numpy as np expr adata.raw[:, CD3D].X.toarray().flatten() print(np.percentile(expr, [0, 10, 50, 90, 99, 100]))如果中位数是0而99分位数也很低说明基因确实低表达如果99分位数和最大值差异巨大说明有极端离群值。根据分位数输出结果再设vmax就合理多了。还有一点容易忽略有些流程在跑完sc.pp.scale之后adata.X会被替换成scale之后的矩阵而adata.raw仍然是log1p的归一化矩阵。如果你的代码里不小心用了adata.X而不是adata.raw的值画图得到的UMAP图看起来就会非常奇怪基因表达负值多、色块不连续。这就是为什么我一直强调画基因表达图要认准表达矩阵存放在哪里。4.3 同一份数据UMAP图每次跑都不一样UMAP算法本身带有随机性特别是初始化阶段。如果你不固定随机种子今天跑出来的cluster边界和明天跑出来的可能不同虽然大致类群结构类似但具体编号可能对不上。这一点在标注类群、复现结果时非常容易造成混乱。Scanpy里统一固定种子其实很简单sc.tl.umap(adata, random_state42)如果后面想微调图效果可能会用到sc.tl.umap的n_epochs参数但常规分析建议保持默认只固定种子。Seurat里的对应做法是set.seed(42) pbmc - RunUMAP(pbmc, reduction pca, dims 1:30, seed.use 42)另外还有一种情况会让UMAP图“不稳定”数据里混杂了非常强的批次效应或技术差异。这种情况下不同批次细胞会各自聚成一坨而不是按照细胞类型掺杂在一起UMAP图形也会随参数剧烈变化。如果遇到这种情况先解决批次效应比如用数据整合的方法跑完再重新聚类、重新UMAP否则后面所有基于UMAP的视觉判断都不可靠。4.4 一个基因在多个类群里都表达是不是注释错了这很常见不一定代表注释错误。以PBMC为例CD3D在CD4 T和CD8 T细胞里都高表达NKG7在NK细胞里高表达同时也在部分细胞毒性T细胞里表达。这是正常的生物学现象因为细胞类型之间本来就不是完全离散的而是一条连续的谱系。关键不在于“这个基因是否只在某群表达”而在于“这个基因在哪群表达最高、在哪群几乎不表达”。所以每次遇到多群阳性时我都会回到定量分析计算每个类群的阳性细胞比例和平均表达量用数字补充图上的视觉判断。如果无论如何都无法用已有marker把某个cluster和新注释对应起来那就应该重新审视这个cluster是不是一个之前没有预期到的新细胞亚型而不是急着硬套已知类型。这也是单细胞分析里最有趣的环节之一——你可能会发现一个新类群。4.5 常见问题速查表我把高频问题和对应处理思路整理成一张表方便遇到问题时快速对照。现象可能原因处理方式基因名报错找不到物种/大小写错误基因被过滤检查物种检查raw.var_names图全蓝/全浅色基因表达过低色标范围不合理调整vmin/vmax查表达量分位数图中出现大片负值用了scale.data或scale后的矩阵确认使用log归一化矩阵或adata.rawUMAP每次图都不同未固定随机种子批次效应强固定random_state先做批次整合基因在多个群里阳性基因本身谱系分布宽看定量表达比例不要单凭视觉定论加了细胞类型注释后图空白obs映射字典键不对检查leiden的category和map字典key5. 把这步做顺的几条实战心得5.1 自动化批量扫marker别一个个手点当你要注释一批细胞类型时逐个基因画UMAP效率很低。我更推荐的方式是准备一个marker基因列表一次性批量出图然后快速翻看结果。Scanpy里可以直接用sc.pl.dotplot或者sc.pl.matrixplot把一组marker在不同cluster里的表达量浓缩成一张点图或热图。点图里点的大小和颜色分别代表阳性细胞比例和平均表达量一眼就能看出哪个cluster富集了哪类marker。sc.pl.dotplot(adata, var_namesmarker_genes, groupbyleiden, standard_scalevar, figsize(10, 6))这张图配合UMAP上的单基因表达图一起看注释效率会高很多。但点图并不能完全替代在UMAP上直接看分布因为有些marker的表达模式在点图里“看起来差不多”实际在UMAP空间里有明显的梯度或亚群结构只有回到UMAP图上才能发现。所以我的流程是先用点图快速锁定方向再用UMAP单基因图精确定位边界。5.2 给图增加“图层信息”大小、透明度、分面在细胞数很大的数据集里UMAP图细节很容易被点重叠掩盖。常规的建议是size控制在10到30之间alpha在0.5左右如果你的数据量超过5万细胞可以考虑分成几个panel局部放大看。Seurat里可以用DimPlot(pbmc, split.by sample)来按样本分面查看这在确认不同样本来源的细胞是否混合均匀时非常有用。Scanpy里对应的参数是sc.pl.umap(adata, colorleiden, ncols2)。还有一个小技巧如果某个cluster的细胞数太多、盖住了另一个簇的边界可以用sc.pl.umap的连接边参数edgesTrue看UMAP图内部结构。连接边其实就是k近邻图的可视化能帮你判断哪些细胞彼此距离更近、是否有清晰的间隙把两群分开。注意edges只适合在小数据集或局部区域打开否则会非常占绘制资源。5.3 把注释结果从“视觉印象”变成“可交付结论”很多同学停在“图看得差不多”就结束了但如果这个结果要用于报告或论文建议再多做两步。第一步是把注释结果写回obs并单独导出一张带细胞类型标注的UMAP图分辨率至少300dpi。第二步是补充一个marker表达的小提琴图或者点图表格作为类群注释的定量证据。这一步不是强制要求但审稿人或者导师大概率会问“你怎么确定这个cluster是B细胞”有一张数字支撑的图比口说更有说服力。最后再分享一个我自己的小习惯。每次标注完类群我不会马上改obs里的注释列而是先在旁边新增一列比如cell_type_manual保留原始的leiden编号结果。这样后续如果注释方案有调整或者需要回溯原始cluster做亚群分析就不需要重新聚类重新跑一遍了。这个习惯看起来不起眼但在项目后期反复调整注释时帮我节省了大量时间。
