scanpy-scrna-seq

Scanpy 单细胞 RNA-seq 分析

概览

Scanpy 是基于 AnnData 格式构建的可扩展 Python 工具包，用于分析单细胞 RNA-seq 数据。本技能涵盖端到端标准工作流：质量控制、归一化、高变基因筛选、降维、聚类、标记基因识别和细胞类型注释。它会生成带注释的数据集和达到发表质量的可视化结果。

适用场景

• 分析单细胞 RNA-seq 计数矩阵（10X Genomics、h5ad、CSV、loom）
• 对 scRNA-seq 数据集进行质量控制筛选（线粒体比例、基因计数）
• 执行降维：PCA、UMAP、t-SNE
• 通过 Leiden 社区检测识别细胞簇
• 查找每个簇的差异表达标记基因（Wilcoxon、t 检验、逻辑回归）
• 根据已知标记基因面板注释细胞类型
• 进行轨迹推断和伪时间分析（PAGA、扩散伪时间）
• 生成达到发表质量的单细胞图表（点图、热图、堆叠小提琴图）
• 比较不同实验条件下各细胞类型中的基因表达
• 对导出的表格绘制通用图表（火山图／热图／柱状图／KM），请使用 omics-plotting SKILL；单细胞视图仍使用 scanpy sc.pl.*
• 若要在现有 R 工作流中分析 scRNA-seq，或需要 Seurat 专属的 assay 类型，请改用 Seurat（R/Bioconductor）

前提条件

• Python 包：scanpy>=1.10、leidenalg、igraph、anndata
• 数据要求：基因表达计数矩阵（细胞 x 基因）。常见格式：10X Cell Ranger 输出（filtered_feature_bc_matrix/）、.h5ad、.h5、.csv
• 环境：Python 3.9+；细胞数超过 50k 时建议使用 16GB+ RAM

pip install "scanpy[leiden]" anndata

预先访谈

编写任何分析代码之前，先与用户确认以下事项。

decisions:
  - id: D1
    param: pct_counts_mt
    kind: required
    source: data
    ask: "Dying cells leak cytoplasmic RNA and read as high-mitochondrial - where should the cut sit for this tissue?"
    default: "20% (PBMC ~5%, solid tumour ~20%; propose from the observed distribution)"

  - id: D2
    param: batch_key
    kind: required
    source: data
    ask: "Which variable separates samples that were processed apart (donor, run, 10X lane), so its effect is not read as biology?"
    default: null
    skip_if: "all cells come from one sample and one processing run"

  - id: D3
    param: resolution
    kind: optional
    source: user
    ask: "How finely should cells be split - broad lineages, or subtypes within them?"
    default: "0.8 (0.1-0.3 broad lineages, 1.0-2.0 fine subtypes)"

  - id: D4
    param: n_top_genes
    kind: optional
    source: user
    ask: "How many variable genes should drive the embedding? More captures subtle states but adds noise."
    default: 2000

  - id: D5
    param: n_pcs
    kind: optional_conditional
    source: data
    depends_on: [D4]
    ask: "How many principal components carry real structure, by the variance-ratio elbow?"
    default: 40

  - id: D6
    param: min_genes, min_cells
    kind: optional_conditional
    source: data
    ask: "Where should empty droplets and undetected genes be cut off?"
    default: "cells with <200 genes, genes in <3 cells"

  - id: D6b
    param: annotationStrategy
    kind: required
    source: user
    ask: "After clustering, should clusters be named from canonical markers by hand, or handed to a reference-based annotator?"
    default: "manual marker-based within this skill; see `single-cell-annotation-guide` to choose"

  - id: D6c
    param: markerPanel
    kind: required
    source: literature
    depends_on: [D6b]
    ask: "Which canonical markers define the populations expected in this tissue and state?"
    default: null
    skip_if: "annotation delegated to a reference-based tool"

  - id: D7
    param: method (rank_genes_groups)
    kind: optional
    source: user
    ask: "Which test should rank marker genes per cluster?"
    default: "wilcoxon"

  - id: D8
    param: target_sum, flavor
    kind: optional
    source: user
    ask: "Should the standard library-size normalization and Seurat v3 highly-variable-gene selection be used?"
    default: "target_sum 1e4, flavor seurat_v3"

  - id: D9
    param: random_state
    kind: never_ask
    source: data
    reason: "Fixes the draw for reproducibility; does not change what the data supports"
    default: 0

D5 取决于 D4，因为肘部位置是根据选定的变异基因计算出的 PCA 来判断的——基因数量不同，肘部位置也会变化。D2 会影响整个流程，而不仅仅是单次调用：指定批次后，必须在邻居图构建之前增加校正阶段（harmony-batch-correction 或 ComBat）。

工作流

步骤 1：设置和数据加载

配置 scanpy 设置，并将计数矩阵读入 AnnData 对象。确保基因名称唯一。

import scanpy as sc
import numpy as np
import pandas as pd

# Configure settings
sc.settings.verbosity = 3  # errors=0, warnings=1, info=2, hints=3
sc.settings.set_figure_params(dpi=80, facecolor="white")

# Load data — pick the reader matching your format
adata = sc.read_10x_mtx(
    "path/to/filtered_feature_bc_matrix/",
    var_names="gene_symbols",
    cache=True,
)
# Alternatives:
# adata = sc.read_h5ad("data.h5ad")
# adata = sc.read_10x_h5("data.h5")

adata.var_names_make_unique()
print(f"Loaded: {adata.n_obs} cells x {adata.n_vars} genes")
print(f"Sparsity: {1 - adata.X.nnz / (adata.n_obs * adata.n_vars):.1%}")

步骤 2：质量控制

标注线粒体／核糖体基因，计算质量控制指标，可视化其分布，并筛选低质量细胞。

# Annotate gene groups
adata.var["mt"] = adata.var_names.str.startswith("MT-")    # human; use "mt-" for mouse
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))

# Calculate QC metrics
sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt", "ribo"], percent_top=None, log1p=False, inplace=True
)

# Visualize QC distributions
sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True,
)
sc.pl.scatter(adata, x="total_counts", y="pct_counts_mt")
sc.pl.scatter(adata, x="total_counts", y="n_genes_by_counts")

# Filter — adjust thresholds based on the QC plots above
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, :]  # remove potential doublets
adata = adata[adata.obs.pct_counts_mt < 20, :]         # remove dying cells

print(f"After QC: {adata.n_obs} cells x {adata.n_vars} genes")

步骤 3：归一化和特征筛选

归一化文库大小，进行对数转换，并识别高变基因（HVG）。保存原始计数，以便后续进行差异表达分析。

# Preserve raw counts in a layer
adata.layers["counts"] = adata.X.copy()

# Normalize to 10,000 counts per cell, then log-transform
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# Save normalized data as .raw for visualization
adata.raw = adata

# Identify highly variable genes
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat_v3", layer="counts")
sc.pl.highly_variable_genes(adata)

# Subset to HVGs
adata = adata[:, adata.var.highly_variable]
print(f"HVG subset: {adata.n_obs} cells x {adata.n_vars} genes")

步骤 4：缩放和回归

回归去除混杂因素，并缩放基因表达值，为 PCA 做准备。

# Regress out unwanted sources of variation
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])

# Scale to unit variance, clip extreme values
sc.pp.scale(adata, max_value=10)

print(f"Scaled matrix: mean={adata.X.mean():.4f}, std={adata.X.std():.4f}")

步骤 5：降维

计算 PCA，查看方差比肘部图，然后构建邻域图并使用 UMAP 进行嵌入。

# PCA
sc.tl.pca(adata, svd_solver="arpack", n_comps=50)
sc.pl.pca_variance_ratio(adata, log=True, n_pcs=50)

# Determine n_pcs from elbow plot (typically 30-50)
n_pcs = 40

# Compute k-nearest neighbor graph
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=n_pcs)

# UMAP embedding
sc.tl.umap(adata)
sc.pl.umap(adata, color=["n_genes_by_counts", "total_counts", "pct_counts_mt"])

print(f"UMAP shape: {adata.obsm['X_umap'].shape}")

步骤 6：聚类

使用多个分辨率运行 Leiden 聚类，并通过目视检查选出最佳粒度。

# Test multiple resolutions
for res in [0.3, 0.5, 0.8, 1.0, 1.2]:
    sc.tl.leiden(adata, resolution=res, key_added=f"leiden_res{res}", flavor="igraph", n_iterations=2)

# Compare resolutions side-by-side
sc.pl.umap(
    adata,
    color=["leiden_res0.3", "leiden_res0.5", "leiden_res0.8", "leiden_res1.0"],
    ncols=2,
    legend_loc="on data",
)

# Pick final resolution based on biological plausibility
adata.obs["leiden"] = adata.obs["leiden_res0.5"]
n_clusters = adata.obs["leiden"].nunique()
print(f"Selected resolution=0.5: {n_clusters} clusters")

步骤 7：标记基因识别

使用原始计数，通过 Wilcoxon 秩和检验查找每个簇的差异表达基因。

# Rank genes per cluster
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon", use_raw=True)

# Visualization
sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False)
sc.pl.rank_genes_groups_dotplot(adata, n_genes=5)
sc.pl.rank_genes_groups_heatmap(adata, n_genes=10, use_raw=True, show_gene_labels=True)

# Extract as DataFrame
markers_df = sc.get.rank_genes_groups_df(adata, group=None)
top_markers = markers_df[markers_df["pvals_adj"] < 0.05].groupby("group").head(10)
print(f"Significant markers: {len(top_markers)} genes across {top_markers['group'].nunique()} clusters")
print(top_markers[["group", "names", "logfoldchanges", "pvals_adj"]].head(15))

步骤 8：细胞类型注释和导出

根据经典标记基因分配细胞类型，进行可视化，并保存所有结果。

# Define canonical markers (example: human PBMC)
marker_genes = {
    "CD4 T cells":      ["CD3D", "CD3E", "IL7R"],
    "CD8 T cells":      ["CD8A", "CD8B", "GZMK"],
    "B cells":          ["MS4A1", "CD79A", "CD79B"],
    "NK cells":         ["NKG7", "GNLY", "KLRD1"],
    "CD14+ Monocytes":  ["CD14", "LYZ", "S100A9"],
    "FCGR3A+ Monocytes":["FCGR3A", "MS4A7"],
    "Dendritic cells":  ["FCER1A", "CST3"],
    "Platelets":        ["PPBP", "PF4"],
}

# Dot plot — rows = clusters, columns = markers
sc.pl.dotplot(adata, var_names=marker_genes, groupby="leiden", use_raw=True)

# Map clusters → cell types (adjust based on your dot plot)
cluster_to_celltype = {
    "0": "CD4 T cells",
    "1": "CD14+ Monocytes",
    "2": "B cells",
    "3": "CD8 T cells",
    "4": "NK cells",
    "5": "FCGR3A+ Monocytes",
    "6": "Dendritic cells",
    "7": "Platelets",
}
adata.obs["cell_type"] = adata.obs["leiden"].map(cluster_to_celltype).fillna("Unknown")
sc.pl.umap(adata, color="cell_type", legend_loc="on data", frameon=False)

# Save results
adata.write("results/annotated_data.h5ad", compression="gzip")
adata.obs.to_csv("results/cell_metadata.csv")
markers_df.to_csv("results/marker_genes.csv", index=False)
print(f"Saved: {adata.n_obs} cells with {adata.obs['cell_type'].nunique()} cell types")

关键参数

参数	默认值	范围 / 选项	作用
min_genes (filter_cells)	200	100-500	每个细胞检测到的最少基因数；值越低，保留的细胞越多
min_cells (filter_genes)	3	3-10	表达某个基因的最少细胞数；值越高，筛选越严格
pct_counts_mt 阈值	20	5-30	线粒体比例上限；取决于组织类型（PBMC 为 5%，肿瘤为 20%）
n_top_genes (HVG)	2000	1000-5000	高度可变基因的数量；数量越多，捕捉到的细微变化越多
flavor (HVG)	"seurat_v3"	"seurat"、"seurat_v3"、"cell_ranger"	HVG 检测算法
n_pcs (neighbors)	40	20-50	邻居图使用的 PC 数量；请查看肘部图
n_neighbors	15	5-50	k-NN 图中的 k；值越低，越强调局部结构
resolution (leiden)	0.5	0.1-2.0	聚类粒度；值越高 → 聚类数量越多
method (rank_genes)	"wilcoxon"	"wilcoxon"、"t-test"、"logreg"	DE 检验；推荐在发表研究中使用 Wilcoxon
target_sum (normalize)	1e4	1e4-1e5	每个细胞的目标文库大小

常见方案

方案：使用 ComBat 进行批次校正

适用场景：多个样本或批次在 UMAP 上呈现明显的批次效应。

# Requires 'batch' column in adata.obs
sc.pp.combat(adata, key="batch")
# Re-run PCA, neighbors, UMAP, clustering after correction
sc.tl.pca(adata, svd_solver="arpack")
sc.pp.neighbors(adata, n_pcs=40)
sc.tl.umap(adata)
sc.pl.umap(adata, color=["batch", "leiden"])

方案：使用 PAGA + Diffusion Pseudotime 推断轨迹

适用场景：细胞呈现发育连续谱，而非离散的聚类。

# PAGA graph abstraction
sc.tl.paga(adata, groups="leiden")
sc.pl.paga(adata, color="leiden", threshold=0.03)

# Reinitialize UMAP with PAGA layout
sc.tl.umap(adata, init_pos="paga")

# Diffusion pseudotime — set root cell
adata.uns["iroot"] = np.flatnonzero(adata.obs["leiden"] == "0")[0]
sc.tl.diffmap(adata)
sc.tl.dpt(adata)
sc.pl.umap(adata, color=["dpt_pseudotime", "leiden"])

方案：绘制适用于发表的图表

适用场景：为论文或演示文稿生成图表。

sc.settings.set_figure_params(dpi=300, frameon=False, figsize=(5, 5))

# UMAP with custom palette
sc.pl.umap(
    adata, color="cell_type",
    palette="Set2",
    legend_loc="on data",
    legend_fontsize=10,
    legend_fontoutline=2,
    title="",
    save="_celltype_publication.pdf",
)

# Stacked violin plot of top markers
flat_markers = [g for genes in marker_genes.values() for g in genes[:2]]
sc.pl.stacked_violin(
    adata, var_names=flat_markers, groupby="cell_type",
    use_raw=True, swap_axes=True,
    save="_markers_violin.pdf",
)

方案：比较不同条件下的差异表达

适用场景：在特定细胞类型中比较处理组与对照组。

# Subset to cell type of interest
adata_t = adata[adata.obs["cell_type"] == "CD4 T cells"].copy()

# DE between conditions
sc.tl.rank_genes_groups(adata_t, groupby="condition", groups=["treated"], reference="control", method="wilcoxon", use_raw=True)
sc.pl.rank_genes_groups_volcano(adata_t, groups=["treated"])

de_results = sc.get.rank_genes_groups_df(adata_t, group="treated")
sig_genes = de_results[de_results["pvals_adj"] < 0.05]
print(f"Significant DE genes: {len(sig_genes)}")

预期输出

• results/annotated_data.h5ad — 完整的 AnnData，包含嵌入（PCA、UMAP）、聚类标签、细胞类型注释和原始计数层
• results/cell_metadata.csv — 按细胞排列的表格：条形码、聚类、细胞类型、QC 指标（n_genes、total_counts、pct_mt）
• results/marker_genes.csv — 每个聚类的 DE 基因：基因名称、对数倍数变化、校正后的 p 值
• QC 图表：小提琴图（n_genes、total_counts、pct_mt）、散点图
• 嵌入图表：按聚类、细胞类型、QC 指标、伪时间着色的 UMAP
• 标记基因图表：点图、热图、经典标记基因的堆叠小提琴图

故障排除

问题	原因	解决方案
ModuleNotFoundError: leidenalg	未安装 Leiden 包	pip install leidenalg igraph
PCA 期间出现 MemoryError	稠密矩阵超出 RAM 容量	使用 sc.pp.pca(adata, chunked=True, chunk_size=20000)
UMAP 显示单个团块	过滤过度或 HVG 数量过少	放宽 QC 阈值；将 n_top_genes 增加到 3000-5000
聚类过多或过少	分辨率不匹配	将 resolution 从 0.1 调整到 2.0，并比较 UMAP 图
找不到 KeyError: 'MT-'	物种前缀错误	小鼠使用 mt-（小写），人类使用 MT-
批次效应主导 UMAP	未校正的技术变异	应用 ComBat 方案，或使用 Harmony/scVI 进行整合
出现 adata.raw is None 错误	子集筛选前忘记保存 raw	归一化后、筛选 HVG 子集前，设置 adata.raw = adata
标记基因缺乏特异性	分辨率过低，导致细胞类型合并	提高分辨率，或使用 logistic_regression 方法
regress_out 运行缓慢	细胞数量较大	跳过 regress_out；仅使用 sc.pp.scale()（通常已足够）
rank_genes_groups 中出现 ValueError	某个聚类中的细胞过少	在进行 DE 前移除细胞数少于 10 的聚类：adata = adata[adata.obs.groupby('leiden').filter(lambda x: len(x)>=10).index]

随附资源

此技能包含用于深入查询的参考文件。当主 SKILL.md 需要更多详情时，按需读取这些文件。

references/api_reference.md

按模块（sc.pp.、sc.tl.、sc.pl.*）整理的所有 scanpy 函数速查表，包含完整签名、常用参数和 AnnData 结构参考。适用场景：需要查找特定操作的确切函数名称或参数时。

references/plotting_guide.md

全面的可视化指南：质控图、UMAP 样式设置、标记基因图（点图、热图、堆叠小提琴图）、轨迹图、多面板图、出版级自定义设置和配色方案建议。适用场景：为出版物或演示文稿制作图表时。

references/standard_workflow.md

针对各个流程阶段的详细分步参考，包含决策点、参数选择依据和替代方法（基于 MAD 的过滤、scran 归一化、自动注释）。适用场景：从头开始执行完整分析，且需要比上方“工作流”部分更深入的指导时。

参考资料

• Scanpy 文档 — 官方 API 参考和教程（BSD-3）
• Scanpy PBMC3k 教程 — 经典演示（BSD-3）
• Galaxy Training：使用 Scanpy 对 3K PBMC 进行聚类 — 分步教程（CC-BY）
• Luecken & Theis（2019） — “单细胞 RNA-seq 分析的当前最佳实践：教程”，Mol Syst Biol
• Heumos et al.（2023） — “跨模态单细胞分析的最佳实践”，Nat Rev Genet
• scverse 生态系统 — 相关工具：anndata、scvi-tools、squidpy、cellrank