📚 课程目录第2章 Scanpy入门

Scanpy 入门：安装、Anndata 与预处理流程

好，我们直接进入正题。这一章我会带你从零搭好 Scanpy 环境，再拆开 Anndata 这个核心数据结构，最后走一遍预处理管线。我刚开始搞 scRNA-seq 分析的那会儿，光装包就踩了三个坑，后面你们看得到。😏

1. 安装 Scanpy 与运行环境

Scanpy 的安装本身不复杂，但依赖链条有点长。我个人习惯用 conda 创建独立环境，因为之前用 pip 装的时候，NumPy 和 SciPy 版本冲突搞得我差点重装系统（夸张了，但确实折腾了一下午）。

conda create -n scanpy_env python=3.9
conda activate scanpy_env
pip install scanpy[leiden]

小提示：[leiden] 这个 extra 会帮你装上 leidenalg 和 igraph，跑聚类用。如果你在 M1/M2 Mac 上，建议用 conda install -c conda-forge leidenalg 替代。

装完验证一下：

python -c "import scanpy as sc; print(sc.__version__)"

嗯，看到版本号输出就成功了。我记得第一次跑这条命令的时候，终端卡了十秒，我还以为死机了……

注意：别在 base 环境里直接 pip install scanpy。我之前图省事这么干过，结果跟系统中已有的 R 包产生符号链接冲突，报了一堆 “ImportError: libtbb.so.2”。老老实实建虚拟环境。

2. 理解 Anndata：单细胞数据的“集装箱”

Scanpy 的所有数据都装在一个叫 AnnData 的对象里。说白了，它就是一个能同时存矩阵、观测信息、特征信息和非结构化元数据的容器。

属性名	维度 / 类型	存什么？
`.X`	cells × genes	表达矩阵（通常稀疏矩阵）
`.obs`	cells × 元数据	细胞注释、批次、样本ID等
`.var`	genes × 元数据	基因名、染色体、高变基因标记等
`.obsm`	cells × 嵌入坐标	PCA、UMAP、t-SNE 坐标
`.varm`	genes × 特征载荷	PCA 载荷、基因权重
`.uns`	dict	自定义元数据、调色板、聚类树
`.layers`	dict of matrices	归一化前、后、校正后等多层表达

你想想看，一个对象就能装下所有分析中间结果，不用来回写 CSV。我在项目里就吃过纯用 CSV 的亏：20 个样本的 UMAP 坐标和聚类标签经常对不上号，后来全改用 h5ad 格式，一步到位。

3. 数据加载与预处理流程

拿到原始表达矩阵后，预处理一般按这个顺序走。我把它画成一张流程图，方便你对照着看。

流程看起来很线性，但有几个关键点一定要搞明白。

3.1 读取数据

最常见的是 10x Genomics 的 HDF5 格式。直接一行搞定：

ada = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
ada.var_names_make_unique()  # 防止重复基因名

如果你拿到的是 mtx 文件（三个文件），可以用：

ada = sc.read_10x_mtx("path/to/matrix.mtx", var_names='gene_symbols')

  避坑指南：我曾经遇到读进来的 .X 是普通 ndarray，占用几十 GB 内存。建议后续手动转稀疏：ada.X = scipy.sparse.csr_matrix(ada.X)。

3.2 质量控制

我习惯先算几个核心指标：

sc.pp.calculate_qc_metrics(ada, percent_top=None, log1p=False, inplace=True)

然后根据库大小 / 基因数 / 线粒体比例过滤。看你的数据分布，通常这样写：

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

别一刀切：不同组织、不同测序深度，阈值差很多。我做过一个胰腺数据，线粒体比例在 15% 左右才是正常，按 20% 切反而把好细胞滤掉了。先画个小提琴图看看分布。

3.3 标准化与对数变换

标准做法是文库大小归一化 + log1p：

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

这一步把表达值拉到类似正态分布，方便后续 PCA 和聚类。其实 log1p 本质是 ln(1+x)，你想想看，零表达不会变成负无穷，保留了稀疏特性。

3.4 高变基因选择

选 2000 个左右的高变基因，能大大降低噪音和计算量：

sc.pp.highly_variable_genes(ada, n_top_genes=2000, flavor='seurat_v3')
ada = ada[:, ada.var.highly_variable].copy()

个人经验：flavor 参数默认是 'seurat'，如果数据比较脏（比如 dropout 率高），试试 'seurat_v3' 会更稳定。我曾经在 PBMC 数据上对比过，'seurat_v3' 选出的基因在下游差异分析中富集更好。

4. 完整预处理代码模板

把上面串起来，就得到一个可以直接用的函数：

import scanpy as sc
import scipy.sparse

def preprocess_scrna(path, sample_id="sample"):
    # 读取
    ada = sc.read_10x_h5(path)
    ada.var_names_make_unique()
    ada.obs["sample"] = sample_id

    # 转稀疏
    if not isinstance(ada.X, scipy.sparse.csr_matrix):
        ada.X = scipy.sparse.csr_matrix(ada.X)

    # QC
    ada.var["mt"] = ada.var_names.str.startswith("MT-")
    sc.pp.calculate_qc_metrics(ada, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)
    sc.pp.filter_cells(ada, min_genes=200)
    sc.pp.filter_genes(ada, min_cells=3)
    ada = ada[ada.obs.pct_counts_mt < 20, :].copy()

    # 标准化 + 对数
    sc.pp.normalize_total(ada, target_sum=1e4)
    sc.pp.log1p(ada)

    # 高变基因
    sc.pp.highly_variable_genes(ada, n_top_genes=2000, flavor='seurat_v3')
    ada = ada[:, ada.var.highly_variable].copy()

    print(f"预处理完成：{ada.n_obs} 细胞，{ada.n_vars} 基因")
    return ada

注意：这段代码默认线粒体基因前缀是 MT-（人类）。如果你做小鼠，记得改成 mt-。我之前没改，跑了半天发现 QC 全废了。

嗯，到这里你已经能完整地把 10x 原始数据读进来，做完标准化和高变基因选择，得到一个干净的 AnnData 对象。下一件事就是降维 + 聚类，但我们先让这一步走通。记住：预处理的质量直接决定后面批次校正的结果，千万别凑合。