企业官网建设流程全解析

讲真，最近在生信交流群里看到不少朋友把 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 python=3.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 -y

2.3 加速方案与国内安装体验

如果你在国内，conda 默认源速度会非常感人。可以换清华源或阿里源，配置~/.condarc，把default_channels和custom_channels都指到国内镜像。pip 这边同样可以用清华的 PyPI 镜像，这样下载 scanpy 和它那一堆依赖包的速度能快一个数量级。

我自己的习惯是：conda-forge 装重型依赖（igraph、leidenalg、scanpy 本身），pip 装更新较快的纯 python 包（如 harmonypy）。这样不容易卡在某个包上。

3. AnnData：scanpy 的核心数据结构

理解了 AnnData，scanpy 就学会了一半。它跟 Seurat 不一样，Seurat 的分析对象是Seurat对象，一个 S4 对象里塞了多套 assay；而 scanpy 的核心对象叫AnnData，设计思路更接近“带注释的数据矩阵”。

3.1 AnnData 结构逐层拆解

一个典型的 AnnData 对象长这样：

这里有一个新手容易懵的地方：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() # 存原始 counts

3.2 从 10x 数据读入 AnnData

scanpy 支持多种读取方式，最常用的是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_names="gene_symbols", make_unique=True)

注意read_10x_mtx的make_unique=True参数：如果基因名有重复（比如某些版本会同时输出基因符号和基因 ID），必须设置去重，否则后面.var_names重复会导致一系列报错。另外，读取后要检查一下adata维度是不是符合预期，常见的 cellranger 输出过滤后矩阵一般细胞数在几千到几万，基因数在 2 万左右。

3.3 和 Seurat 对象互转的实用经验

如果团队里 R 和 python 两边都有人，互转是绕不开的。scanpy 官方推荐的方式是通过sceasy库转成.rds：

import sceasy # anndata 转 Seurat（需要本地装了 R 的 Seurat） sceasy.convert(adata, to="seurat", outFile="adata.rds") # Seurat 转 anndata sceasy.convert(adata, to="anndata", outFile="adata.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_top=None, inplace=True) sc.pl.violin(adata, keys=["n_genes_by_counts", "total_counts", "pct_counts_mt"], multi_panel=True)

阈值怎么定？Seurat 经典教程里一般用 nFeature_RNA > 200 且 < 2500，percent.mt < 5%。但实际项目里千万别死板套用：

• 如果是组织样本（比如肿瘤组织），细胞的线粒体比例天然比细胞系高，8% 甚至 10% 都可能是正常细胞群。
• 如果是冷冻组织解离后的数据，基因数普遍偏低，阈值可以适当放宽。

我的做法是先画小提琴图看分布，再结合总 UMI 数的双峰分布来定。QC 过滤这步用 sc.pp.filter_cells 和 sc.pp.filter_genes 即可：

# 先粗略过滤，pac 掉明显低质量的 sc.pp.filter_cells(adata, min_genes=200) sc.pp.filter_genes(adata, min_cells=3) # 再根据 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_rate=0.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）与 log1p

scanpy 的数据归一化流程和 Seurat 有概念上的差异。Seurat 的NormalizeData默认是 log1p(CPM/100)，也就是 log1p(normalize total=1e4)。scanpy 里对应的操作是：

sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata)

这两行代码合起来就等于 Seurat 的默认NormalizeData。但scanpy 不会把归一化后的数据存到scale.data，它直接把adata.X给覆盖了。所以我前面才强调，必须在归一化之前先把原始 counts 存到adata.layers["counts"]，否则你想用 MAST 或 wilcox 做差异分析时，没有原始 counts 会很被动。

如果你之后要做sc.pp.scale（Z-score 标准化，目的是让每个基因的均值 0、方差 1），注意 scanpy 的scale默认是zero_center=True，运行后adata.X会变成稠密矩阵，内存可能暴涨。对于大矩阵，建议显式设置max_value=10截断，并且分析完马上把adata.X重新赋值为稀疏矩阵。

4.3 高变基因（HVG）与 PCA：关键参数的取舍

高变基因这一步，Seurat 默认的FindVariableFeatures是选 2000 个。scanpy 里sc.pp.highly_variable_genes默认是flavor="seurat"，同样是选 2000 个，所以两边流程能对得上。

sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat") adata.var["highly_variable"].value_counts() # 后续 PCA 只用到高变基因 adata = adata[:, adata.var.highly_variable].copy() sc.pp.scale(adata, max_value=10) sc.tl.pca(adata, svd_solver="arpack", n_comps=50) sc.pl.pca_variance_ratio(adata, n_pcs=50)

关于svd_solver，我建议你用arpack而不是randomized。randomized在超大矩阵时更快，但很多版本下结果不够稳定，PCA 结果在不同次运行之间会有细微差异，影响下游聚类可重复性。arpack是精确求前几个主成分，虽然慢一点，但结果更可解释。

主成分数量怎么选？Seurat 教程里常见取 10-30，但这跟数据本身高度相关。我的习惯是看pca_variance_ratio的拐点，同时结合邻居图的稳定性综合判断，而不是死扣某一个方差解释率阈值。实际项目里，PBMC 10x 数据取 20 问题不大，但如果是细胞异质性极高的肿瘤组织数据，有时候要取到 40-50 才能保住稀有细胞群。

4.4 邻居图、Leiden 聚类与 UMAP/tSNE

scanpy 的聚类核心是sc.pp.neighbors和sc.tl.leiden：

sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30) sc.tl.umap(adata, min_dist=0.5) sc.tl.leiden(adata, resolution=1.0, key_added="leiden_1_0", flavor="igraph", n_iterations=2) sc.pl.umap(adata, color="leiden_1_0")

这里有个跟 Seurat 很微妙的差异：Seurat 里FindClusters的resolution参数和 scanpy 里的resolution虽然数值相近，但内部算法实现不一样（Seurat 用 Louvain 为主，scanpy 用 Leiden 为主），所以两边跑出的 cluster 数量和细胞分群结果不会完全一致。这不是 bug，是算法差异。实际项目中如果想对结果做跨工具验证，我更倾向于比较 marker gene 表达的模式是否一致，而不是死磕 cluster 编号是否对得上。

如果不指定flavor，scanpy 的sc.tl.leiden在不同版本间默认值有变化。新版默认是 igraph 的 Leiden 实现，速度更快，但如果你装了旧版 scanpy，可能默认走的是 leidenalg 的老接口。为了可复现，代码里显式写flavor="igraph"这个参数是我强烈建议的。

UMAP 的min_dist参数也是经验值。默认 0.5 偏稳健，适合看大局；调成 0.1 会分得更开，适合检查稀有亚群；但别调太低，否则容易把连续分化的细胞强行切碎。我在多篇项目里实测下来，0.3-0.5 是大多数转录组数据最靠谱的区间。

4.5 marker gene 与细胞类型注释：别只盯着 top marker

聚类完成后的下一步就是注释细胞类型。scanpy 里找 marker 基因最常用的是sc.tl.rank_genes_groups：

sc.tl.rank_genes_groups(adata, "leiden_1_0", method="wilcoxon", use_raw=False) sc.pl.rank_genes_groups_heatmap(adata, n_genes=20, groupby="leiden_1_0")

注意，use_raw=False意味着直接使用adata.X做检验。如果你前面跑过sc.pp.scale，那就不能用这个模式，因为 scale 之后的表达值包含负数和截断，做差异检验会产生大量假阳性。所以我在项目中的惯例是：

• 在归一化后、scale 之前保存一个scaled=False的副本，比如adata_clean = adata.copy()。差异检验用这个副本。
• scale 之后的adata只用来做 PCA/UMAP/聚类。

用method="wilcoxon"还是method="t-test"？scanpy 默认是t-test，但实际项目中我首选 wilcoxon，因为它是非参数检验，对单细胞数据常见的分布偏态和离群值更稳健。Seurat 里默认是 wilcoxon rank sum test，两者逻辑类似。你若追求速度，大数据集上method="t-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_name="CD8_score")

这里踩过的坑是：不同批次数据之间 marker 基因的检出率差异巨大。比如有的样本因为建库原因，CD8A 表达整体偏低，你不能只因为 CD8_score 低于某个绝对阈值，就粗暴地把这群细胞注释成 CD4 T。这个时候需要看一下 CD4 的表达是不是也低，如果都低，那大概率是检测灵敏度问题，而不是真的没有 CD8 T。实操中我一般会画一个 featureplot 把 CD4/CD8A/CD8B 同时投影到 UMAP 上，用人眼看整体的空间分布，而不是只看一个数字。

5. 与 Seurat 流程的关键差异与迁移技巧

如果你是 Seurat 老手，直接看这一节就够了。我把最关键的差异点整理成一张对照表，下面再补充几个我踩过坑的迁移细节。

5.1 代码迁移的几个易错点

第一，基因名大小写和特殊字符。human 数据里 Seurat 一般会自动把基因名转成首字母大写（取决于你读入的 matrix），scanpy 不会主动改。如果你的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 在默认情况下很多操作会原地修改对象（inplace=True），这跟 R 的 copy-on-modify 语义完全不同。想保留中间结果就.copy()，否则你后面会发现跑着跑着原始数据没了。内存方面，如果你在sc.pp.scale之后发现内存爆了，赶紧把zero_center=False试试，或者只对高变基因做 scale 而不是全基因 scale。

5.2 批次整合的选型建议：Harmony 还是 BBKNN？

批次整合是单细胞分析里的一个大头。Seurat 那边常用IntegrateData（CCA/MNN），scanpy 生态里有几种方案：harmony（通过 harmonypy）、bbknn、scVI，以及外部的scanorama。

我的选型逻辑是这样的：

• 数据来源差异不大、只是不同批次/不同样本：直接用 Harmony，速度快，参数少，和 PCA embedding 配合融洽。
• 数据来源差异大（比如不同平台、不同物种）：用 scVI 这样的深度生成模型效果更好，但训练时间取决于 GPU。
• 就想快速跑个整合看看：BBKNN 可以，但本质是改邻居图，不是校正数据本身，后续差异分析时需要注意。

import scanpy.external as sce sc.pp.pca(adata, n_comps=50) sce.pp.harmony_integrate(adata, "batch", max_iter_harmony=20) # 整合后结果是 adata.obsm["X_pca_harmony"] sc.pp.neighbors(adata, use_rep="X_pca_harmony")

注意 Harmony 整合的前提是 PCA 已经跑完，并且你输入的 batch 列在adata.obs中存在。harmony_integrate默认会生成一个新的X_pca_harmony表示，之后邻居计算要显式指定use_rep="X_pca_harmony"，否则默认用的还是原来的 PCA，等于白跑。

6. 常见问题与排查技巧实录

最后这一节，我把实操中反复遇到的几个问题整理成速查表，都是踩过坑之后才总结出来的经验。

这里再说一个我踩过最深的坑：高变基因筛选时机。如果你读取的是多批次合并后的数据，直接用全部基因做 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 接口）对特定版本有依赖。我在生产环境里用的是 scanpy==1.9.*，如果是新项目可以考虑 1.10 或更高，但团队已有的脚本如果跑得好好的，就别轻易动。数据分析和软件工程一样，最怕的就是“升级一时爽，排查火葬场”。