讲真最近在生信交流群里看到不少朋友把 Seurat 跑得飞起但一提到 python 生态就头大。这个系列前面已经写了 Seurat 做单细胞转录组的标准流程这一篇我专门来讲讲 python 这边的玩法核心是 scanpy。如果你正在纠结要不要从 R 迁移到 python或者老板突然让你把 Seurat 的流程换成 scanpy 重跑一遍那这篇文章就是给你准备的。我会直接用一套真实项目中跑通的代码流程把 scanpy 从环境搭建、AnnData 数据结构到质控、归一化、聚类、marker 注释的完整链路过一遍并且会重点对比 scanpy 和 Seurat 在处理同一批数据时到底哪里不一样、有哪些坑。这样不管是刚入门的新手还是已经熟悉 Seurat 想换生态的老手都能拿着文章里的代码直接改着用。1. 先想清楚为什么还要用 scanpy不要觉得这是重复造轮子。Seurat 和 scanpy 各自背后是 R 和 python 两套完全不同的生态而这个差异在实际项目中会直接影响你的分析效率、算法选择甚至最后发文章的审稿观感。1.1 Seurat 和 scanpy 的生态定位差异Seurat 是 2015 年左右从 Satija 实验室出来的早期核心优势是把单细胞分析流程封装得极其友好一行NormalizeData、一行FindMarkers对湿实验出身的朋友来说非常“傻瓜”。但 R 的内存管理在大规模数据上确实吃力1 个 10x 的 10 万细胞样本跑下来动辄几十 G 内存稍微叠加多个样本就变得很痛苦。scanpy 是德国那边实验室主导开发的构建在 anndata 和 numpy/scipy 之上底层向量化做得更好处理大规模数据时内存占用明显比 R 系友好而且 scanpy 和 python 的机器学习库比如 sklearn、umap-learn衔接天然通畅。你如果在做 atlas 级别几十万到上百万细胞的项目scanpy 的concat和批量处理能力比 Seurat 顺手很多。1.2 scanpy 的优势场景与适用人群我见过三类人特别适合转向 scanpy第一类是计划做算法开发或深度学习的。现在很多单细胞大模型如 scGPT、Geneformer的输入就是 AnnData 或 h5ad 格式用 scanpy 做预处理是最顺的路线绕开 R 的数据交换问题。第二类是处理超大数据的。scanpy 的sc.pp.pca用的随机 SVD计算速度在百万细胞级别依然可接受配合scanpy.external里的 harmony 集成多批次整合也比 R 那边配置环境更省心。第三类是纯 python 写代码习惯的人。对这些人来说每写一行代码都切回 R 的代价太高scanpy 可以让他们在同一个 notebook 里完成预处理、可视化、机器学习建模。当然Seurat 也一直在进步比如SCTransform、bridge这些方法依然很强。我的观点是不是让你二选一而是让你有能力双持。数据量小、团队全是 R 用户你用 Seurat 没问题但如果你要扩大分析规模、或者准备接 python 生态的下游工具那 scanpy 必须得会。2. 环境准备先把 python scanpy 跑起来这一步看着简单但我实测下来是坑最多的环节。很多人的第一步不是死在代码逻辑上而是死在 conda 环境、依赖冲突、还有 igraph 或 leidenalg 装不上这类问题。2.1 python 环境与安装方式强烈建议不要直接往系统 python 里装 scanpy这跟直接在 R 里装几百个包然后互相冲突是一个道理。你用 conda 或 mamba 建一个独立环境后面想怎么折腾都行。我是用 mamba 的因为 conda 在处理 scanpy 这种依赖树很深的包时慢得让人抓狂。# 建一个干净的环境指定 python 版本 mamba create -n scanpy python3.10 -y conda activate scanpy # 安装 scanpy 和常用周边 mamba install -c conda-forge scanpy python-igraph leidenalg -y pip install anndata pip install harmonypy # harmony 批次整合 pip install scanpy[tools]python 3.10 是目前兼容性和库支持最稳的版本python 3.12 虽然新但一些老代码里的 numba 或 anndata 旧版本可能有兼容问题。生产环境建议老老实实用 3.10。2.2 scanpy 安装与依赖注意scanpy 最关键的依赖是anndata、numpy、scipy、pandas、matplotlib、scikit-learn、umap-learn、leidenalg或python-igraph。其中leidenalg 是聚类绕不开的环节但它依赖python-igraph这两个包之间版本如果不匹配会直接报ModuleNotFoundError或TypeError。我踩过最典型的一个坑用 pip 单独装leidenalg时它自己带了一个igraph但跟你环境里已经存在的python-igraph冲突导致sc.tl.leiden直接报错。解决办法是用 mamba 同时装这两个包让 conda 帮你解析依赖关系mamba install -c conda-forge python-igraph leidenalg -y2.3 加速方案与国内安装体验如果你在国内conda 默认源速度会非常感人。可以换清华源或阿里源配置~/.condarc把default_channels和custom_channels都指到国内镜像。pip 这边同样可以用清华的 PyPI 镜像这样下载 scanpy 和它那一堆依赖包的速度能快一个数量级。我自己的习惯是conda-forge 装重型依赖igraph、leidenalg、scanpy 本身pip 装更新较快的纯 python 包如 harmonypy。这样不容易卡在某个包上。3. AnnDatascanpy 的核心数据结构理解了 AnnDatascanpy 就学会了一半。它跟 Seurat 不一样Seurat 的分析对象是Seurat对象一个 S4 对象里塞了多套 assay而 scanpy 的核心对象叫AnnData设计思路更接近“带注释的数据矩阵”。3.1 AnnData 结构逐层拆解一个典型的 AnnData 对象长这样属性作用对应 Seurat 的位置adata.X主表达矩阵通常是稀疏矩阵行是细胞列是基因GetAssayData(object, slotcounts)adata.obs细胞 metadata 的 DataFrame行名是细胞 barcodeobjectmeta.dataadata.var基因 metadata 的 DataFrame行名是基因名objectassays$RNAmeta.featuresadata.obsm降维结果PCA、UMAP、tSNE 等字典类型objectreductionsadata.varm基因在降维空间的加载矩阵Loadings(object, reductionpca)adata.uns非结构化数据存聚类结果、marker 列表等objectmiscadata.layers可存储多个表达矩阵比如 raw counts 和 normalizedSlot 里不同矩阵这里有一个新手容易懵的地方adata.X到底是 counts 还是 normalized 数据scanpy 不像 Seurat 那样把RNAassay 下的counts、data、scale.data分开存好。adata.X的内容完全由你自己控制所以你的分析流程必须严格注意在哪个步骤修改了adata.X否则后面跑完归一化可能再也拿不到原始 counts 了。我的建议是一开始就把原始 counts 备份在adata.layers[counts]里import scanpy as sc adata sc.read_10x_h5(filtered_feature_bc_matrix.h5) adata.layers[counts] adata.X.copy() # 存原始 counts3.2 从 10x 数据读入 AnnDatascanpy 支持多种读取方式最常用的是sc.read_10x_h5和sc.read_10x_mtx。import scanpy as sc # 从 h5 文件读取 adata sc.read_10x_h5(filtered_feature_bc_matrix.h5) # 从 mtx 目录读取10x 的 cellranger 标准输出 adata sc.read_10x_mtx(filtered_gene_bc_matrices/hg19/, var_namesgene_symbols, make_uniqueTrue)注意read_10x_mtx的make_uniqueTrue参数如果基因名有重复比如某些版本会同时输出基因符号和基因 ID必须设置去重否则后面.var_names重复会导致一系列报错。另外读取后要检查一下adata维度是不是符合预期常见的 cellranger 输出过滤后矩阵一般细胞数在几千到几万基因数在 2 万左右。3.3 和 Seurat 对象互转的实用经验如果团队里 R 和 python 两边都有人互转是绕不开的。scanpy 官方推荐的方式是通过sceasy库转成.rdsimport sceasy # anndata 转 Seurat需要本地装了 R 的 Seurat sceasy.convert(adata, toseurat, outFileadata.rds) # Seurat 转 anndata sceasy.convert(adata, toanndata, outFileadata.h5ad)不过在实际项目中我更喜欢直接操作 h5ad 文件adata.write(data.h5ad)之后R 那边用SeuratDisk的Convert和LoadH5Seurat函数读进来有时候比 sceasy 更稳尤其是对象特别大的时候。互转之后记得检查基因名大小写问题scanpy 默认保留原始大小写R 读入后可能变成首字母大写不统一就会丢基因。4. 标准流程实操从 raw counts 到聚类注释这一节是整篇最核心的部分。我会按真实项目顺序把 scanpy 的完整分析流程走一遍每一步都会对比 Seurat 的对应操作方便你从 R 迁移时能快速对上号。4.1 质量控制QC线粒体、核糖体、双细胞QC 的逻辑跟 Seurat 完全一致只是代码表现不同。核心指标就三个每个细胞的基因数n_genes_by_counts、总 UMI 数total_counts、线粒体基因比例pct_counts_mt。import scanpy as sc import matplotlib.pyplot as plt adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], percent_topNone, inplaceTrue) sc.pl.violin(adata, keys[n_genes_by_counts, total_counts, pct_counts_mt], multi_panelTrue)阈值怎么定Seurat 经典教程里一般用 nFeature_RNA 200 且 2500percent.mt 5%。但实际项目里千万别死板套用如果是组织样本比如肿瘤组织细胞的线粒体比例天然比细胞系高8% 甚至 10% 都可能是正常细胞群。如果是冷冻组织解离后的数据基因数普遍偏低阈值可以适当放宽。我的做法是先画小提琴图看分布再结合总 UMI 数的双峰分布来定。QC 过滤这步用 sc.pp.filter_cells 和 sc.pp.filter_genes 即可# 先粗略过滤pac 掉明显低质量的 sc.pp.filter_cells(adata, min_genes200) sc.pp.filter_genes(adata, min_cells3) # 再根据 QC 指标过滤 adata adata[adata.obs.n_genes_by_counts 6000, :].copy() adata adata[adata.obs.pct_counts_mt 20, :].copy()双细胞过滤方面scanpy 生态里常用scrublet跟 Seurat 的DoubletFinder是同一个思路import scrublet as scr counts_matrix adata.X.T # scrublet 需要 cell × gene 矩阵 scrub scr.Scrublet(counts_matrix, expected_doublet_rate0.06) doublet_scores, predicted_doublets scrub.scrub_doublets() adata.obs[doublet_score] doublet_scores adata.obs[predicted_doublet] predicted_doublets实操中注意不同样本的 doublet rate 会不一样expected_doublet_rate不要统一用 0.06。比如 10x 的 8000 细胞捕获量doublet rate 可能到 5%-8%但如果是低捕获量的样本比如 3000 细胞doublet rate 通常不到 3%。4.2 归一化Normalization与 log1pscanpy 的数据归一化流程和 Seurat 有概念上的差异。Seurat 的NormalizeData默认是 log1p(CPM/100)也就是 log1p(normalize total1e4)。scanpy 里对应的操作是sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata)这两行代码合起来就等于 Seurat 的默认NormalizeData。但scanpy 不会把归一化后的数据存到scale.data它直接把adata.X给覆盖了。所以我前面才强调必须在归一化之前先把原始 counts 存到adata.layers[counts]否则你想用 MAST 或 wilcox 做差异分析时没有原始 counts 会很被动。如果你之后要做sc.pp.scaleZ-score 标准化目的是让每个基因的均值 0、方差 1注意 scanpy 的scale默认是zero_centerTrue运行后adata.X会变成稠密矩阵内存可能暴涨。对于大矩阵建议显式设置max_value10截断并且分析完马上把adata.X重新赋值为稀疏矩阵。4.3 高变基因HVG与 PCA关键参数的取舍高变基因这一步Seurat 默认的FindVariableFeatures是选 2000 个。scanpy 里sc.pp.highly_variable_genes默认是flavorseurat同样是选 2000 个所以两边流程能对得上。sc.pp.highly_variable_genes(adata, n_top_genes2000, flavorseurat) adata.var[highly_variable].value_counts() # 后续 PCA 只用到高变基因 adata adata[:, adata.var.highly_variable].copy() sc.pp.scale(adata, max_value10) sc.tl.pca(adata, svd_solverarpack, n_comps50) sc.pl.pca_variance_ratio(adata, n_pcs50)关于svd_solver我建议你用arpack而不是randomized。randomized在超大矩阵时更快但很多版本下结果不够稳定PCA 结果在不同次运行之间会有细微差异影响下游聚类可重复性。arpack是精确求前几个主成分虽然慢一点但结果更可解释。主成分数量怎么选Seurat 教程里常见取 10-30但这跟数据本身高度相关。我的习惯是看pca_variance_ratio的拐点同时结合邻居图的稳定性综合判断而不是死扣某一个方差解释率阈值。实际项目里PBMC 10x 数据取 20 问题不大但如果是细胞异质性极高的肿瘤组织数据有时候要取到 40-50 才能保住稀有细胞群。4.4 邻居图、Leiden 聚类与 UMAP/tSNEscanpy 的聚类核心是sc.pp.neighbors和sc.tl.leidensc.pp.neighbors(adata, n_neighbors15, n_pcs30) sc.tl.umap(adata, min_dist0.5) sc.tl.leiden(adata, resolution1.0, key_addedleiden_1_0, flavorigraph, n_iterations2) sc.pl.umap(adata, colorleiden_1_0)这里有个跟 Seurat 很微妙的差异Seurat 里FindClusters的resolution参数和 scanpy 里的resolution虽然数值相近但内部算法实现不一样Seurat 用 Louvain 为主scanpy 用 Leiden 为主所以两边跑出的 cluster 数量和细胞分群结果不会完全一致。这不是 bug是算法差异。实际项目中如果想对结果做跨工具验证我更倾向于比较 marker gene 表达的模式是否一致而不是死磕 cluster 编号是否对得上。如果不指定flavorscanpy 的sc.tl.leiden在不同版本间默认值有变化。新版默认是 igraph 的 Leiden 实现速度更快但如果你装了旧版 scanpy可能默认走的是 leidenalg 的老接口。为了可复现代码里显式写flavorigraph这个参数是我强烈建议的。UMAP 的min_dist参数也是经验值。默认 0.5 偏稳健适合看大局调成 0.1 会分得更开适合检查稀有亚群但别调太低否则容易把连续分化的细胞强行切碎。我在多篇项目里实测下来0.3-0.5 是大多数转录组数据最靠谱的区间。4.5 marker gene 与细胞类型注释别只盯着 top marker聚类完成后的下一步就是注释细胞类型。scanpy 里找 marker 基因最常用的是sc.tl.rank_genes_groupssc.tl.rank_genes_groups(adata, leiden_1_0, methodwilcoxon, use_rawFalse) sc.pl.rank_genes_groups_heatmap(adata, n_genes20, groupbyleiden_1_0)注意use_rawFalse意味着直接使用adata.X做检验。如果你前面跑过sc.pp.scale那就不能用这个模式因为 scale 之后的表达值包含负数和截断做差异检验会产生大量假阳性。所以我在项目中的惯例是在归一化后、scale 之前保存一个scaledFalse的副本比如adata_clean adata.copy()。差异检验用这个副本。scale 之后的adata只用来做 PCA/UMAP/聚类。用methodwilcoxon还是methodt-testscanpy 默认是t-test但实际项目中我首选 wilcoxon因为它是非参数检验对单细胞数据常见的分布偏态和离群值更稳健。Seurat 里默认是 wilcoxon rank sum test两者逻辑类似。你若追求速度大数据集上methodt-test会明显快但代价是假阳性率偏高对下游 marker 筛选不友好。细胞注释这件事没有银弹。我的工作流是先跑一个 broad 级别的注释用经典 marker 把 T cell、B cell、Myeloid、NK、Epithelial 这些大类分开再用sc.tl.score_genes对细化的 signatures 做打分进一步细分亚群。# 以 CD8 T 细胞为例 tcell_markers [CD3D, CD3E] ctla4_markers [CTLA4, FOXP3, IL2RA] adata.obs[CD8_score] sc.tl.score_genes(adata, gene_list[CD8A, CD8B], score_nameCD8_score)这里踩过的坑是不同批次数据之间 marker 基因的检出率差异巨大。比如有的样本因为建库原因CD8A 表达整体偏低你不能只因为 CD8_score 低于某个绝对阈值就粗暴地把这群细胞注释成 CD4 T。这个时候需要看一下 CD4 的表达是不是也低如果都低那大概率是检测灵敏度问题而不是真的没有 CD8 T。实操中我一般会画一个 featureplot 把 CD4/CD8A/CD8B 同时投影到 UMAP 上用人眼看整体的空间分布而不是只看一个数字。5. 与 Seurat 流程的关键差异与迁移技巧如果你是 Seurat 老手直接看这一节就够了。我把最关键的差异点整理成一张对照表下面再补充几个我踩过坑的迁移细节。步骤Seurat 代码scanpy 代码注意点读取 10xRead10X()CreateSeuratObject()sc.read_10x_h5()/sc.read_10x_mtx()scanpy 不区分 assay直接是 AnnDataQCsubset(nFeature_RNA 200)sc.pp.filter_cells(adata, min_genes200)阈值需结合数据分布调整归一化NormalizeData(normalization.methodLogNormalize)sc.pp.normalize_totalsc.pp.log1p必须先备份原始 counts高变基因FindVariableFeatures(nfeatures2000)sc.pp.highly_variable_genes(n_top_genes2000)默认 flavor 就是 seuratPCARunPCA(npcs50)sc.tl.pca(adata, n_comps50)推荐svd_solverarpack聚类FindNeighborsFindClusters(resolution1.0)sc.pp.neighborssc.tl.leiden(resolution1.0)两者 cluster 编号不对应非线性降维RunUMAPsc.tl.umapmin_dist参数需多尝试markerFindAllMarkers()sc.tl.rank_genes_groups()注意use_raw的设置数据导出saveRDS()adata.write()保存为 h5ad建议保留原始数据备份5.1 代码迁移的几个易错点第一基因名大小写和特殊字符。human 数据里 Seurat 一般会自动把基因名转成首字母大写取决于你读入的 matrixscanpy 不会主动改。如果你的adata.var_names是全小写或带.比如 ENSEMBL ID后期跟外部注释文件合并时总是对不齐所以读入后最好统一成一种格式。adata.var_names [gene.upper() for gene in adata.var_names] # 转大写第二Seurat 的SCTransform没有完美平替。scanpy 生态里有类似的 normalize 方法比如sc.external.pp.scanorama做整合或者直接用sc.pp.normalize_total。但 SCT 在消除文库深度影响这件事上还是有自己的优势。如果一定要对比 SCT 的处理我建议用 scanpy 跑完标准流程后再用 Seurat 的 SCT 流程跑一遍比较 marker 基因的稳定性而不是强行在 scanpy 里模拟出 100% 一样的结果。第三对象大小与内存管理。AnnData 在默认情况下很多操作会原地修改对象inplaceTrue这跟 R 的 copy-on-modify 语义完全不同。想保留中间结果就.copy()否则你后面会发现跑着跑着原始数据没了。内存方面如果你在sc.pp.scale之后发现内存爆了赶紧把zero_centerFalse试试或者只对高变基因做 scale 而不是全基因 scale。5.2 批次整合的选型建议Harmony 还是 BBKNN批次整合是单细胞分析里的一个大头。Seurat 那边常用IntegrateDataCCA/MNNscanpy 生态里有几种方案harmony通过 harmonypy、bbknn、scVI以及外部的scanorama。我的选型逻辑是这样的数据来源差异不大、只是不同批次/不同样本直接用 Harmony速度快参数少和 PCA embedding 配合融洽。数据来源差异大比如不同平台、不同物种用 scVI 这样的深度生成模型效果更好但训练时间取决于 GPU。就想快速跑个整合看看BBKNN 可以但本质是改邻居图不是校正数据本身后续差异分析时需要注意。import scanpy.external as sce sc.pp.pca(adata, n_comps50) sce.pp.harmony_integrate(adata, batch, max_iter_harmony20) # 整合后结果是 adata.obsm[X_pca_harmony] sc.pp.neighbors(adata, use_repX_pca_harmony)注意 Harmony 整合的前提是 PCA 已经跑完并且你输入的 batch 列在adata.obs中存在。harmony_integrate默认会生成一个新的X_pca_harmony表示之后邻居计算要显式指定use_repX_pca_harmony否则默认用的还是原来的 PCA等于白跑。6. 常见问题与排查技巧实录最后这一节我把实操中反复遇到的几个问题整理成速查表都是踩过坑之后才总结出来的经验。现象可能原因排查方向与解法leiden报错ModuleNotFoundError: No module named igraph环境里没装python-igraph或版本冲突重新用 mamba 安装python-igraphleidenalg别混用 pip 和 condasc.pp.neighbors跑得极慢数据矩阵是稠密矩阵或者n_neighbors设置太大检查adata.X的稀疏性用scipy.sparse.issparse()确认n_neighbors一般 10-20 就够UMAP 图杂乱无章、无分群用了未 scale 的数据跑 PCA或者内部邻居参数不合适确认已经sc.pp.scale检查n_pcs是否取太少/太多尝试调大min_dist和 Seurat 结果 cluster 数量不一致算法本身差异不是 bug比较 marker 表达一致性别强求 cluster ID 完全一致保存 h5ad 后再次读取时基因名乱码编码问题或 gene symbol 没去重用make_uniqueTrue统一编码为 UTF-8画小提琴图/特征图时中文乱码或字体警告matplotlib 字体问题配置 matplotlib 字体plt.rcParams[font.sans-serif] [Arial]sc.pp.scale后内存崩溃生成稠密矩阵数据量太大只对高变基因 scale或zero_centerFalse或分块处理这里再说一个我踩过最深的坑高变基因筛选时机。如果你读取的是多批次合并后的数据直接用全部基因做 HVG 再跑 PCA容易被高表达基因比如线粒体基因、核糖体基因主导。我之前处理一个肿瘤数据集时HVG 列表里全是核糖体蛋白基因RPL/RPS 家族PCA 之后所有细胞被核糖体表达差异拉开真正的免疫细胞群完全没有分开。后来我在 HVG 之前先剔除核糖体基因和线粒体基因再跑流程T 细胞亚群马上就分出来了。# 剔除核糖体和线粒体基因 ribo_genes adata.var_names.str.startswith((RPL, RPS)) mt_genes adata.var_names.str.startswith(MT-) adata adata[:, ~(ribo_genes | mt_genes)].copy()另外一个经验是单细胞分析流程一定要做成脚本而不是纯 notebook 手动点。我用 snakemake 把 QC、归一化、聚类、marker 全部串成流水线每个中间结果都输出 h5ad 存档。这样不仅方便复现排查问题时也能很快定位是哪一步出的问题。scanpy 在这方面比 Seurat 有天然优势因为所有中间结果都能很方便地存成 h5ad而且完全基于文本的脚本能让 diff 非常直观。最后再说一个关于 scanpy 版本的小建议尽量锁定 scanpy 版本别跟着最新版一路升。scanpy 还在快速迭代API 偶有变动有些外部工具比如 CellChat、monocle 的 python 接口对特定版本有依赖。我在生产环境里用的是 scanpy1.9.*如果是新项目可以考虑 1.10 或更高但团队已有的脚本如果跑得好好的就别轻易动。数据分析和软件工程一样最怕的就是“升级一时爽排查火葬场”。