做单细胞数据分析的人多半绕不开 Seurat。我前面两篇教程把 R 语言里 Seurat 的完整流程写过一遍从数据读入、质控、标准化到聚类和 marker 基因鉴定一套跑下来确实顺手。但这几年我自己的项目里出现了一个很现实的情况样本量越来越大、数据动辄几十万细胞R 的内存开销开始吃紧同时分析流程里需要跟 Python 生态里的深度学习模型、可视化库、scanpy 之外的其它工具做配合用 R 做完整流程就感觉两头受堵。于是我开始把整套单细胞转录组测序数据分析流程往 python 的 scanpy 上迁。这篇教程就是给那些已经在用或者准备用 scanpy 的朋友准备的尤其是从 Seurat 转过来的那一批人。我不打算把官方文档抄一遍而是站在一个实际跑过大量样本的角度把 scanpy 的完整分析流程拆开来讲环境怎么搭、AnnData 对象怎么理解、每一步参数为什么要这么设、跟 Seurat 的哪些函数能对上号以及我踩过的那些坑。内容覆盖从 10x 原始数据到 marker 基因输出的全流程适合已经会一点 Python、但还没系统跑过单细胞分析的人参考。1. 从 Seurat 转向 scanpy 的真实理由与生态对比1.1 为什么你该考虑在 Python 里做单细胞分析先说结论我不是劝你放弃 Seurat。如果你是纯粹做单细胞分析、不做太多定制化开发Seurat 依然是很成熟的选择文档清楚、社区活跃、画图也漂亮。但当你的任务开始变复杂的时候Python 这套生态的优势就出来了。第一个理由是规模。Seurat 的对象设计本身是没有问题的但它默认把大部分数据加载进内存表达矩阵加上降维结果、邻接图这些一个十万细胞的数据集跑完 UMAP 和聚类内存随便就吃掉二三十个 G。scanpy 底层用 AnnData 对象基于 h5py 和 scipy 稀疏矩阵能够更有效地存储和处理大规模数据配合 out-of-core 读取方式处理百万级细胞也不是稀罕事。第二个理由是生态联动。单细胞数据分析从来不只停在聚类和 marker现在大家跑到最后都要接拟时序、细胞通讯、转录因子分析。拟时序常用的 Monocle 在 R 里但 Python 里也有 scVeloRNA velocity、CellRank 这些更往动力学方向走的工具。你要做深度学习批处理、做细胞类型注释Python 这边直接可以调用不需要跨语言搬运数据。第三个理由很实际我的合作方越来越多人交出来的数据预处理结果都是 h5ad 格式跟 Python 的 scanpy 无缝衔接。如果在 R 里跑完 Seurat想转出去给别人做下游还要过一层格式转换麻烦不说还有可能丢失一部分元数据信息。1.2 scanpy 与 Seurat 的定位差异scanpy 的设计哲学和 Seurat 不太一样。Seurat 是把一整套流程封装得比较完整你按着标准流程走UI 式的体验默认参数基本能给出一个可接受的结果。scanpy 更像是一个工具箱它提供的是模块化的函数每一步你都可以单独调用也可以自由组合。这种差异在调试的时候体现得很明显。用 Seurat 的时候如果一段流程报错你往往需要去理解它内部完整的 S4 对象结构看它是怎么存储和传递数据的。而 scanpy 的函数大多是接收 AnnData、返回修改后的 AnnData你对每一步之间数据是怎么流动的心里是很清楚的。还有一个差异在可视化上。Seurat 的画图函数内置度很高各种 DimPlot、FeaturePlot、VlnPlot 都给你配好了。scanpy 的绘图功能其实也不弱而且它是基于 matplotlib 的所以你对图表的控制力更强——你想调整坐标轴、字体、颜色、布局直接在 matplotlib 层面上做就好不用去破解 Seurat 的画图函数参数。但也因为这种灵活性你有时候需要自己多写几行代码来调样式不能像 Seurat 那样开箱即用。2. 环境搭建与安装避坑2.1 用 conda 建环境比直接用 pip 省心十倍很多刚开始玩 python 单细胞分析的朋友第一步就栽在环境上。直接在系统 Python 里 pip install scanpy多半会把依赖搞乱尤其是 numpy、scipy、pandas 这些科学计算库的版本稍有不慎就会冲突。我自己用的方式是 conda 建独立环境一步到位干净隔离。创建环境的命令很常规但有几个小地方要留意conda create -n scanpy python3.9 -y conda activate scanpy我推荐用 Python 3.9 而不是最新版本。不是说我保守而是 scanpy 生态里很多依赖库在最新 Python 版本下可能会出现预编译包缺失的问题conda 要现场编译的话会很痛苦。3.9 是当前兼容性最好的版本大部分生物信息学工具都能在上面稳定运行。接着装核心包conda install -c conda-forge scanpy python-igraph leidenalg -y这里有一个很关键的细节leidenalg 和 python-igraph 一定要一起装。leiden 聚类算法是 scanpy 做聚类时的默认选择比早期的 louvain 算法在社区划分质量上更好。但 leidenalg 的安装经常出问题尤其是当你只通过 pip 安装的时候很可能会遇到编译错误。用 conda-forge 渠道安装它会自动帮你匹配好 igraph 和 leidenalg 的版本省去很多麻烦。安装完之后建议顺手装几个下游常用包pip install scrublet pip install scvi-toolsscrublet 是用来预测双细胞的scvi-tools 是深度学习的批次整合工具。当然这些不是必须要装的看你项目需求来定。装完之后用python -c import scanpy as sc; print(sc.__version__)验证一下能出来版本号就说明环境没问题。2.2 关于国内镜像源的取舍这一步纯粹是实际体验问题。conda 和 pip 默认的源在下载速度上确实不稳定国内的朋友经常会遇到下载到一半超时的情况。我可以明确地跟你说配置国内镜像源是合规且非常常见的做法能显著提升安装效率。比如 pip 源我长期用清华的 PyPI 镜像pip config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simpleconda 源也是一样在~/.condarc里配置一下清华的 conda-forge 和 defaults 镜像。但要注意不要同时在.condarc里混合太多不同来源的 channel不然 conda 在解析依赖时会非常慢甚至出现 channel 优先级冲突导致装不上包的问题。还有一个小建议装包的时候尽量用mamba替代condamamba 用的是 C 写的求解器依赖解析速度比 conda 快很多尤其是环境里要装的包比较多的时候差距是数量级的。3. scanpy 完整分析流程实操3.1 数据读入与 AnnData 对象初识拿到单细胞转录组测序数据之后第一步肯定是把数据读进来。scanpy 支持很多格式但我实际项目里碰到最多的还是 10x Genomics 的 Cell Ranger 输出。如果是 Cell Ranger 的过滤矩阵输出目录通常包含 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz 三个文件读取方式如下import scanpy as sc import numpy as np adata sc.read_10x_mtx( data/filtered_feature_bc_matrix/, var_namesgene_symbols, make_uniqueTrue )如果是单个的 h5 文件就用sc.read_10x_h5()。读取之后adata 就是一个 AnnData 对象它的核心结构是adata.X存的是细胞×基因的表达矩阵。这个矩阵在 Python 里是稀疏矩阵格式这也是为什么 scanpy 能处理大规模数据的关键之一——大部分基因在大部分细胞里表达量是 0稀疏矩阵只存非零值内存占用可以压缩到原来的几十分之一。拿到数据之后第一步不要急着跑分析先看一眼数据的基本面貌print(adata.shape) print(adata.var.head()) print(adata.obs.head())adata.obs存的是细胞层面的元数据adata.var存的是基因层面的元数据。后面每一步分析的结果比如质控指标、聚类标签、UMAP 坐标都会存到这两个 DataFrame 里。理解 AnnData 这个结构后面的逻辑就顺了。3.2 质控过滤线粒体基因与双细胞质控是整个流程里最影响结果的一步这一枪打得不准后面再怎么调参都白搭。我一般会做三件事过滤低质量细胞、过滤低表达基因、预测并去除双细胞。先计算质控指标adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue)这一步会在adata.obs里生成n_genes_by_counts、total_counts、pct_counts_mt这几个列。我习惯先把它们画出来看看分布再决定阈值而不是上来就套一个固定数值。sc.pl.violin(adata, keys[n_genes_by_counts, total_counts, pct_counts_mt], multi_panelTrue)关于阈值不同组织、不同建库方式差异很大但一般可以参考细胞检测到的基因数太少比如小于 200的极可能是空液滴或者破损细胞线粒体基因比例过高比如超过 20%的细胞通常胞浆 RNA 已经大量流失只剩线粒体还在表达。实际过滤的时候多试几个阈值结合小提琴图的分布来定不要一刀切用网上的默认值。sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata adata[adata.obs.pct_counts_mt 20, :] adata adata[adata.obs.n_genes_by_counts 6000, :]高表达基因数超过某个上限的细胞一般是双细胞或者测序异常导致的也需要过滤。这个上界同样要结合数据分布来定。双细胞预测是目前很多项目里容易被跳过的一步。如果你用的是 10x 数据双细胞的比例通常在 0.8% 到 8% 之间样本质量差的时候更高。双细胞的存在会直接影响下游聚类因为它会形成一些看起来像真实细胞类型、实际上却是两个细胞混合的簇。我一般用 scrublet 做这个事import scrublet as scr scrub scr.Scrublet(adata.X, expected_doublet_rate0.06) doublet_scores, predicted_doublets scrub.scrub_doublets(min_counts2, min_cells3, min_gene_variability_pctl85) adata.obs[doublet_score] doublet_scores adata.obs[predicted_doublet] predicted_doublets adata adata[~adata.obs[predicted_doublet], :].copy()注意expected_doublet_rate这个参数不是拍脑袋设的它要根据你上机时的细胞加载量来估算。10x 官方文档里有一个经验表格比如目标回收 5000 个细胞时双细胞率大约是 0.4%回收 10000 个时大约 0.8%回收 20000 个时可以达到 1.6%。如果做的是超高加载量的话这个比例还会更高。3.3 标准化、对数化与高变基因筛选质控完成之后表达矩阵需要先标准化再进入后续分析。scanpy 默认的做法是normalize_totallog1p也就是先让每个细胞的总表达量一致再取对数。sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata)这里的target_sum1e4含义是把每个细胞的所有基因表达量总和归一化到 10000这样不同细胞之间的测序深度差异就被去掉了。取对数则是为了压缩数据动态范围避免少数高表达基因主导后续的降维和距离计算。但这里有一个很容易被忽视的点normalize_total默认的 target_sum 是 1e4这个值实际是 Seurat 里NormalizeData的 scale.factor 默认值的延续。理论上你设置成 1e4、1e5 甚至 1e6对最终结果影响不大因为后续 PCA 之前还会有一步 scale真正的生物学信号不是靠这个绝对值体现的。接下来是高变基因筛选sc.pp.highly_variable_genes(adata, n_top_genes2000, flavorseurat)这里我特意用了flavorseurat它会用 Seurat 里的 vst 方法挑选高变基因这样跟 R 流程出来的结果可比性更强。默认的flavorseurat_v3需要原始 counts 数据如果你已经做了对数化就会报错。高变基因的意义在于基因组里大约有两万个蛋白编码基因但不是每个基因都携带细胞类型分群的信号。大部分基因在不同细胞里表达量都差不多属于背景噪音。选 2000 个左右高变基因相当于把最有区分度的特征挑出来做后续降维计算量小很多聚类效果反而更好。一定要记得选完高变基因之后后续的 PCA、邻居图、聚类都要基于高变基因来做不是直接在全部基因上跑。这一步可以在 AnnData 里做好标记adata.raw adata adata adata[:, adata.var.highly_variable].copy()保存adata.raw很重要因为后面找 marker 基因的时候需要在所有基因的背景下计算而不是局限在高变基因里。3.4 PCA 降维与批次效应处理高变基因筛好之后数据维度还是很高直接拿来做距离计算依然有维数灾难的问题。PCA 把高维数据投影到低维空间同时尽量保留数据中的方差信息。sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack, n_comps50)注意 PCA 之前需要先 scale 数据就是让每个基因的表达量在不同细胞间变成均值为 0、标准差为 1 的分布。这一步必须做否则高表达基因会主导 PCA 结果。max_value10是为了截断极端值防止少数异常高的表达量拉偏主成分。关于主成分数量的选择我现在很少通过碎石图来人工判断了因为样本多的时候这么做很累。我的习惯是直接设一个相对充足的数量比如 30 或 50然后在邻居图构建的时候通过n_neighbors和n_pcs参数来调节。scanpy 默认的流程是 50 个主成分全部拿来构建邻居图但实际上很多成分包含的是无关的噪音。批次效应是单细胞数据分析里绕不开的问题。如果你一个项目里有多个样本或者同一批样本分了好几次上机那不同批次之间的技术差异会混入生物学差异有可能导致聚类时不同批次的细胞各自抱团而不是真正的细胞类型决定聚类。最简单的检查方式是画图的时候按批次上色sc.pl.pca(adata, color[sample_id, cell_type], components[1,2])如果sample_id在 PCA 图上显示出明显的分离模式那就说明存在批次效应需要用工具来校正。常见的选择是 harmony 或 scVI。harmony 的原理是在 PCA 空间里做迭代的软聚类和校正速度快容易上手scVI 是基于深度生成模型的效果通常会更好但对 GPU 和调参要求高一些。harmony 在 scanpy 里的用法也不复杂import scanpy.external as sce sce.pp.harmony_integrate(adata, sample_id)这个函数会把校正后的数据写到adata.obsm[X_pca_harmony]里后面构建邻居图的时候指定用这个嵌入就行。3.5 聚类与 UMAP 可视化PCA 做完之后需要构建细胞之间的邻居图然后基于图结构做聚类和可视化嵌入。这一步在 scanpy 里是一条线走下来的sc.pp.neighbors(adata, n_neighbors15, n_pcs30, use_repX_pca) sc.tl.umap(adata, min_dist0.5) sc.tl.leiden(adata, resolution0.5, flavorigraph, n_iterations2, directedFalse)n_neighbors决定每个细胞跟多少个其它细胞建立连接关系。这个值默认是 15当数据量很大的时候可以适当调低比如 10当样本量少的时候可以调高比如 20 或 30让图结构更稳健。min_dist控制 UMAP 图里点的紧凑程度。值越小不同类别的边界越清晰但簇内会显得很松散值越大簇内越紧凑但不同簇之间可能糊在一起。实际做可视化的时候我通常先用 0.5 看全局再降到 0.1 看局部结构。resolution是聚类粒度的核心参数。它越大得到的簇越多。0.5 是普遍的经验起点。实际使用中我通常会跑一串 resolution 值比如 0.1 到 2然后结合细胞类型 marker 基因表达情况来选最合理的分辨率。没有绝对正确的分辨率只有跟你的生物学问题匹配不匹配。聚类完成后直接画 UMAPsc.pl.umap(adata, color[leiden, sample_id])看到的分群结果如果跟样本批次高度重叠优先回头检查批次校正这一步有没有做到位。3.6 marker 基因鉴定与可视化聚类分完之后必须要回答一个生物学问题每个簇是什么细胞类型这一步靠 marker 基因来完成。scanpy 里用rank_genes_groups来做差异表达分析sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon, n_genes50)这里的method我默认选 wilcoxon 秩和检验它对不同簇的细胞数量差异不敏感比 t 检验更稳健。实际项目里如果每个簇的细胞数很不均衡用 wilcoxon 是更安全的选择。跑完之后结果存在adata.uns[rank_genes_groups]里可以用便捷函数查看每个簇的 top markersc.pl.rank_genes_groups(adata, n_genes20, shareyFalse)也可以用 dotplot 看关键 marker 基因在各簇的表达情况sc.pl.dotplot(adata, var_names[CD3D, CD79A, LYZ, MZB1], groupbyleiden)这一步是经验活。不要机械地只看 top 基因列表要去对照文献和已知的细胞类型 marker。比如 T 细胞的 CD3D、B 细胞的 CD79A、单核细胞的 LYZ、浆细胞的 MZB1这些经典 marker 在每个簇里如果有明确的倾向性表达簇的注释基本就不会错。还有一个小技巧如果某个簇的 marker 全是线粒体基因或者核糖体蛋白基因那大概率这个簇是低质量细胞残留下来的而不是真实细胞类型。遇到这种情况回到质控那一步看看过滤阈值是不是太松了。4. Seurat 与 scanpy 的流程对照4.1 常用函数映射速查从 Seurat 转过来的朋友最痛苦的就是找函数。我整理了一个常用函数的对照表贴在我显示器边上用了很久Seurat (R)scanpy (Python)功能说明CreateSeuratObject / Read10Xsc.read_10x_mtx / sc.read_10x_h5读取数据NormalizeDatasc.pp.normalize_total sc.pp.log1p标准化 对数化FindVariableFeaturessc.pp.highly_variable_genes高变基因筛选ScaleDatasc.pp.scale数据缩放RunPCAsc.tl.pcaPCA 降维RunUMAPsc.tl.umapUMAP 降维/可视化FindNeighborssc.pp.neighbors构建邻居图FindClusterssc.tl.leiden / sc.tl.louvain聚类分群FindAllMarkerssc.tl.rank_genes_groups各组 marker 基因鉴定FeaturePlot / DimPlotsc.pl.umap / sc.pl.embedding可视化VlnPlotsc.pl.violin小提琴图DotPlotsc.pl.dotplot点图展示 marker 表达subsetadata[adata.obs[leiden].isin(...)]选取指定细胞亚群这里特别提醒一个容易搞混的点Seurat 的ScaleData和NormalizeData是两步分开的scanpy 虽然也是分开的函数但很多教程喜欢连着写。如果忘了在 PCA 前做sc.pp.scalePCA 结果会完全变味甚至聚类结果不可信。4.2 结果一致性与差异说明很多人关心的一个问题是同一份数据Seurat 和 scanpy 跑出来的结果会不会完全一致我的经验是核心的细胞分群大体一致但细节上有差异。差异的来源主要有几个方面。第一是高变基因的筛选逻辑不同Seurat 的 vst 方法和 scanpy 的 seurat 风味在算法实现细节上有差异选出来的基因集合略有不同。第二是 PCA 的随机初始化问题如果设了random.seed结果可复现但不设的话每次跑出来的 UMAP 图可能会旋转、镜像。Umap 的嵌入本身带有随机性这是正常现象不影响簇的划分。第三是 leiden 聚类的随机性igrah 和 leidenalg 底层实现不同可能导致部分边界细胞被分到不同簇里。所以我的建议是不要在两个工具必须出一样结果上钻牛角尖。关键验证方式是看细胞类型注释是否一致看关键 marker 基因的富集模式是否一致。如果这两点一致分析结论就是稳健的。项目里要做正式分析的时候还是固定用一套流程不要中途混合切工具这样能保证结果的可复现性。5. 踩坑记录与性能优化实录5.1 版本不兼容与常见报错scanpy 的版本更新比较快即使你照着文档写代码也可能因为版本差异遇到莫名其妙的报错。我遇到过最典型的报错是ValueError: could not convert string to float: Gene Expression这个一般发生在sc.read_10x_mtx读取 feature 文件的时候。原因是新版 Cell Ranger 输出的 features.tsv.gz 文件第一列包含了 gene ID第二列是 gene symbol如果你是直接拿第一列做基因名或者文件里混入了非基因表达的 feature比如抗体捕获的蛋白标签就会触发。解决方法是读取的时候指定var_namesgene_symbols然后在读取后过滤掉不需要的 feature 类型adata adata[:, adata.var[feature_types] Gene Expression].copy()另一个高频报错是TypeError: NoneType object is not subscriptable这个多半是在调用sc.pl.umap或者sc.pl.rank_genes_groups时AnnData 对象里没有对应的 UMAP 坐标或差异分析结果。一定记得按顺序跑流程先sc.tl.umap再sc.pl.umap先sc.tl.rank_genes_groups再sc.pl.rank_genes_groups。还有一个新手特别容易踩的坑用adata.obs[leiden].cat.reorder_categories()重排因子顺序之后直接保存 h5ad再读进来的时候报错说 category 顺序不对。解决方法是保存前先把类型转成字符串adata.obs[leiden] adata.obs[leiden].astype(str)5.2 几十万细胞数据的性能优化心得当你的数据量到了二三十万细胞以上scanpy 默认的设置跑起来就开始吃力了。我有几个实际用过有效优化策略分享出来。第一个是控制高变基因的数量。默认n_top_genes2000在十万细胞级别速度还能接受到二三十万细胞的时候每一步都变慢。可以适当减少到 1500 甚至 1000聚类结果差别不大但速度能快不少。第二个是用sc.pp.neighbors的methodpynndescent配合n_neighbors10在高维空间中构建邻居图时比默认方法快得多实测在二十万细胞的数据上能快三到五倍。第三个是尽量使用稀疏矩阵。AnnData 默认从 10x 读取的矩阵就是稀疏的但有些操作如果你不小心转成了稠密矩阵内存会瞬间爆炸。比如有人习惯用.toarray()去看数据看完之后那个变量还留在内存里后面所有操作就都慢下来了。建议养成好习惯不要随便 materialize 稀疏矩阵。如果条件允许可以考虑把数据迁移到 GPU 上跑用 rapids-singlecell 这个库它实现了 scanpy 大部分核心函数的 GPU 版本。我自己在单张 RTX 4090 上跑过二十万细胞的完整流程从 PCA 到 UMAP 到 Leiden 聚类耗时不到原来 CPU 版本的十分之一。但注意 GPU 方案对显存的要求不低建议至少 12G 显存以上再尝试。5.3 数据格式兼容与保存的细节多步骤分析做完一定要及时保存中间结果。scanpy 保存的标准格式是 h5adadata.write(analysis_result.h5ad)读取的时候用sc.read_h5ad(analysis_result.h5ad)就能完整恢复。但有几个细节要提醒adata.raw默认不会被 saven 到 h5ad 之外单独保存但它本身会跟 h5ad 一起写入所以不用额外操心。如果数据量很大可以设置压缩参数adata.write(result.h5ad, compressiongzip)文件体积能缩小不少代价是后续读取稍微慢一点。跟 Seurat 互转的时候推荐使用sce库里的sce.pp.read_seurat或者手动导出稀疏矩阵再读。要是直接拿 R 里 saveRDS 的文件给 Python 用是行不通的。如果你需要配合 Seurat 的流程对比结果可以先用 scanpy 导出关键结果表格adata.obs[[leiden, cell_type]].to_csv(cluster_annotation.csv) result sc.get.rank_genes_groups_df(adata, groupNone) result.to_csv(marker_genes.csv, indexFalse)这样 R 和 Python 两端都能用同一份注释信息方便后续交叉验证。6. 最后分享两个实用经验第一个是关于分析的记录习惯。我现在会在每个分析项目里写一份环境说明包含 scanpy、numpy、pandas、python-igraph 这几个核心包的版本号。不要觉得这是形式主义过两个月你自己回来看代码没有版本记录的话可能连报错都复现不出来。用 conda env export 导出环境配置是最省事的。第二个是关于学习路径的建议。如果你刚开始接触 scanpy不要急着把整个流程每一个参数都背下来。先用一套公开数据集比如 PBMC 3k把标准流程完整跑一遍跑通之后再回过来逐个理解参数的意义。先跑通、再调优这条路径比我一开始逐行研究文档要有效得多也是我在这个领域踩了一圈坑之后最想告诉你的经验。
做单细胞数据分析的人多半绕不开 Seurat。我前面两篇教程把 R 语言里 Seurat 的完整流程写过一遍从数据读入、质控、标准化到聚类和 marker 基因鉴定一套跑下来确实顺手。但这几年我自己的项目里出现了一个很现实的情况样本量越来越大、数据动辄几十万细胞R 的内存开销开始吃紧同时分析流程里需要跟 Python 生态里的深度学习模型、可视化库、scanpy 之外的其它工具做配合用 R 做完整流程就感觉两头受堵。于是我开始把整套单细胞转录组测序数据分析流程往 python 的 scanpy 上迁。这篇教程就是给那些已经在用或者准备用 scanpy 的朋友准备的尤其是从 Seurat 转过来的那一批人。我不打算把官方文档抄一遍而是站在一个实际跑过大量样本的角度把 scanpy 的完整分析流程拆开来讲环境怎么搭、AnnData 对象怎么理解、每一步参数为什么要这么设、跟 Seurat 的哪些函数能对上号以及我踩过的那些坑。内容覆盖从 10x 原始数据到 marker 基因输出的全流程适合已经会一点 Python、但还没系统跑过单细胞分析的人参考。1. 从 Seurat 转向 scanpy 的真实理由与生态对比1.1 为什么你该考虑在 Python 里做单细胞分析先说结论我不是劝你放弃 Seurat。如果你是纯粹做单细胞分析、不做太多定制化开发Seurat 依然是很成熟的选择文档清楚、社区活跃、画图也漂亮。但当你的任务开始变复杂的时候Python 这套生态的优势就出来了。第一个理由是规模。Seurat 的对象设计本身是没有问题的但它默认把大部分数据加载进内存表达矩阵加上降维结果、邻接图这些一个十万细胞的数据集跑完 UMAP 和聚类内存随便就吃掉二三十个 G。scanpy 底层用 AnnData 对象基于 h5py 和 scipy 稀疏矩阵能够更有效地存储和处理大规模数据配合 out-of-core 读取方式处理百万级细胞也不是稀罕事。第二个理由是生态联动。单细胞数据分析从来不只停在聚类和 marker现在大家跑到最后都要接拟时序、细胞通讯、转录因子分析。拟时序常用的 Monocle 在 R 里但 Python 里也有 scVeloRNA velocity、CellRank 这些更往动力学方向走的工具。你要做深度学习批处理、做细胞类型注释Python 这边直接可以调用不需要跨语言搬运数据。第三个理由很实际我的合作方越来越多人交出来的数据预处理结果都是 h5ad 格式跟 Python 的 scanpy 无缝衔接。如果在 R 里跑完 Seurat想转出去给别人做下游还要过一层格式转换麻烦不说还有可能丢失一部分元数据信息。1.2 scanpy 与 Seurat 的定位差异scanpy 的设计哲学和 Seurat 不太一样。Seurat 是把一整套流程封装得比较完整你按着标准流程走UI 式的体验默认参数基本能给出一个可接受的结果。scanpy 更像是一个工具箱它提供的是模块化的函数每一步你都可以单独调用也可以自由组合。这种差异在调试的时候体现得很明显。用 Seurat 的时候如果一段流程报错你往往需要去理解它内部完整的 S4 对象结构看它是怎么存储和传递数据的。而 scanpy 的函数大多是接收 AnnData、返回修改后的 AnnData你对每一步之间数据是怎么流动的心里是很清楚的。还有一个差异在可视化上。Seurat 的画图函数内置度很高各种 DimPlot、FeaturePlot、VlnPlot 都给你配好了。scanpy 的绘图功能其实也不弱而且它是基于 matplotlib 的所以你对图表的控制力更强——你想调整坐标轴、字体、颜色、布局直接在 matplotlib 层面上做就好不用去破解 Seurat 的画图函数参数。但也因为这种灵活性你有时候需要自己多写几行代码来调样式不能像 Seurat 那样开箱即用。2. 环境搭建与安装避坑2.1 用 conda 建环境比直接用 pip 省心十倍很多刚开始玩 python 单细胞分析的朋友第一步就栽在环境上。直接在系统 Python 里 pip install scanpy多半会把依赖搞乱尤其是 numpy、scipy、pandas 这些科学计算库的版本稍有不慎就会冲突。我自己用的方式是 conda 建独立环境一步到位干净隔离。创建环境的命令很常规但有几个小地方要留意conda create -n scanpy python3.9 -y conda activate scanpy我推荐用 Python 3.9 而不是最新版本。不是说我保守而是 scanpy 生态里很多依赖库在最新 Python 版本下可能会出现预编译包缺失的问题conda 要现场编译的话会很痛苦。3.9 是当前兼容性最好的版本大部分生物信息学工具都能在上面稳定运行。接着装核心包conda install -c conda-forge scanpy python-igraph leidenalg -y这里有一个很关键的细节leidenalg 和 python-igraph 一定要一起装。leiden 聚类算法是 scanpy 做聚类时的默认选择比早期的 louvain 算法在社区划分质量上更好。但 leidenalg 的安装经常出问题尤其是当你只通过 pip 安装的时候很可能会遇到编译错误。用 conda-forge 渠道安装它会自动帮你匹配好 igraph 和 leidenalg 的版本省去很多麻烦。安装完之后建议顺手装几个下游常用包pip install scrublet pip install scvi-toolsscrublet 是用来预测双细胞的scvi-tools 是深度学习的批次整合工具。当然这些不是必须要装的看你项目需求来定。装完之后用python -c import scanpy as sc; print(sc.__version__)验证一下能出来版本号就说明环境没问题。2.2 关于国内镜像源的取舍这一步纯粹是实际体验问题。conda 和 pip 默认的源在下载速度上确实不稳定国内的朋友经常会遇到下载到一半超时的情况。我可以明确地跟你说配置国内镜像源是合规且非常常见的做法能显著提升安装效率。比如 pip 源我长期用清华的 PyPI 镜像pip config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simpleconda 源也是一样在~/.condarc里配置一下清华的 conda-forge 和 defaults 镜像。但要注意不要同时在.condarc里混合太多不同来源的 channel不然 conda 在解析依赖时会非常慢甚至出现 channel 优先级冲突导致装不上包的问题。还有一个小建议装包的时候尽量用mamba替代condamamba 用的是 C 写的求解器依赖解析速度比 conda 快很多尤其是环境里要装的包比较多的时候差距是数量级的。3. scanpy 完整分析流程实操3.1 数据读入与 AnnData 对象初识拿到单细胞转录组测序数据之后第一步肯定是把数据读进来。scanpy 支持很多格式但我实际项目里碰到最多的还是 10x Genomics 的 Cell Ranger 输出。如果是 Cell Ranger 的过滤矩阵输出目录通常包含 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz 三个文件读取方式如下import scanpy as sc import numpy as np adata sc.read_10x_mtx( data/filtered_feature_bc_matrix/, var_namesgene_symbols, make_uniqueTrue )如果是单个的 h5 文件就用sc.read_10x_h5()。读取之后adata 就是一个 AnnData 对象它的核心结构是adata.X存的是细胞×基因的表达矩阵。这个矩阵在 Python 里是稀疏矩阵格式这也是为什么 scanpy 能处理大规模数据的关键之一——大部分基因在大部分细胞里表达量是 0稀疏矩阵只存非零值内存占用可以压缩到原来的几十分之一。拿到数据之后第一步不要急着跑分析先看一眼数据的基本面貌print(adata.shape) print(adata.var.head()) print(adata.obs.head())adata.obs存的是细胞层面的元数据adata.var存的是基因层面的元数据。后面每一步分析的结果比如质控指标、聚类标签、UMAP 坐标都会存到这两个 DataFrame 里。理解 AnnData 这个结构后面的逻辑就顺了。3.2 质控过滤线粒体基因与双细胞质控是整个流程里最影响结果的一步这一枪打得不准后面再怎么调参都白搭。我一般会做三件事过滤低质量细胞、过滤低表达基因、预测并去除双细胞。先计算质控指标adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, log1pFalse, inplaceTrue)这一步会在adata.obs里生成n_genes_by_counts、total_counts、pct_counts_mt这几个列。我习惯先把它们画出来看看分布再决定阈值而不是上来就套一个固定数值。sc.pl.violin(adata, keys[n_genes_by_counts, total_counts, pct_counts_mt], multi_panelTrue)关于阈值不同组织、不同建库方式差异很大但一般可以参考细胞检测到的基因数太少比如小于 200的极可能是空液滴或者破损细胞线粒体基因比例过高比如超过 20%的细胞通常胞浆 RNA 已经大量流失只剩线粒体还在表达。实际过滤的时候多试几个阈值结合小提琴图的分布来定不要一刀切用网上的默认值。sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) adata adata[adata.obs.pct_counts_mt 20, :] adata adata[adata.obs.n_genes_by_counts 6000, :]高表达基因数超过某个上限的细胞一般是双细胞或者测序异常导致的也需要过滤。这个上界同样要结合数据分布来定。双细胞预测是目前很多项目里容易被跳过的一步。如果你用的是 10x 数据双细胞的比例通常在 0.8% 到 8% 之间样本质量差的时候更高。双细胞的存在会直接影响下游聚类因为它会形成一些看起来像真实细胞类型、实际上却是两个细胞混合的簇。我一般用 scrublet 做这个事import scrublet as scr scrub scr.Scrublet(adata.X, expected_doublet_rate0.06) doublet_scores, predicted_doublets scrub.scrub_doublets(min_counts2, min_cells3, min_gene_variability_pctl85) adata.obs[doublet_score] doublet_scores adata.obs[predicted_doublet] predicted_doublets adata adata[~adata.obs[predicted_doublet], :].copy()注意expected_doublet_rate这个参数不是拍脑袋设的它要根据你上机时的细胞加载量来估算。10x 官方文档里有一个经验表格比如目标回收 5000 个细胞时双细胞率大约是 0.4%回收 10000 个时大约 0.8%回收 20000 个时可以达到 1.6%。如果做的是超高加载量的话这个比例还会更高。3.3 标准化、对数化与高变基因筛选质控完成之后表达矩阵需要先标准化再进入后续分析。scanpy 默认的做法是normalize_totallog1p也就是先让每个细胞的总表达量一致再取对数。sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata)这里的target_sum1e4含义是把每个细胞的所有基因表达量总和归一化到 10000这样不同细胞之间的测序深度差异就被去掉了。取对数则是为了压缩数据动态范围避免少数高表达基因主导后续的降维和距离计算。但这里有一个很容易被忽视的点normalize_total默认的 target_sum 是 1e4这个值实际是 Seurat 里NormalizeData的 scale.factor 默认值的延续。理论上你设置成 1e4、1e5 甚至 1e6对最终结果影响不大因为后续 PCA 之前还会有一步 scale真正的生物学信号不是靠这个绝对值体现的。接下来是高变基因筛选sc.pp.highly_variable_genes(adata, n_top_genes2000, flavorseurat)这里我特意用了flavorseurat它会用 Seurat 里的 vst 方法挑选高变基因这样跟 R 流程出来的结果可比性更强。默认的flavorseurat_v3需要原始 counts 数据如果你已经做了对数化就会报错。高变基因的意义在于基因组里大约有两万个蛋白编码基因但不是每个基因都携带细胞类型分群的信号。大部分基因在不同细胞里表达量都差不多属于背景噪音。选 2000 个左右高变基因相当于把最有区分度的特征挑出来做后续降维计算量小很多聚类效果反而更好。一定要记得选完高变基因之后后续的 PCA、邻居图、聚类都要基于高变基因来做不是直接在全部基因上跑。这一步可以在 AnnData 里做好标记adata.raw adata adata adata[:, adata.var.highly_variable].copy()保存adata.raw很重要因为后面找 marker 基因的时候需要在所有基因的背景下计算而不是局限在高变基因里。3.4 PCA 降维与批次效应处理高变基因筛好之后数据维度还是很高直接拿来做距离计算依然有维数灾难的问题。PCA 把高维数据投影到低维空间同时尽量保留数据中的方差信息。sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack, n_comps50)注意 PCA 之前需要先 scale 数据就是让每个基因的表达量在不同细胞间变成均值为 0、标准差为 1 的分布。这一步必须做否则高表达基因会主导 PCA 结果。max_value10是为了截断极端值防止少数异常高的表达量拉偏主成分。关于主成分数量的选择我现在很少通过碎石图来人工判断了因为样本多的时候这么做很累。我的习惯是直接设一个相对充足的数量比如 30 或 50然后在邻居图构建的时候通过n_neighbors和n_pcs参数来调节。scanpy 默认的流程是 50 个主成分全部拿来构建邻居图但实际上很多成分包含的是无关的噪音。批次效应是单细胞数据分析里绕不开的问题。如果你一个项目里有多个样本或者同一批样本分了好几次上机那不同批次之间的技术差异会混入生物学差异有可能导致聚类时不同批次的细胞各自抱团而不是真正的细胞类型决定聚类。最简单的检查方式是画图的时候按批次上色sc.pl.pca(adata, color[sample_id, cell_type], components[1,2])如果sample_id在 PCA 图上显示出明显的分离模式那就说明存在批次效应需要用工具来校正。常见的选择是 harmony 或 scVI。harmony 的原理是在 PCA 空间里做迭代的软聚类和校正速度快容易上手scVI 是基于深度生成模型的效果通常会更好但对 GPU 和调参要求高一些。harmony 在 scanpy 里的用法也不复杂import scanpy.external as sce sce.pp.harmony_integrate(adata, sample_id)这个函数会把校正后的数据写到adata.obsm[X_pca_harmony]里后面构建邻居图的时候指定用这个嵌入就行。3.5 聚类与 UMAP 可视化PCA 做完之后需要构建细胞之间的邻居图然后基于图结构做聚类和可视化嵌入。这一步在 scanpy 里是一条线走下来的sc.pp.neighbors(adata, n_neighbors15, n_pcs30, use_repX_pca) sc.tl.umap(adata, min_dist0.5) sc.tl.leiden(adata, resolution0.5, flavorigraph, n_iterations2, directedFalse)n_neighbors决定每个细胞跟多少个其它细胞建立连接关系。这个值默认是 15当数据量很大的时候可以适当调低比如 10当样本量少的时候可以调高比如 20 或 30让图结构更稳健。min_dist控制 UMAP 图里点的紧凑程度。值越小不同类别的边界越清晰但簇内会显得很松散值越大簇内越紧凑但不同簇之间可能糊在一起。实际做可视化的时候我通常先用 0.5 看全局再降到 0.1 看局部结构。resolution是聚类粒度的核心参数。它越大得到的簇越多。0.5 是普遍的经验起点。实际使用中我通常会跑一串 resolution 值比如 0.1 到 2然后结合细胞类型 marker 基因表达情况来选最合理的分辨率。没有绝对正确的分辨率只有跟你的生物学问题匹配不匹配。聚类完成后直接画 UMAPsc.pl.umap(adata, color[leiden, sample_id])看到的分群结果如果跟样本批次高度重叠优先回头检查批次校正这一步有没有做到位。3.6 marker 基因鉴定与可视化聚类分完之后必须要回答一个生物学问题每个簇是什么细胞类型这一步靠 marker 基因来完成。scanpy 里用rank_genes_groups来做差异表达分析sc.tl.rank_genes_groups(adata, leiden, methodwilcoxon, n_genes50)这里的method我默认选 wilcoxon 秩和检验它对不同簇的细胞数量差异不敏感比 t 检验更稳健。实际项目里如果每个簇的细胞数很不均衡用 wilcoxon 是更安全的选择。跑完之后结果存在adata.uns[rank_genes_groups]里可以用便捷函数查看每个簇的 top markersc.pl.rank_genes_groups(adata, n_genes20, shareyFalse)也可以用 dotplot 看关键 marker 基因在各簇的表达情况sc.pl.dotplot(adata, var_names[CD3D, CD79A, LYZ, MZB1], groupbyleiden)这一步是经验活。不要机械地只看 top 基因列表要去对照文献和已知的细胞类型 marker。比如 T 细胞的 CD3D、B 细胞的 CD79A、单核细胞的 LYZ、浆细胞的 MZB1这些经典 marker 在每个簇里如果有明确的倾向性表达簇的注释基本就不会错。还有一个小技巧如果某个簇的 marker 全是线粒体基因或者核糖体蛋白基因那大概率这个簇是低质量细胞残留下来的而不是真实细胞类型。遇到这种情况回到质控那一步看看过滤阈值是不是太松了。4. Seurat 与 scanpy 的流程对照4.1 常用函数映射速查从 Seurat 转过来的朋友最痛苦的就是找函数。我整理了一个常用函数的对照表贴在我显示器边上用了很久Seurat (R)scanpy (Python)功能说明CreateSeuratObject / Read10Xsc.read_10x_mtx / sc.read_10x_h5读取数据NormalizeDatasc.pp.normalize_total sc.pp.log1p标准化 对数化FindVariableFeaturessc.pp.highly_variable_genes高变基因筛选ScaleDatasc.pp.scale数据缩放RunPCAsc.tl.pcaPCA 降维RunUMAPsc.tl.umapUMAP 降维/可视化FindNeighborssc.pp.neighbors构建邻居图FindClusterssc.tl.leiden / sc.tl.louvain聚类分群FindAllMarkerssc.tl.rank_genes_groups各组 marker 基因鉴定FeaturePlot / DimPlotsc.pl.umap / sc.pl.embedding可视化VlnPlotsc.pl.violin小提琴图DotPlotsc.pl.dotplot点图展示 marker 表达subsetadata[adata.obs[leiden].isin(...)]选取指定细胞亚群这里特别提醒一个容易搞混的点Seurat 的ScaleData和NormalizeData是两步分开的scanpy 虽然也是分开的函数但很多教程喜欢连着写。如果忘了在 PCA 前做sc.pp.scalePCA 结果会完全变味甚至聚类结果不可信。4.2 结果一致性与差异说明很多人关心的一个问题是同一份数据Seurat 和 scanpy 跑出来的结果会不会完全一致我的经验是核心的细胞分群大体一致但细节上有差异。差异的来源主要有几个方面。第一是高变基因的筛选逻辑不同Seurat 的 vst 方法和 scanpy 的 seurat 风味在算法实现细节上有差异选出来的基因集合略有不同。第二是 PCA 的随机初始化问题如果设了random.seed结果可复现但不设的话每次跑出来的 UMAP 图可能会旋转、镜像。Umap 的嵌入本身带有随机性这是正常现象不影响簇的划分。第三是 leiden 聚类的随机性igrah 和 leidenalg 底层实现不同可能导致部分边界细胞被分到不同簇里。所以我的建议是不要在两个工具必须出一样结果上钻牛角尖。关键验证方式是看细胞类型注释是否一致看关键 marker 基因的富集模式是否一致。如果这两点一致分析结论就是稳健的。项目里要做正式分析的时候还是固定用一套流程不要中途混合切工具这样能保证结果的可复现性。5. 踩坑记录与性能优化实录5.1 版本不兼容与常见报错scanpy 的版本更新比较快即使你照着文档写代码也可能因为版本差异遇到莫名其妙的报错。我遇到过最典型的报错是ValueError: could not convert string to float: Gene Expression这个一般发生在sc.read_10x_mtx读取 feature 文件的时候。原因是新版 Cell Ranger 输出的 features.tsv.gz 文件第一列包含了 gene ID第二列是 gene symbol如果你是直接拿第一列做基因名或者文件里混入了非基因表达的 feature比如抗体捕获的蛋白标签就会触发。解决方法是读取的时候指定var_namesgene_symbols然后在读取后过滤掉不需要的 feature 类型adata adata[:, adata.var[feature_types] Gene Expression].copy()另一个高频报错是TypeError: NoneType object is not subscriptable这个多半是在调用sc.pl.umap或者sc.pl.rank_genes_groups时AnnData 对象里没有对应的 UMAP 坐标或差异分析结果。一定记得按顺序跑流程先sc.tl.umap再sc.pl.umap先sc.tl.rank_genes_groups再sc.pl.rank_genes_groups。还有一个新手特别容易踩的坑用adata.obs[leiden].cat.reorder_categories()重排因子顺序之后直接保存 h5ad再读进来的时候报错说 category 顺序不对。解决方法是保存前先把类型转成字符串adata.obs[leiden] adata.obs[leiden].astype(str)5.2 几十万细胞数据的性能优化心得当你的数据量到了二三十万细胞以上scanpy 默认的设置跑起来就开始吃力了。我有几个实际用过有效优化策略分享出来。第一个是控制高变基因的数量。默认n_top_genes2000在十万细胞级别速度还能接受到二三十万细胞的时候每一步都变慢。可以适当减少到 1500 甚至 1000聚类结果差别不大但速度能快不少。第二个是用sc.pp.neighbors的methodpynndescent配合n_neighbors10在高维空间中构建邻居图时比默认方法快得多实测在二十万细胞的数据上能快三到五倍。第三个是尽量使用稀疏矩阵。AnnData 默认从 10x 读取的矩阵就是稀疏的但有些操作如果你不小心转成了稠密矩阵内存会瞬间爆炸。比如有人习惯用.toarray()去看数据看完之后那个变量还留在内存里后面所有操作就都慢下来了。建议养成好习惯不要随便 materialize 稀疏矩阵。如果条件允许可以考虑把数据迁移到 GPU 上跑用 rapids-singlecell 这个库它实现了 scanpy 大部分核心函数的 GPU 版本。我自己在单张 RTX 4090 上跑过二十万细胞的完整流程从 PCA 到 UMAP 到 Leiden 聚类耗时不到原来 CPU 版本的十分之一。但注意 GPU 方案对显存的要求不低建议至少 12G 显存以上再尝试。5.3 数据格式兼容与保存的细节多步骤分析做完一定要及时保存中间结果。scanpy 保存的标准格式是 h5adadata.write(analysis_result.h5ad)读取的时候用sc.read_h5ad(analysis_result.h5ad)就能完整恢复。但有几个细节要提醒adata.raw默认不会被 saven 到 h5ad 之外单独保存但它本身会跟 h5ad 一起写入所以不用额外操心。如果数据量很大可以设置压缩参数adata.write(result.h5ad, compressiongzip)文件体积能缩小不少代价是后续读取稍微慢一点。跟 Seurat 互转的时候推荐使用sce库里的sce.pp.read_seurat或者手动导出稀疏矩阵再读。要是直接拿 R 里 saveRDS 的文件给 Python 用是行不通的。如果你需要配合 Seurat 的流程对比结果可以先用 scanpy 导出关键结果表格adata.obs[[leiden, cell_type]].to_csv(cluster_annotation.csv) result sc.get.rank_genes_groups_df(adata, groupNone) result.to_csv(marker_genes.csv, indexFalse)这样 R 和 Python 两端都能用同一份注释信息方便后续交叉验证。6. 最后分享两个实用经验第一个是关于分析的记录习惯。我现在会在每个分析项目里写一份环境说明包含 scanpy、numpy、pandas、python-igraph 这几个核心包的版本号。不要觉得这是形式主义过两个月你自己回来看代码没有版本记录的话可能连报错都复现不出来。用 conda env export 导出环境配置是最省事的。第二个是关于学习路径的建议。如果你刚开始接触 scanpy不要急着把整个流程每一个参数都背下来。先用一套公开数据集比如 PBMC 3k把标准流程完整跑一遍跑通之后再回过来逐个理解参数的意义。先跑通、再调优这条路径比我一开始逐行研究文档要有效得多也是我在这个领域踩了一圈坑之后最想告诉你的经验。
