实验 4：单细胞 RNA-seq 入门（PBMC 3k · scanpy）

预计用时：8–12 小时。
对应学习指南「项目三：单细胞 RNA-seq 分析」。批量 RNA-seq 测的是一管细胞的平均值；单细胞测序则给每个细胞
各测一份表达谱，所以能回答"这管血里有哪几种细胞"这类问题。
数据和流程都是领域内最经典的入门案例，结果可以和 scanpy、Seurat 的官方教程逐项对照。

一、问题与数据

问题：一份健康人的外周血单个核细胞（PBMC）里，有哪些免疫细胞类型？各占多少？每种细胞靠哪些基因认出来？

数据：10x Genomics 公开的 PBMC 3k 数据集，约 2700 个细胞 × 32738 个基因。

course fetch 4 下载约 7.6 MB 的压缩包，并自动解压到 data/pbmc3k/
矩阵非常稀疏：每个细胞通常只测到 1000 个左右的基因，其余全是 0。所以要用稀疏矩阵（scipy.sparse）存储

二、目标

理解单细胞数据的特点：稀疏、细胞之间测序深度差异大、会混进空液滴、双细胞和死细胞
手写质控指标和标准化，并和 scanpy 的结果对上
走通标准流程：高变基因 → PCA → 邻居图 → Leiden 聚类 → UMAP
用标记基因把"簇 0、簇 1……"翻译成细胞类型，并理解为什么要先做 z-score
读懂 UMAP 图，知道它能说明什么、不能说明什么

三、前置

概念	一句话	在代码里
10x 液滴测序	每个细胞被包进一个带条形码的油滴里，测序后按条形码把读段分回各个细胞	`barcodes.tsv`
UMI	分子标签，能去掉 PCR 扩增造成的重复，计数更准	矩阵里的整数
AnnData	scanpy 的数据容器：`.X` 是表达矩阵（细胞×基因），`.obs` 是细胞表，`.var` 是基因表，`.obsm` 存降维结果	`adata`
空液滴 / 碎片	检测到的基因太少（< 200 个）	`n_genes` 太低
双细胞（doublet）	两个细胞进了同一个液滴，检测到的基因会异常多	`n_genes` 太高
濒死细胞	细胞膜破了，胞质里的 mRNA 漏掉，线粒体 RNA 的比例就会偏高	`pct_mt` 太高
高变基因	在细胞之间变化最大的基因，最能区分细胞类型；只用它们可以降噪、提速	`highly_variable`
邻居图 + Leiden	在 PCA 空间里给每个细胞连上它最近的 10 个邻居，再在这张图上找"抱团"的社区	`leiden` 列
UMAP	把邻居图摊平到二维，用来可视化。两团之间的距离和团的大小都没有定量含义	`X_umap`
标记基因	只在某一类细胞里高表达的基因，例如 B 细胞表达 MS4A1（也叫 CD20）	`MARKERS`

PBMC 里常见的 8 类细胞：CD4 T、CD8 T、NK（自然杀伤细胞）、B、CD14 单核细胞、FCGR3A 单核细胞、树突状细胞（DC）、血小板。

四、任务

先运行 course fetch 4 和 course start 4，然后在 交付/pbmc.py 里完成。

任务 1：质控指标 `qc_metrics`

在稀疏矩阵上手算：每个细胞的基因数、总计数、线粒体比例。
测试会拿 sc.pp.calculate_qc_metrics 对答案。总计数为 0 的细胞不能出现除以 0 的错误。

任务 2：过滤 `qc_filter`

先过滤基因（至少在 3 个细胞里出现），再算质控指标，然后过滤细胞：基因数在 200 到 2500 之间，线粒体比例低于 5%。
验收：剩 2638 个细胞 × 13714 个基因，和 scanpy、Seurat 官方教程完全一致。

任务 3：标准化 `normalize_log`

每个细胞缩放到总计数 1 万，再取 log1p。在稀疏矩阵上完成，不能把矩阵变成稠密的，也不能改动输入。
测试会拿 sc.pp.normalize_total 加 sc.pp.log1p 的结果对答案。

任务 4：降维与聚类 `cluster`

按 docstring 的顺序调用 scanpy。重点是理解每一步的输入和输出，以及 adata.raw 为什么要在挑高变基因之前保存。
验收：1838 个高变基因，7 到 12 个簇。

任务 5：注释 `cluster_marker_means` / `annotate`

标记基因的表达量级相差很大：LYZ 在所有细胞里都很高，FCER1A 普遍很低。所以要先在各簇之间做 z-score，再打分。
验收：8 种细胞类型全部找到，各类数量落在文献范围内。例如 B 细胞 320–360 个，血小板 9–16 个。

任务 6：`course run 4`

交付/ 里会生成：qc_violin.png、umap_clusters.png、umap_celltypes.png、markers_dotplot.png、marker_genes_top10.csv、cell_type_counts.csv、结果摘要.md。

思考题（写进 `交付/思考题.md`）

把 max_pct_mt 从 5 改成 20，会多出多少个细胞？这些细胞在 UMAP 上落在哪里？它们可能是什么？
把 resolution 改成 0.5 和 2.0 各跑一次，簇的数量怎么变？CD4 T 细胞会不会被拆成好几个簇？"簇"和"细胞类型"是一回事吗？
UMAP 上 B 细胞离 T 细胞很远，离单核细胞近一点。能不能据此说"B 细胞和单核细胞更相似"？为什么？
这份数据只有一个人。如果要比较"病人和健康人的 CD14 单核细胞"，实验设计上要注意什么？提示：把实验 2 的"样本"和这里的"细胞"区分开，再查一下 pseudobulk。

选做

用 sc.tl.rank_genes_groups 找每个簇的前 5 个标记基因，和 MARKERS 对照，看有没有意外的发现。
用 scrublet（sc.pp.scrublet）做双细胞检测，比较它和"基因数 < 2500"这条粗过滤的差别。

五、验收标准

检查点	通过标准	对应测试
质控	和 scanpy 完全一致；空细胞不报错	`任务1_*`
过滤	2638 × 13714	`任务2_*`
标准化	和 scanpy 误差 < 1e-5；稀疏；不改动输入	`任务3_*`
聚类	1838 个高变基因；7–12 个簇；有 UMAP	`任务4_*`
注释	z-score 例子正确；8 类细胞齐全；数量在范围内；B 细胞的标记基因含 CD79A 或 MS4A1	`任务5_*`
交付	`course check 4` 全部通过；7 个文件齐全；思考题已写	`任务6_*`

六、常见错误

现象	原因	解决
内存暴涨、运行很慢	把稀疏矩阵 `.toarray()` 成了稠密矩阵	全程用稀疏矩阵运算
`X.sum(axis=1)` 的形状是 (n, 1)	稀疏矩阵求和返回的是矩阵	`np.asarray(...).ravel()`
过滤后细胞数不是 2638	先过滤细胞、后过滤基因，或者边界条件写成了 `<=`	按 docstring 的顺序和不等号来
过滤后基因数是 32738	没有过滤"少于 3 个细胞"的基因	先做基因过滤
找标记基因时说某个基因不存在	那个基因不是高变基因，已经被切掉了	用 `adata.raw`，在切片之前保存
DC 被标成 CD14 Mono	没做 z-score，被 LYZ 这类高表达基因带偏了	先按列做 z-score
每次运行簇的编号都不同	没固定随机种子	统一传 `random_state=seed`

参考实现：pbmc

卡住超过 20 分钟再看。对照着看自己是哪一步想岔了 —— 抄一遍没用。

展开参考解

看答案

"""实验 4 · 单细胞 RNA-seq 入门（参考解）

数据：10x Genomics PBMC 3k——一位健康供者的外周血单个核细胞，约 2700 个细胞。
问题：血液里有哪些细胞类型？各占多少？每种细胞靠哪些基因认出来？
流程参照 scanpy 官方教程（Wolf et al. 2018）和 Seurat 的 PBMC3k 教程。
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
import pandas as pd
import scipy.sparse as sp

HERE = Path(__file__).resolve().parent
ROOT = HERE.parents[2]
sys.path.insert(0, str(ROOT / "tools"))
import datasets  # noqa: E402

# 经典 PBMC 标记基因（Seurat PBMC3k 教程）。每个细胞类型挑 2–4 个
MARKERS: dict[str, list[str]] = {
    "CD4 T": ["IL7R", "CD3E", "CD3D", "LDHB"],
    "CD8 T": ["CD8A", "CD8B", "CD3E", "CCL5"],
    "NK": ["GNLY", "NKG7", "KLRD1", "PRF1"],
    "B": ["MS4A1", "CD79A", "CD79B"],
    "CD14 Mono": ["CD14", "LYZ", "S100A8"],
    "FCGR3A Mono": ["FCGR3A", "MS4A7", "LST1"],
    "DC": ["FCER1A", "CLEC10A", "CST3"],
    "Platelet": ["PPBP", "PF4"],
}

PARAMS = dict(min_genes=200, min_cells=3, max_genes=2500, max_pct_mt=5.0,
              n_pcs=40, n_neighbors=10, resolution=1.0, seed=0)


def load_pbmc3k():
    import scanpy as sc

    adata = sc.read_10x_mtx(datasets.pbmc3k_matrix_dir(), var_names="gene_symbols", cache=False)
    adata.var_names_make_unique()
    return adata


# ---------------------------------------------------------------- 任务 1：质控指标
def qc_metrics(X, var_names) -> pd.DataFrame:
    """X：细胞×基因 的原始计数（scipy 稀疏矩阵）。返回每个细胞的 n_genes / total_counts / pct_mt。"""
    X = sp.csr_matrix(X)
    total = np.asarray(X.sum(axis=1)).ravel()
    n_genes = np.asarray((X > 0).sum(axis=1)).ravel()
    mt = np.asarray([str(g).startswith("MT-") for g in var_names])
    mt_counts = np.asarray(X[:, mt].sum(axis=1)).ravel()
    pct = np.divide(mt_counts, total, out=np.zeros_like(total, dtype=float), where=total > 0) * 100
    return pd.DataFrame({"n_genes": n_genes, "total_counts": total, "pct_mt": pct})


# ---------------------------------------------------------------- 任务 2：过滤
def qc_filter(adata, min_genes=200, min_cells=3, max_genes=2500, max_pct_mt=5.0):
    """先去掉几乎空的液滴和几乎不表达的基因，再去掉疑似双细胞（基因太多）和濒死细胞（线粒体比例高）。"""
    counts_per_gene = np.asarray((sp.csr_matrix(adata.X) > 0).sum(axis=0)).ravel()
    adata = adata[:, counts_per_gene >= min_cells].copy()
    qc = qc_metrics(adata.X, adata.var_names)
    qc.index = adata.obs_names
    adata.obs[["n_genes", "total_counts", "pct_mt"]] = qc
    keep = (qc["n_genes"] >= min_genes) & (qc["n_genes"] < max_genes) & (qc["pct_mt"] < max_pct_mt)
    return adata[keep.to_numpy()].copy()


# ---------------------------------------------------------------- 任务 3：标准化
def normalize_log(X, target_sum: float = 1e4):
    """每个细胞缩放到总数 target_sum，再 log1p。输入输出都是 CSR 稀疏矩阵（不改原矩阵）。"""
    X = sp.csr_matrix(X, dtype=np.float64, copy=True)
    totals = np.asarray(X.sum(axis=1)).ravel()
    scale = np.divide(target_sum, totals, out=np.zeros_like(totals), where=totals > 0)
    X = sp.diags(scale) @ X
    X.data = np.log1p(X.data)
    return sp.csr_matrix(X)


# ---------------------------------------------------------------- 任务 4：降维与聚类
def cluster(adata, n_pcs=40, n_neighbors=10, resolution=1.0, seed=0):
    """adata.X 已是 log 标准化表达。挑高变基因 → 缩放 → PCA → 邻居图 → Leiden → UMAP。
    返回只含高变基因的新 AnnData（raw 里保留全部基因的 log 表达，用于找标记基因和画图）。"""
    import scanpy as sc

    sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
    adata.raw = adata
    hv = adata[:, adata.var["highly_variable"]].copy()
    sc.pp.scale(hv, max_value=10)
    sc.tl.pca(hv, n_comps=50, svd_solver="arpack", random_state=seed)
    sc.pp.neighbors(hv, n_neighbors=n_neighbors, n_pcs=n_pcs, random_state=seed)
    sc.tl.leiden(hv, resolution=resolution, random_state=seed, flavor="igraph", n_iterations=2, directed=False)
    sc.tl.umap(hv, random_state=seed)
    return hv


# ---------------------------------------------------------------- 任务 5：注释
def cluster_marker_means(adata, cluster_key="leiden") -> pd.DataFrame:
    """簇 × 标记基因 的平均 log 表达（用 raw 里的全部基因）。数据里没有的标记基因跳过。"""
    genes = [g for g in dict.fromkeys(g for gs in MARKERS.values() for g in gs) if g in adata.raw.var_names]
    X = adata.raw[:, genes].X
    X = X.toarray() if sp.issparse(X) else np.asarray(X)
    return pd.DataFrame(X, columns=genes, index=adata.obs_names).groupby(adata.obs[cluster_key], observed=True).mean()


def annotate(means: pd.DataFrame, markers: dict[str, list[str]] = MARKERS) -> pd.Series:
    """每个标记基因先在各簇之间做 z-score（让高表达基因和低表达基因可比），
    每个细胞类型的得分 = 它的标记基因 z-score 的平均；每个簇取得分最高的类型。"""
    z = (means - means.mean(axis=0)) / means.std(axis=0, ddof=0).replace(0, np.nan)
    z = z.fillna(0)
    scores = pd.DataFrame({t: z[[g for g in gs if g in z.columns]].mean(axis=1) for t, gs in markers.items()})
    return scores.idxmax(axis=1).rename("cell_type")


# ---------------------------------------------------------------- 流程与交付
def run_pipeline(params=PARAMS):
    raw = load_pbmc3k()
    adata = qc_filter(raw, params["min_genes"], params["min_cells"], params["max_genes"], params["max_pct_mt"])
    adata.layers["counts"] = adata.X.copy()
    adata.X = normalize_log(adata.X)
    hv = cluster(adata, params["n_pcs"], params["n_neighbors"], params["resolution"], params["seed"])
    labels = annotate(cluster_marker_means(hv))
    hv.obs["cell_type"] = hv.obs["leiden"].map(labels).astype("category")
    return raw, hv, labels


def main(out_dir) -> None:
    import warnings

    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    import scanpy as sc

    warnings.filterwarnings("ignore")
    out = Path(out_dir)
    out.mkdir(parents=True, exist_ok=True)
    sc.settings.figdir = out

    raw, hv, labels = run_pipeline()
    print(f"细胞：{raw.n_obs} → 质控后 {hv.n_obs}；基因：{raw.n_vars} → 高变基因 {hv.n_vars}")

    # 质控小提琴图（过滤后）
    fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
    for ax, col in zip(axes, ["n_genes", "total_counts", "pct_mt"]):
        ax.violinplot(hv.obs[col], showmedians=True)
        ax.set_title(col)
        ax.set_xticks([])
    fig.tight_layout()
    fig.savefig(out / "qc_violin.png", dpi=150)
    plt.close(fig)

    for key, name in [("leiden", "umap_clusters.png"), ("cell_type", "umap_celltypes.png")]:
        fig = sc.pl.umap(hv, color=key, legend_loc="on data", legend_fontsize=7, frameon=False, show=False, return_fig=True)
        fig.savefig(out / name, dpi=150, bbox_inches="tight")
        plt.close(fig)

    sc.tl.rank_genes_groups(hv, "cell_type", method="wilcoxon", use_raw=True)
    top = pd.DataFrame(hv.uns["rank_genes_groups"]["names"]).head(10)
    top.to_csv(out / "marker_genes_top10.csv", index=False, encoding="utf-8-sig")

    genes = [g for g in dict.fromkeys(g for gs in MARKERS.values() for g in gs[:2]) if g in hv.raw.var_names]
    dp = sc.pl.dotplot(hv, genes, groupby="cell_type", use_raw=True, show=False, return_fig=True)
    dp.savefig(out / "markers_dotplot.png", dpi=150, bbox_inches="tight")
    plt.close("all")

    comp = hv.obs["cell_type"].value_counts().rename("n_cells").to_frame()
    comp["percent"] = (comp["n_cells"] / comp["n_cells"].sum() * 100).round(1)
    comp.to_csv(out / "cell_type_counts.csv", encoding="utf-8-sig")

    lines = ["# 实验 4 结果摘要（PBMC 3k）\n",
             f"- 细胞 {raw.n_obs} → 质控后 {hv.n_obs}；高变基因 {hv.n_vars}；Leiden 簇 {hv.obs['leiden'].nunique()} 个\n",
             "## 簇 → 细胞类型\n", "| cluster | cell type | n |", "|---|---|---:|"]
    sizes = hv.obs["leiden"].value_counts()
    lines += [f"| {c} | {t} | {sizes[c]} |" for c, t in labels.sort_index(key=lambda s: s.astype(int)).items()]
    lines += ["", "## 细胞类型组成\n", "| cell type | n | % |", "|---|---:|---:|"]
    lines += [f"| {t} | {int(n)} | {p} |" for t, n, p in zip(comp.index, comp["n_cells"], comp["percent"])]
    (out / "结果摘要.md").write_text("\n".join(lines) + "\n", encoding="utf-8")
    print("\n".join(lines))

← 实验 3：富集分析（ORA 超几何检验 · BH 校正 · GSEA）目录

实验 4：单细胞 RNA-seq 入门（PBMC 3k · scanpy）

一、问题与数据

二、目标

三、前置

四、任务

任务 1：质控指标 qc_metrics

任务 2：过滤 qc_filter

任务 3：标准化 normalize_log

任务 4：降维与聚类 cluster

任务 5：注释 cluster_marker_means / annotate

任务 6：course run 4

思考题（写进 交付/思考题.md）

选做