企业官网建设流程全解析

两篇讲完 Seurat 之后，后台一直有人问同样的问题：R 的流程学会了，是不是还得学 Python 的 Scanpy？说实话，现在做单细胞转录组测序数据分析，Seurat 和 Scanpy 就是绕不开的两套工具，两侧都有自己的用户群。这篇文章作为系列的第三篇，我准备把同样的分析思路用 scanpy 完整重写一遍，包括环境搭建、数据读取、QC 过滤、聚类、marker 鉴定和细胞注释，以及很多人真正关心的问题——这两个工具到底差在哪，分析结果能不能对上，数据怎么互转。你要是已经跟着前两篇用 Seurat 跑通过一遍，这篇文章能帮你说清楚“Seurat 里那套逻辑在 Python 里对应的是谁”；要是你本来就是写 Python 的，这篇文章可以直接照着敲代码，跑通全流程。

1. 为什么在 Seurat 之后还要学 Scanpy：两个生态的差异与迁移价值

1.1 两者到底有什么区别

很多初学者会觉得 Seurat 和 Scanpy 是“同类工具选一个就行”，其实这个理解有点片面。Seurat 是 R 语言生态里的老牌主力，封装程度非常高，从读数据到出图，十几个函数就能跑完整条流程；Scanpy 则是 Python 生态里的事实标准，底层基于 AnnData、numpy、scipy 这些科学计算库，API 更细粒度，灵活度也更高。

这个差异造成的实际影响是：Seurat 适合快速上手、交互式探索，尤其是你还在熟悉单细胞分析逻辑的时候，它的“保姆式”封装能让你少走很多弯路。但一旦分析需求超出常规流程，比如自定义一个过滤规则、把数据接进某个深度学习模型、写一套可以反复跑的自动化 pipeline，Python 这边的优势就很明显了。单细胞转录组测序数据分析发展到今天，早已不是“跑完聚类就行”的阶段，下游的轨迹推断、细胞通讯、RNA velocity、空间转录组、单细胞大模型，绝大多数新工具的第一个发布版本都是 Python 环境，这就是你绕不开 Scanpy 的根本原因。

1.2 从 Seurat 迁移过来需要理解的三个思维变化

第一，对象模型不一样。Seurat 里你操作的是 Seurat 对象，meta.data、assay、reductions 都挂在同一个对象上；Scanpy 里你操作的是 AnnData 对象，观测（细胞）、变量（基因）、降维坐标、非结构化信息分别放在 obs、var、obsm、uns 这些槽位里，概念上并不难对应。

第二，行为模式不一样。Seurat 的函数大多是“方法”，比如NormalizeData(seurat_obj)，传入对象然后返回处理后的对象；Scanpy 的很多函数默认是原地修改，sc.pp.normalize_total(adata)之后 adata 本身就被改了，不需要再赋值。

第三，找函数的方式不一样。Seurat 里功能相近的函数可能分散在不同包里，而 Scanpy 的函数前缀很有规律：sc.pp是预处理，sc.tl是工具类分析，sc.pl是绘图。搞清楚这三个前缀，你基本就能自己在文档里找到想要的功能。

1.3 两个工具适用的场景对照

这个对照不是要你站队，而是说单细胞分析本身是一套方法论，语言和工具只是载体。两边都掌握之后，你拿到一个新数据集、一个新算法时，就不再被工具限制。

2. 环境准备：Python、Scanpy 安装与依赖避坑

2.1 创建独立环境，别往 base 里直接怼

开始之前先把 Python 环境准备好。我推荐用 conda 或者 mamba 新建一个独立环境，不要直接装在系统 Python 或 conda base 里，否则后面装不同版本依赖时很容易把环境搞坏。Python 版本建议 3.9 到 3.11，Scanpy 1.10.x 在这些版本下兼容性最好；Python 3.12 目前也能装，但个别依赖还没完全跟上，遇到问题排查成本高。

创建环境的命令很简单：

conda create -n scanpy_env python=3.10 -y conda activate scanpy_env

如果你的网络访问 conda 默认源比较慢，可以在创建环境前先配置国内镜像源，比如清华或阿里云的 conda 镜像，速度会快很多。这一步属于老生常谈，但确实能省下不少等待时间。

2.2 安装 Scanpy 与配套算法包

Scanpy 的主程序可以直接用 pip 安装：

pip install scanpy[leiden]

注意这个[leiden]，它会顺带把 leidenalg 和 python-igraph 装上，这两个是跑 Leiden 聚类算法必需的。如果你用的是 conda-forge 源，也可以这样装，更适合 Windows 用户：

conda install -c conda-forge scanpy python-igraph leidenalg

装完之后建议顺手把 jupyter 或者 jupyterlab 也装上，因为单细胞分析高度依赖交互式环境，边写边看结果比写完整脚本再跑要高效得多。

验证安装是否成功：

import scanpy as sc print(sc.__version__)

能正常输出版本号，说明环境基本没问题。如果在这里就报错，先别急着找 scanpy 的麻烦，九成是依赖包版本冲突。

2.3 安装过程中最常见的三个报错

第一个是 numba 和 llvmlite 版本冲突。Scanpy 的很多计算会调用 numba 加速，而 numba 对 llvmlite 的版本要求很严格，经常出现ImportError: cannot import name 'xxx' from 'numba'这类错误。解决办法是把 numba 和 llvmlite 绑定重装，比如：

pip uninstall numba llvmlite -y pip install "numba==0.58.1" "llvmlite==0.41.1"

第二个是 Windows 下装 igraph 失败。这个在老版本里很常见，新版本 python-igraph 已经提供预编译的 wheel，但我仍然建议 Windows 用户优先走 conda-forge 安装路线，它会把 igraph、leidenalg 的二进制依赖一起处理掉，少很多折腾。

第三个是 pip 和 conda 混用导致的环境混乱。尽量不要先 conda 装一半再用 pip 补，先 pip 报错又用 conda 覆盖，这种操作会让依赖关系变得非常不可控。我的习惯是：环境创建用 conda，包管理全部走 pip，如果一定要用 conda 装，就同一批包都用 conda 装。

3. 数据接入与分析对象：AnnData 结构与 10X 数据读取

3.1 AnnData 到底长什么样

接触 Scanpy 的第一道坎就是理解 AnnData。简单说，AnnData 是一个容器，里面主要装着四类东西：

• adata.X：表达量矩阵，行是细胞，列是基因，通常是稀疏矩阵，因为单细胞数据里绝大多数元素是 0。
• adata.obs：细胞注释信息，DataFrame 格式，每一行对应一个细胞，可以存批次、样本、聚类结果等。
• adata.var：基因注释信息，DataFrame 格式，每一行对应一个基因，可以存基因名、是否高变基因、线粒体基因标记等。
• adata.obsm：降维结果，比如 PCA 坐标存在adata.obsm['X_pca']，UMAP 坐标存在adata.obsm['X_umap']。

和 Seurat 对象对比就很容易理解：obs对应meta.data，var对应rownames的基因注释，obsm对应reductions。理解了这个映射关系，你其实已经掌握了 80% 的数据结构基础。

3.2 读取 10X 标准输出

绝大多数单细胞转录组测序数据分析项目拿到的原始数据是 10X Genomics 的 Cell Ranger 输出，常见有两种格式：filtered_feature_bc_matrix.h5和matrix.mtx三个文件那套。Scanpy 分别有对应的读取函数：

import scanpy as sc # 方式一：读 h5 文件 adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5") # 方式二：读 mtx 目录 adata = sc.read_10x_mtx("filtered_feature_bc_matrix/", var_names="gene_symbols", make_unique=True)

读进来之后，先做两件小事：adata.var_names_make_unique()，确保基因名唯一；再print(adata.shape)看一眼矩阵维度。正常的 PBMC 3K 示例数据读进来应该是约 2700 个细胞、两万多个基因。

还有一个容易忽略的点：read_10x_h5读进来的基因名可能是 Ensembl ID 而不是基因 symbol，要看你的原始文件是哪种。如果读进来发现基因名是一串ENSG开头的编号，后续画 marker 图时非常痛苦，最好先做 ID 转换。反过来，read_10x_mtx可以直接通过var_names="gene_symbols"参数读基因 symbol。

4. 从 QC 到聚类：Scanpy 核心流程与 Seurat 的参数对照

4.1 QC 指标计算与过滤：别一把梭哈

先算质控指标，这一步和 Seurat 里PercentageFeatureSet加subset的流程对应：

# 标记线粒体基因 adata.var["mt"] = adata.var_names.str.startswith("MT-") # 计算质控指标 sc.pp.calculate_qc_metrics( adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True )

运行之后，adata.obs里会出现total_counts、n_genes_by_counts、pct_counts_mt这三列，分别对应 Seurat 里的nCount_RNA、nFeature_RNA、percent.mt。

接下来先画图看分布，再定阈值：

sc.pl.violin(adata, keys=["total_counts", "n_genes_by_counts", "pct_counts_mt"], multi_panel=True)

小提琴图能直观看出哪些细胞明显偏离主体。常见的过滤逻辑是线粒体比例过高说明细胞状态差，基因数过低说明可能是空液滴，基因数过高则可能是双细胞或者多细胞。但阈值没有统一标准，需要结合数据情况看。以 PBMC 为例，我一般先用：

sc.pp.filter_cells(adata, min_genes=200) sc.pp.filter_genes(adata, min_cells=3) adata = adata[adata.obs.pct_counts_mt < 20, :].copy() adata = adata[adata.obs.n_genes_by_counts < 6000, :].copy()

这里有两个细节。第一，filter_genes的min_cells=3是把只在极少数细胞里表达的基因过滤掉，这个操作能显著减少稀疏噪声，但对数据量大、dropout 严重的样本别过滤太狠，不然会丢真实信号。第二，线粒体比例这个阈值，不同组织差异很大，血液样本 20% 可能合理，但某些代谢活跃的组织细胞本身线粒体比例就偏高，20% 反而会误杀太多，建议每次都在小提琴图上确认拐点。另外，如果你已经用 DoubletFinder 或 scrublet 去除了双细胞，这里的高基因数过滤可以适当放宽。

4.2 标准化、高变基因、PCA：每个函数的副作用要清楚

预处理这块要求精细。Scanpy 的步骤拆得很细，每一步都对应 Seurat 的一个环节：

# 1. 标准化：每个细胞 counts 总数归一化到 1 万 sc.pp.normalize_total(adata, target_sum=1e4) # 2. 对数化 sc.pp.log1p(adata) # 3. 找高变基因 sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat") # 4. 用高变基因做一个子集 adata = adata[:, adata.var.highly_variable].copy() # 5. 标准化到均值为 0 方差为 1 sc.pp.scale(adata, max_value=10) # 6. PCA sc.tl.pca(adata, svd_solver="arpack", n_comps=50)

target_sum=1e4对应 Seurat 里NormalizeData默认的 LogNormalize，意思是每个细胞的基因表达量之和都缩放到 1 万，再取对数。flavor="seurat"表示用 Seurat 的 vst 方法选高变基因，如果你更习惯 Cell Ranger 的方式，可以改成flavor="cell_ranger"；针对原始 counts 是 UMI 的数据，flavor="seurat_v3"也是个不错的选项。

这里有一个非常关键、新手经常踩的坑：sc.pp.scale会直接覆盖adata.X，把里面的 counts 变成标准化后的残差。如果你后面想重新看原始表达量，或者重新做标准化，原始数据已经没了。所以务必要在 scale 之前留一个 raw：

adata.raw = adata.copy()

adata.raw会保存当前的对象，后续画 marker 图时能从 raw 里取原始 counts，这是一个很好的习惯。

PCA 跑完之后，用sc.pl.pca_variance_ratio(adata, n_pcs=50)画碎石图，看前多少个主成分能解释大部分方差。PBMC 3K 这种数据通常取前 30 个主成分即可，没必要纠结太多。主成分个数会影响后面的聚类，取少了丢信息，取多了把技术噪声也纳入进来。

4.3 邻居图与聚类：resolution 才是靠手感调的东西

PCA 之后是建邻居图和聚类，Scanpy 对应 Seurat 的FindNeighbors和FindClusters：

# 建邻居图，n_pcs 就是刚才 PCA 保留的维度 sc.pp.neighbors(adata, n_neighbors=10, n_pcs=30) # UMAP 降维 sc.tl.umap(adata) # Leiden 聚类 sc.tl.leiden(adata, resolution=0.5, key_added="leiden")

n_neighbors默认是 10，小数据集可以调到 5 到 8，大数据集可以调到 15 到 30。resolution是聚类分辨率，直接影响分群数量，越大分得越细。这个参数没有一个绝对正确的值，经验上是看 UMAP 图上有没有把同一种细胞切碎、或者把明显不同的细胞糊在一起。PBMC 数据用 0.5 左右通常能得到比较合理的分群。

我个人的习惯是：先在 resolution=1.0 下粗跑一遍，看看大概有多少群，标记一下每个群的 marker，再根据 biology 判断哪些群需要合并或细分，最后用合适的 resolution 重新跑。很多教程直接告诉你“resolution=0.5”，但对真实数据来说，这个值真的需要你自己调。

Leiden 和 Louvain 都是社区发现算法，Scanpy 里默认用 Leiden。相比 Louvain，Leiden 能保证聚类结果内部连通性更好，不会出现“一个群看起来连在一起但实际是断开的”这种问题。如果你确实想用 louvain，sc.tl.louvain(adata, resolution=0.5)也能跑，但我推荐直接上 leiden。

4.4 降维可视化：让数据自己说话

聚类结束后，第一个动作永远是在 UMAP 上看分群：

sc.pl.umap(adata, color="leiden", legend_loc="on data")

legend_loc="on data"会把群标签直接标在对应簇上，方便看每个群的位置。接下来画几个关键 marker 验证分群质量：

sc.pl.umap( adata, color=["CD3D", "MS4A1", "LYZ", "NKG7"], cmap="Reds", ncols=4 )

CD3D 是 T 细胞 marker，MS4A1 是 B 细胞 marker，LYZ 是髓系 marker，NKG7 是 NK 细胞 marker。如果这些基因的表达在 UMAP 上形成清晰、互不重叠的区域，说明聚类结果和细胞类型基本对得上。这时候再去看每个 cluster 是不是一类细胞，心里就有数了。

5. Marker 基因鉴定与细胞注释：找到每个细胞群的身份

5.1 rank_genes_groups：怎么找出每个群的 marker

聚类只是把细胞分成了若干群，但这些群到底是什么细胞，还得靠差异表达来找 marker。Scanpy 里的函数是rank_genes_groups：

sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon", pts=True)

method="wilcoxon"用的是 Wilcoxon 秩和检验，也是单细胞差异表达分析里最常用的方法。它会把每个 cluster 的每个基因跟其他所有 cluster 做比较，输出 log fold change、调整后 p 值等统计量。pts=True表示同时计算每个基因在每个 cluster 中的表达比例，这个信息在过滤“只在少数细胞中表达的 marker”时很有用。

结果提取成 DataFrame：

result = sc.get.rank_genes_groups_df(adata, group=None) result.to_csv("markers_leiden.csv", index=False)

输出的表里大致包含group、names、logfoldchanges、pvals_adj、pct_nz_group、pct_nz_reference几列，对应 SeuratFindAllMarkers的结果。注意，这里logfoldchanges是自然对数变换后的差异倍数，不是 log2，如果要跟 Seurat 的数值对比，记得统一换算。

筛选 marker 时，我一般看重两个标准：pvals_adj < 0.05且logfoldchanges > 0.5。同时要求这个基因在目标 cluster 中的表达比例明显高于其他 cluster，这样可以避免选出那种只有零散几个细胞表达、偶然显著升高的基因。

5.2 细胞注释的实操路线

拿到每个群的 top marker 之后，最直接的做法是把已知 marker 基因画成 dotplot：

marker_genes = { "T cell": ["CD3D", "CD3E"], "B cell": ["MS4A1", "CD79A"], "Monocyte": ["LYZ", "CD14"], "NK": ["NKG7", "GNLY"], } sc.pl.dotplot(adata, var_names=marker_genes, groupby="leiden")

dotplot 里点的大小代表表达该基因的细胞比例，颜色深浅代表平均表达量。如果某个 cluster 的 dot 模式跟你预期的 marker 对上了，就可以给它注释细胞类型。

对于没有明确 marker 的群，可以用score_genes给每个细胞打一个“基因集得分”：

sc.tl.score_genes(adata, gene_list=["CD3D", "CD3E"], score_name="T_score") sc.pl.umap(adata, color="T_score")

这个得分的原理很简单，就是计算目标基因集的平均表达量减去随机对照基因集的平均表达量，用来判断细胞倾向于哪类。它虽然不是严谨的差异检验，但作为快速验证手段非常好用。

现在也有很多自动注释工具，比如 CellTypist、SingleR、scCATCH。我的经验是：先手动注释一遍，用已知 marker 建立对数据的直觉，再用自动注释工具做交叉验证。直接拿自动注释结果当最终结论，很容易被数据库偏好带偏，尤其是样本来自稀有细胞类型或特定组织时。

6. 用同样的数据跑一遍：Scanpy 与 Seurat 结果差异分析

6.1 同一份 PBMC 3K，两边的结果有哪些不一样

我拿 10X 官方的 PBMC 3K 示例数据分别用 Seurat 和 Scanpy 跑过一遍，结论是：整体分群结果高度一致，T、B、NK、单核细胞这些主要类型两边都能分出来，marker 基因的排名也基本吻合。但细节上确实有差异：

第一，聚类数量不一定完全一样。Seurat 默认的resolution=0.8和 Scanpy 常用的resolution=0.5如果不对齐，分群数量自然会差几个。这不代表谁对谁错，纯粹是参数差异。第二，少量边界细胞的分群归属可能不同，因为两个工具在邻居图构建、PCA 计算时的默认参数实现略有差异。第三，UMAP 坐标不要直接比较，用完全相同的数据跑出来的 UMAP 图，大概率坐标长得不一样，这是 UMAP 算法本身的随机性和实现差异导致的，你比较的应该是聚类标签和 marker 基因，而不是坐标位置。

这里也提醒一个新手容易踩的误区：很多人看到两边 UMAP 长得不完全一样，就以为某个工具算错了。实际上，就算是只用 Seurat，每次跑 UMAP 图都不可能完全一样，因为 UMAP 自带随机种子。只要细胞类型的注释结果一致，就不需要恐慌。

6.2 两个工具之间的数据互转

实际项目里，团队里可能有人习惯 R，有人习惯 Python，数据互转几乎是绕不开的需求。Scanpy 的h5ad格式现在基本是单细胞数据交换的事实标准，Seurat 那边也有成熟的转换路径。

保存 Scanpy 对象：

adata.write("pbmc3k_scanpy.h5ad")

想在 R 里读这个文件，可以用zellkonverter或sceasy。sceasy 的用法大概是这样：

library(sceasy) # h5ad 转 Seurat sceasy::convertFormat( "pbmc3k_scanpy.h5ad", from = "anndata", to = "seurat", outFile = "pbmc3k_seurat.rds" )

反过来，从 Seurat 对象转成 h5ad：

library(sceasy) sceasy::convertFormat( pbmc_seurat, from = "seurat", to = "anndata", outFile = "pbmc3k_from_seurat.h5ad" )

转换过程中最常见的问题是基因名重复、assay 名不匹配、稀疏矩阵格式不兼容。转换完建议马上print(dim(adata))确认维度没变，再看一眼adata.obs里的分组标签有没有完整带过来。我这边的经验是，最稳妥的做法是只转原始 counts 和必要的 meta，降维坐标和聚类结果可以在转换后重新跑一遍，不要迷信转换工具能把所有分析结果都无损带过去。

7. 实际跑数据时最容易踩的坑和调试经验

7.1 索引错位的坑：过滤之后记得检查 index

Scanpy 的过滤操作，比如adata = adata[adata.obs.pct_counts_mt < 20, :].copy()，本质上是按布尔条件筛选行。如果你连续做了多次过滤，又使用了 concat 合并多个样本，obs_names可能出现重复或者索引错乱。这个问题通常不会立刻报错，而是会在后面画图、分组时出现莫名其妙的结果，比如某个 cluster 的 marker 跟实际表达对不上。

解决方式很简单，在预处理一开始就做好两件事：

adata.obs_names_make_unique() adata.var_names_make_unique()

这两个函数保证细胞名和基因名全局唯一。凡是做过 subset、concat、merge 之后，也建议重新跑一遍，成本很低但能防住很多诡异 bug。

7.2 scale 覆盖原始数据的坑

前面已经强调过sc.pp.scale会覆盖adata.X，这里再单独拎出来讲，是因为我见过太多人在这一步栽跟头。有些版本 Scanpy 在scale之后会默认保留一份 raw，但把你自己的原始 counts 备份出来永远是最稳的：

adata.raw = adata.copy() sc.pp.scale(adata, max_value=10)

有了这份 raw，后面无论什么时候想看原始表达，都能用sc.pl.umap(adata, color="CD3D")这种写法自动从 raw 取值。如果没有备份，你只能重新读原始文件再跑一遍整个预处理，浪费的时间没法估量。

7.3 依赖与内存问题

Scanpy 对内存的优化总体做得不错，因为底层用的scipy.sparse稀疏矩阵，默认只存非零元素。但对于几十万细胞、几万个基因的数据，PCA 和邻居图构建仍然会很占内存。如果内存吃紧，我的习惯是用sc.pp.downsample_counts先做一轮下采样，确认流程没问题再上全量数据。

另外，环境稳定性方面，Linux 服务器跑 Scanpy 的体验最稳，Mac 也问题不大，Windows 上偶尔会遇到python-igraph、leidenalg的编译问题。如果 Windows 下实在装不上，我建议直接用 WSL 装一个 Linux 环境，很多依赖问题直接消失。

7.4 batch 效应处理不能少

如果你的数据包含多个样本、多个批次，千万别直接跳过整合那一步。Scanpy 提供了几个常用方案：

# 方案一：Harmony（需要先跑 PCA） sc.external.pp.harmony_integrate(adata, key="batch") sc.pp.neighbors(adata, use_rep="X_pca_harmony") # 方案二：BBKNN（直接用，不需要 harmony 步骤） import bbknn sc.pp.neighbors(adata, n_neighbors=10, use_rep="X_pca") bbknn.bbknn(adata, batch_key="batch")

Harmony 是目前最常用的选项，特点是快、稳、默认参数结果通常就很好。判断是否需要整合，最简单的办法是画 UMAP 时同时按样本和聚类着色，如果每个聚类里样本比例严重失衡，或者 UMAP 上明显按样本分成几块，就说明存在批次效应。

最后分享一个我现在的习惯：每次拿到新数据，我都先用 Scanpy 把从读取到注释的完整流程存成一个脚本，所有参数放在文件头部，换数据时只改路径和阈值，十分钟内就能拿到一版初步结果。这个习惯让我从重复劳动里省出了大量时间，也让我对每个参数的影响有了更深的感觉。单细胞转录组测序数据分析的复杂度不低，但把流程固化成自己的工具链之后，真正的精力就能用来解决生物学问题了。