scrna-seq-pipeline

scRNA-seq 流程（Scanpy）

概述

此技能使用 Scanpy 和 AnnData 在 Python 中生成可复现的 scRNA-seq 分析。涵盖标准工作流：加载计数 → QC 过滤 → 标准化/对数转换 → HVG → 缩放/PCA → 邻居/UMAP → Leiden → 差异表达，以及使用 Harmony 或 scVI 进行批次校正。默认参数遵循 Scanpy PBMC3k 教程的惯例，可作为大多数 10X Genomics 数据集的安全起点。

输出为单个 .h5ad 文件，其中所有中间嵌入、聚类标签和 DE 结果都存储在 AnnData 对象中。

用法

安装依赖项：

pip install "scanpy>=1.10" anndata leidenalg python-igraph harmonypy scvi-tools

端到端运行流程：

import scanpy as sc
import scanpy.external as sce

sc.settings.verbosity = 2
sc.settings.set_figure_params(dpi=100, facecolor="white")

# 1. Load
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
# Alternatives:
# adata = sc.read_10x_mtx("filtered_feature_bc_matrix/", var_names="gene_symbols", cache=True)
# adata = sc.read_h5ad("input.h5ad")
adata.var_names_make_unique()

# 2. QC
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)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata = adata[adata.obs.n_genes_by_counts < 5000, :]
adata = adata[adata.obs.pct_counts_mt < 20, :].copy()

# 3. Normalize + HVG
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata  # freeze log-normalized full gene set BEFORE HVG subset
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat")
adata = adata[:, adata.var.highly_variable].copy()

# 4. Scale + PCA + neighbors + UMAP + Leiden
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50, svd_solver="arpack")
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

# 5. Differential expression per cluster
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")

# 6. Save
adata.write("result.h5ad", compression="gzip")

批次校正 — Harmony

# Assumes adata.obs["batch"] exists and PCA is already computed.
sce.pp.harmony_integrate(adata, key="batch")  # writes adata.obsm["X_pca_harmony"]
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30, use_rep="X_pca_harmony")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

批次校正 — scVI

scVI 需要原始计数。请在标准化之前运行，或保留一个原始计数层。

import scvi

# adata_raw holds unnormalized integer counts; adata.obs["batch"] is the batch key.
scvi.model.SCVI.setup_anndata(adata_raw, batch_key="batch")
model = scvi.model.SCVI(adata_raw, n_latent=30)
model.train(max_epochs=400, early_stopping=True)
adata.obsm["X_scVI"] = model.get_latent_representation()

sc.pp.neighbors(adata, n_neighbors=15, use_rep="X_scVI")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

响应格式

• 先给出用户需要的命令或代码，再进行说明
• 结构为：确认输入 → 可运行的代码 → 关键参数说明 → 注意事项
• 每项任务提供一个完整可运行的示例；不要展示所有替代方案
• 代码注释保持精简且实用（说明用途，而非存在原因）
• 目标：50–100 行代码，附简短说明

核心概念

• AnnData 布局。 adata.X 是当前表达矩阵（细胞 × 基因）。adata.obs 存储每个细胞的元数据（QC 指标、聚类标签、批次）。adata.var 存储每个基因的元数据（HVG 标记、mt）。嵌入存储在 adata.obsm（X_pca、X_umap、X_scVI）中。邻居图和 DE 结果存储在 adata.uns 中。
• .raw 快照。 adata.raw = adata 会在 HVG 子集筛选之前保存完整的对数标准化矩阵。下游绘图函数（sc.pl.umap(..., color="GENE")、sc.pl.rank_genes_groups）会从 .raw 读取数据，因此即使将 adata 限制为 HVG，也可以可视化任意基因。
• QC 阈值取决于数据集。 对于类似人类 10X PBMC 的数据，n_genes ∈ [200, 5000] 和 pct_mt < 20 是合理的默认值。确定阈值前先检查分布（sc.pl.violin）。人类使用 MT-，小鼠使用 mt-。
• HVG 风格很重要。 对对数标准化数据使用 flavor="seurat"（默认值）。对原始计数使用 flavor="seurat_v3"——它要求计数为整数，用于对数数据时会报错或产生误导结果。
• 邻居参数。 n_neighbors 控制局部/全局结构（10–50；值越大，UMAP 越平滑）。n_pcs 应能捕捉 sc.pl.pca_variance_ratio 中的拐点（通常为 20–50）。Leiden 的 resolution 控制聚类粒度（0.2–1.5；值越大，聚类越多）。
• 批次校正的选择。 Harmony 速度快，在 PCA 上运行，适用于轻度至中度的批次效应。scVI 是深度生成模型，需要原始计数，并且使用 GPU 可提升速度；它能处理更强的技术变异和多模态设计。

快速参考

加载

输入	调用
10X HDF5	sc.read_10x_h5("filtered_feature_bc_matrix.h5")
10X MTX 目录	sc.read_10x_mtx("path/", var_names="gene_symbols", cache=True)
现有 AnnData	sc.read_h5ad("input.h5ad")
加载后始终执行	adata.var_names_make_unique()

关键参数范围

参数	默认值	常见范围	说明
filter_cells(min_genes=)	200	100–500	去除空液滴/低质量液滴
filter_genes(min_cells=)	3	3–10	去除很少表达的基因
n_genes_by_counts <	5000	2500–8000	上限可标记双细胞
pct_counts_mt <	20	5–25	取决于人类组织
normalize_total(target_sum=)	1e4	1e4	CP10k；标准设置
n_top_genes（HVG）	2000	1000–5000	HVG 越多 → 信号和噪声都越多
scale(max_value=)	10	10	截断极端 z 分数
pca(n_comps=)	50	30–100	为下游 n_pcs 选择拐点
neighbors(n_neighbors=)	15	10–50	值越大 → 流形越平滑
neighbors(n_pcs=)	30	20–50	≤ n_comps
leiden(resolution=)	0.5	0.2–1.5	值越大 → 聚类越多

DE 与检查

sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

# Top-N markers per cluster as a DataFrame
import pandas as pd
result = adata.uns["rank_genes_groups"]
groups = result["names"].dtype.names
markers = pd.DataFrame({g: result["names"][g][:20] for g in groups})

保存 / 加载

adata.write("result.h5ad", compression="gzip")
adata = sc.read_h5ad("result.h5ad")

常见错误

• 错误： 在保存前先筛选高变基因（HVG），再执行 adata.raw = adata 正确： 始终先在对数归一化后的完整矩阵上设置 .raw，再执行 adata = adata[:, adata.var.highly_variable].copy() 原因： 不设置 .raw，就无法在下游分析中绘制或评分非高变基因

• 错误： 在对数归一化数据上使用 flavor="seurat_v3" 正确： 根据数据状态选择匹配的 HVG 算法：seurat_v3 需要原始整数计数；seurat（默认值）需要对数归一化数据 原因： 算法选择错误会导致 HVG 排名毫无意义或报错，进而影响所有下游步骤

• 错误： 加载 10X 数据后跳过 adata.var_names_make_unique() 正确： 加载数据后立即调用 adata.var_names_make_unique() 原因： 重复的基因符号会导致索引、HVG 选择和差异表达（DE）分析出现隐蔽错误

• 错误： 在归一化 / 对数转换后的数据上运行 scVI 正确： 将原始计数 AnnData 输入 scVI，并通过 setup_anndata 设置 batch_key 原因： scVI 直接对计数建模；输入归一化数据会违背其分布假设，并产生错误的潜在空间

• 错误： 对小鼠线粒体基因使用 MT- 前缀（或对人类基因使用 mt-） 正确： 人类基因使用 MT-，小鼠基因使用 mt- 原因： 前缀错误会导致 pct_counts_mt 始终为零，使 QC 筛选形同虚设

• 错误： 对 AnnData 进行布尔切片后不调用 .copy() 正确： 筛选后始终使用 .copy()：adata = adata[mask, :].copy() 原因： 不使用 .copy()，结果将是一个视图；后续原地操作会发出警告或失败

参考资料

• Scanpy：Wolf 等，Genome Biol，2018，https://doi.org/10.1186/s13059-017-1382-0
• Scanpy 教程：https://scanpy.readthedocs.io/en/stable/tutorials.html
• scVI：Gayoso 等，Nat Biotechnol，2022，https://doi.org/10.1038/s41587-021-01206-w
• Harmony：Korsunsky 等，Nat Methods，2019，https://doi.org/10.1038/s41592-019-0619-0