我前后帮几个做单细胞的朋友跑过官网那套分析流程最大的体会是前面数据读取、质控、标准化、高变基因筛选这几步大家照着示例抄基本不会错真正开始出问题的是从邻接图往下——聚类数量怎么定、UMAP 图为什么每次跑出来都长得不一样、差异基因表里那一堆数字到底该信哪个这些问题官网示例往往一句“just run it”就带过去了。这篇接着前两篇的进度把 scanpy 单细胞分析官网示例的后半段拆开讲neighbors建图、UMAP 降维、Leiden 聚类、rank_genes_groups差异表达、以及最后的细胞类型注释。这些步骤在代码层面只有五六行但每一步背后的参数取舍、结果解读和踩坑点足够让人反复折腾一整周。我尽量把每个“为什么”都讲透让刚上手的人知道怎么跑通也让跑过几轮的人知道哪里容易翻车。1. 邻接图构建neighbors 这一步最短坑却最密1.1 从 PCA 结果到 KNN 图中间到底算了几件事很多人以为sc.pp.neighbors(adata)就是“算个邻居关系”点完就跑根本没意识到它内部至少干了三件事。第一件是在 PCA 降维后的低维空间里对每个细胞找出最近的若干邻居构成 KNN 图第二件是把 KNN 图转换成带权重的连接概率让“邻居的邻居”之间也能建立弱连接这一步参考的是 UMAP 的模糊单纯集思想第三件是把每个细胞的局部邻域关系拼成全局的可遍历图结构。这三件事做完结果才会分别落在adata.obsp[distances]和adata.obsp[connectivities]两个稀疏矩阵里。为什么要强调这个结构因为后面几乎所有操作都依赖它。sc.tl.umap读的是这个图sc.tl.leiden读的也是这个图甚至连一些轨迹推断、邻域富集分析的入参都是它。你如果在这一步参数调得不对后面所有下游结果会一起偏而且偏的方向很难察觉——UMAP 图照样画得出来聚类照样有编号只是那些编号对应的生物学意义可能已经乱了。我在第一次认真看这步的时候就是被“代码只写了一行、结果却影响全局”这件事给噎住的。还有一个容易忽略的点neighbors默认会用已计算的所有 PCA 成分。如果前面sc.tl.pca没有显式指定n_comps默认会算到 50 个成分或者数据维度允许的最大值然后neighbors会全部拿来用。这意味着前面 PCA 里那些解释方差极低的成分也被当成有效信息塞进了距离计算里。在 PBMC 这种细胞类型差异明显的样本中这通常问题不大但在细胞类型差异细微、或者批次噪声比较重的数据里噪声主成分会污染邻接图聚类会出现大量碎片。所以我现在的习惯是显式写n_pcs40甚至在 PCA 前先用sc.pp.highly_variable_genes把基因卡严一点从源头减少噪声。1.2 n_neighbors 和 n_pcs 这对参数该怎么配合n_neighbors控制的是每个细胞在图里连多少个邻居。官网示例里用的是 10默认值是 15这两个数看起来差不多实际效果差异不小。n_neighbors小图更局部能刻画细微的亚群结构但容易把连续过渡的细胞群切碎成多个小团n_neighbors大图更平滑主要群体结构清晰但小群体容易被吞进大群体里。这不是玄学而是 KNN 图本身的局部性与全局性权衡。n_pcs控制的是用多少个主成分来算距离。用多了会把噪声带进来用少了会丢失真实的细胞状态差异。经验上n_pcs大致取到 PCA 碎石图里“肘部”之后一点的位置比较稳妥。以常用的那套外周血单个核细胞示例数据为例前 10 个主成分基本覆盖了主要的细胞类型差异取 30 到 40 是一个安全区间。取 50 也不是不行但对这类数据来说增益有限反而增加计算量。这两个参数不是独立的。n_neighbors小的时候n_pcs取大一点反而更容易让图变得不稳定因为每个细胞在更多维空间里找近邻找到的邻居更“散”。我做过一组对照固定n_pcs40把n_neighbors从 5 到 30 扫一遍聚类数量会出现明显波动反过来固定n_neighbors15扫n_pcs波动相对温和。所以调参顺序建议是先把n_pcs定在一个合理范围再微调n_neighbors不要两头同时动否则根本分不清结果变化是谁造成的。1.3 图对象的结构与三个高频报错neighbors跑完之后adata里多了三样东西adata.uns[neighbors]存参数和元信息adata.obsp[distances]存距离矩阵adata.obsp[connectivities]存连接权重矩阵。前两个是稀疏矩阵用print看会显示形状和稀疏度。理解这个结构对后面排查问题很有帮助如果 UMAP 结果明显不对劲先查connectivities是不是几乎全零或者密度异常基本能定位到是建图环节的问题。高频报错里第一个是“数据还没降维你就来建图”典型提示是让你先跑 PCA。这个好解决。第二个是n_pcs超过实际计算的主成分数比如 PCA 只算了 30 个成分你却传了 50会直接报越界。这个坑在拼接多个数据集之后特别常见因为不同批次算的 PCA 维度可能不一致。第三个是稀疏矩阵内存爆掉尤其当细胞数上到十几万、n_neighbors又设得比较大的时候connectivities的稀疏度下降内存占用会快速上涨。我一般会在建图前用adata.n_obs和adata.n_vars估一下规模超过十万细胞就考虑先做子集抽样或者把n_neighbors压到 10 以内。2. UMAP 降维图画得漂亮不等于读得正确2.1 UMAP 凭什么默认排在 t-SNE 前面官网示例在降维可视化这一步只用了 UMAP没再跑 t-SNE。这不是随意的选择。UMAP 在保留全局结构上比 t-SNE 更稳一些细胞群之间的相对位置有一定参考价值t-SNE 更强调局部邻域全局距离基本没有意义。另外 UMAP 的计算速度在细胞数上千到上万这个量级上明显更快配合n_jobs并行跑几分钟就能出图。对刚上手的人来说UMAP 还有一个隐性优势它的结果对随机种子的敏感度相对可控不像 t-SNE 那样稍微动一下种子图就面目全非。但“全局结构相对保留”这句话经常被过度解读。UMAP 图里两个群离得远不代表它们在表达谱上也离得远两个群挨得近也不代表它们生物学上更相似。它保留的是拓扑意义上的邻接关系不是欧氏距离。我见过不少人拿着 UMAP 图说“这两个群距离这么远肯定是不同来源”结果一查差异基因发现只是同一个大类下的两个状态。这个问题在跨样本整合之后尤其突出因为整合本身就会重塑低维空间。所以我的用法是UMAP 图主要用来“看结构”比如有没有明显孤立的群、群内是否均匀、有没有连续的过渡带而不是用来“量距离”。真要比较群之间的关系去翻差异表达结果或者做相关性分析比盯着图看靠谱得多。2.2 随机性与并行计算带来的可重复性问题UMAP 这一步最让人头疼的是“同样的代码换台机器图就不一样”。原因有两个随机数初始化和并行计算的非确定性。sc.tl.umap提供了random_state参数默认是 0理论上固定种子就能复现。但如果n_jobs设成大于 1底层某些并行实现会引入非确定性导致即使种子固定每次的嵌入结果仍有细微差异。这个差异在细胞群交界处会被放大看起来就像两个群“融合”或者“分开”了。我的处理方式是正式出图、要写进报告的那一次把n_jobs显式设成 1random_state固定成一个写死的数跑完立刻把嵌入坐标存下来。日常探索阶段可以开并行图快但别拿探索阶段的图去和正式图对比否则你会怀疑人生。还有一点UMAP 的min_dist和spread也会影响图的“紧致程度”min_dist越小群内越紧凑群间空隙越大视觉上更“干净”但这种干净是参数调出来的不是数据本身的信号别被它骗了。下面是一段我常用的固定写法把可复现性相关参数都显式写出来import scanpy as sc sc.tl.umap( adata, min_dist0.3, spread1.0, n_components2, random_state42, n_jobs1, ) adata.obsm[X_umap_final] adata.obsm[X_umap].copy()把最终嵌入单独存一份到obsm里是个小技巧。后面如果又跑了一次 UMAP 做对比X_umap会被覆盖但X_umap_final还在随时能拿出来画图不会把之前的结论弄丢。2.3 读图时最容易犯的三个误判第一个误判是把“群的大小”当成“细胞数量”。UMAP 里群的面积和细胞数量没有严格对应关系一个只有几十个细胞的稀有群可能因为分布分散而看起来比一个上千细胞的主群还大。判断细胞数量要去看聚类结果的计数表不能看图估。第二个误判是把“群间距离”当成“分化关系”。前面说过全局距离不可靠但很多人还是会下意识地认为相邻的群是“分化来的”。要判断分化关系得用拟时序或者 RNA 速率这类方法UMAP 图只提供线索不提供结论。第三个误判是忽略“群内异质性”。UMAP 把一个群画成一个紧凑的团容易让人以为这个群内部是同质的。实际上很多群内部还有亚结构只是分辨率参数或者降维参数把它压平了。要检查这一点可以把这个群的细胞单独抽出来重新做一次 PCA 和 UMAP或者看这个群内部差异基因的分布。我在一次分析里就吃过这个亏把一个明显包含两个亚群的 T 细胞群当成一个整体去注释结果差异表达分析出来的 marker 基因互相打架排查了半天才发现是聚类粒度过粗。3. Leiden 聚类resolution 才是真正决定聚类数量的旋钮3.1 为什么官网示例换了 Leiden 而不是继续用 Louvain早期版本的示例用的是 Louvain后来换成了 Leiden。这个更换不是跟风Leiden 算法在 Louvain 的基础上加了连通性保证能避免产生“内部不连通”的社区。Louvain 有时候会把两个实际上没有边相连的子图分到同一个社区里这种社区在生物学上很难解释因为它根本不是连通的细胞群。Leiden 通过一个精修步骤把这类问题处理掉了结果更干净。用 Leiden 需要在环境里装leidenalg和igraph两个包。这两个包在部分环境里装起来会有点折腾尤其是 Windows 上版本不匹配时报错信息还不太友好。我的建议是直接用 conda 装别硬用 pip 去凑版本。装好之后sc.tl.leiden调用本身很轻一两行就完事。新版 scanpy 还提供了flavorigraph参数底层换用 igraph 自带的实现速度更快社区划分结果和leidenalg版本略有差异但整体一致。我一般默认用flavorigraph除非要做严格的历史结果对比。3.2 resolution 扫描我一般怎么选合理的聚类粒度resolution是 Leiden 里最重要的参数没有之一。它控制社区划分的粒度值越大社区越多、越细值越小社区越少、越粗。官网示例用的是 1.0这个值在 PBMC 这类数据上大概能分出七八个群对应主要的免疫细胞类型。但 1.0 不是一个通用最优值换一套数据可能就完全不够用。我通常的做法是做一个 resolution 扫描把 0.2 到 1.5 之间取五到六个值每个值跑一遍 Leiden记录聚类数量和各大群的 marker 基因分布然后挑一个“聚类数量合理、群内 marker 一致、群间 marker 区分明显”的值。判断标准不是聚类数量本身而是每个群能不能被一组干净的 marker 基因支持。如果某个 resolution 下出现大量只有几十个细胞、marker 杂乱的群说明粒度过细如果出现一个群里混着好几种 marker 谱系说明粒度过粗。下面这个扫描写法我用了很多次跑完直接看表格对比import pandas as pd res_list [0.2, 0.4, 0.6, 0.8, 1.0, 1.2] records [] for res in res_list: key fleiden_{res} sc.tl.leiden(adata, resolutionres, key_addedkey, random_state0) n_clusters adata.obs[key].nunique() sizes adata.obs[key].value_counts() records.append({ resolution: res, n_clusters: n_clusters, min_size: sizes.min(), max_size: sizes.max(), }) df_res pd.DataFrame(records) print(df_res)扫描的好处是你能直观看到在哪个 resolution 区间聚类数量发生跳变。跳变点往往对应着数据里一个真实的结构层次比如从“大类”跳到“亚类”。这个信息比单看一个 resolution 的结果有价值得多。3.3 聚类编号的随机性与结果对齐Leiden 的聚类编号0、1、2……本身没有生物学含义而且每次跑编号和群的对应关系可能变化尤其是换了随机种子之后。这意味着你不能拿上一次分析里“cluster 3 是 T 细胞”的结论直接套到新结果上必须重新看 marker 基因。为了减少这种困扰我的习惯是一旦确定了一个最终使用的 resolution把random_state固定下来key_added用一个明确的名字比如leiden_final跑完就把这个标签存进obs的一个独立列里后续所有注释、差异分析都基于这一列。这样即使后面又跑了别的聚类做对比也不会互相干扰。另外如果用多个 resolution 的结果做对比图记得在图上把每个 resolution 的编号都标清楚不然读者根本分不清哪列是哪列。还有一个细节sc.tl.leiden默认会把结果写进adata.obs[leiden]如果你之前已经有一列leiden会被覆盖。这个行为在写脚本时很容易埋雷。我现在的做法是所有聚类结果都用key_added显式命名从不依赖默认列名养成习惯之后脚本的可维护性会好很多。4. rank_genes_groups差异表达分析的方法选择与结果解读4.1 t-test、wilcoxon、logreg 三种方法的取舍sc.tl.rank_genes_groups提供了好几种差异检验方法示例里默认用的是 Wilcoxon 秩和检验。这三种方法各有适用场景选错了不会报错但结果的信噪比会差很多。Wilcoxon 是非参数检验不假设表达值服从正态分布对单细胞数据里大量的零值和偏态分布更稳健是默认推荐。t-test 假设近似正态在细胞数多、表达分布比较连续的时候速度快、结果也还行但对稀疏数据不够友好容易给出偏乐观的显著性。logreg 用的是逻辑回归输出的不是传统意义上的 p 值而是每个基因对分类的贡献程度适合做特征筛选不适合做“这个基因是否显著差异表达”的统计判断。我的默认选择是 Wilcoxon除非细胞数特别大、跑起来太慢才考虑用 t-test 加速。t-test_overestim_var这个变体在组间细胞数不平衡时表现比普通 t-test 好一些但也只是相对好。不管用哪种方法都要记得数据必须先做过标准化和 log 变换否则算出来的倍数变化和 p 值都不可靠。4.2 校正 p 值、logFC 与表达比例三联指标一起看差异表达结果表里通常包含几列names、scores、logfoldchanges、pvals、pvals_adj。很多人只看pvals_adj排序挑前几十个基因就完事这是个常见误区。校正后的 p 值只说明“这个差异不太可能是随机产生的”不说明差异有多大、也不说明这个基因在多少细胞里表达。我的筛选习惯是三联指标一起看pvals_adj小于 0.05 作为显著性门槛logfoldchanges绝对值大于 0.5 作为效应量门槛再结合该基因在目标群里的表达细胞比例可以用sc.get.rank_genes_groups_df拿到结果后自己算或者用sc.pl.rank_genes_groups_dotplot看表达比例。三个条件同时满足的基因才是既有统计显著性、又有实际区分能力的 marker。还有一个细节rank_genes_groups默认是“一个群对其余所有群”的比较这种比较方式在群数量多的时候会把信号稀释掉。如果要找某个群相对某个特定群的差异基因用groups和reference参数指定对比组结果会干净很多。比如想找群 0 相对群 1 的差异基因就把groups0、reference1传进去比默认的全群比较更有针对性。4.3 把差异基因结果导出成可交付表格分析做完结果总得交出去。adata.uns[rank_genes_groups]是一个嵌套字典直接打印很乱用sc.get.rank_genes_groups_df转成 DataFrame 最方便df_markers sc.get.rank_genes_groups_df(adata, groupNone) df_markers.to_csv(markers_all_groups.csv, indexFalse) df_top ( df_markers .sort_values([group, pvals_adj]) .groupby(group) .head(50) ) df_top.to_csv(markers_top50_per_group.csv, indexFalse)导出的时候建议同时存一份全量和一份 top N。全量留着备查top N 用于汇报和画图。列名最好在导出前改成人能看懂的名字比如把pvals_adj改成p_adj、logfoldchanges改成log2fc这样交给不熟悉 scanpy 的人也能看懂。还有一点导出的 CSV 里基因名可能包含特殊字符某些表格软件打开会乱码建议统一用 UTF-8 编码存必要时再加个 BOM。5. 细胞类型注释把聚类编号翻译成生物学标签5.1 marker 基因怎么选、怎么验证聚类跑完得到的是一堆数字编号注释就是给这些编号贴上生物学标签。这一步没有捷径只能靠 marker 基因。经典的外周血免疫细胞 marker 大致有这么几组T 细胞看CD3D、CD3EB 细胞看MS4A1、CD79A自然杀伤细胞看GNLY、NKG7单核细胞看LYZ、CD14树突状细胞看FCER1A、CST3血小板看PPBP。这些基因不是随便选的它们各自对应明确的细胞谱系标志在多数外周血样本里表现稳定。选 marker 的原则是“一组基因互相印证”而不是“一个基因定生死”。单个基因可能因为技术噪声、批次效应或者表达波动出现假阳性或假阴性但一组功能相关的基因如果同时在一个群里高表达可信度就高得多。我一般每个大类准备三到五个候选 marker交叉验证之后再下结论。验证的方式有两种。一种是看 dotplot横轴是群、纵轴是基因点的大小表示表达细胞比例、颜色深浅表示平均表达量一眼就能看出哪个群对应哪组 marker。另一种是看每个群的 top 差异基因列表如果排在前面的基因里有已知 marker基本就能确认。两种方式结合注释的准确率会高很多。5.2 dotplot 和 stacked_violin 在注释里的分工sc.pl.dotplot适合看“群和基因的对应关系”信息密度高适合放在正文里做总览。sc.pl.stacked_violin适合看“表达分布的形状”能看出某个基因在一个群里是双峰、单峰还是几乎不表达这对判断群内是否存在亚结构很有帮助。比如一个 T 细胞群里CD8A呈现双峰分布提示这个群可能混着 CD4 和 CD8 两个亚群需要拆开再看。我通常先画 dotplot 做初筛锁定几个可疑的群再对这几个群画 stacked_violin 确认表达分布最后用sc.pl.umap把关键 marker 的表达叠到 UMAP 上看空间分布是否和聚类边界一致。三步走完注释结论基本就稳了。下面是一个常用的可视化组合marker_dict { T cells: [CD3D, CD3E], B cells: [MS4A1, CD79A], NK cells: [GNLY, NKG7], Monocytes: [LYZ, CD14], Dendritic: [FCER1A, CST3], Platelets: [PPBP], } sc.pl.dotplot(adata, marker_dict, groupbyleiden_final, standard_scalevar) sc.pl.stacked_violin(adata, marker_dict, groupbyleiden_final, rotation90)standard_scalevar让每个基因的表达量按自身范围缩放避免高表达基因把低表达基因的信号压掉。这个小参数对读图体验影响很大建议默认加上。5.3 注释结果的保存与复现清单注释做完把标签写回adata.obs是最后一步。可以直接用adata.obs[cell_type] adata.obs[leiden_final].map(label_map)也可以用rename_categories。我更倾向于新建一列cell_type保留原始聚类编号这样后面如果要回头检查某个群的注释依据原始编号还在不会丢线索。保存方面.h5ad格式是首选能把表达矩阵、降维结果、聚类标签、差异分析结果一起存下来下次直接sc.read_h5ad读回来就能继续分析。存之前记得检查一下adata.obs的列名有没有冲突、adata.uns里有没有存了特别大的临时对象。复现清单我一般会写这么几条数据版本、scanpy 版本、关键参数n_neighbors、n_pcs、resolution、random_state、marker 基因列表、注释映射表。这几样齐了换个人、换台机器也能把结果复现出来。最后分享一个我踩过多次坑之后养成的习惯每完成一个关键步骤就把adata另存一个带步骤名的副本比如adata_after_neighbors.h5ad、adata_after_leiden.h5ad。单细胞分析流程长、参数多中间任何一步想回头重跑都很常见有这些副本在手回退成本几乎为零。磁盘空间换来的安心远比事后从头再跑一遍划算。
