📊 scRNA-seq 标准管道总览

把整个分析当成一条 ETL 数据管道：输入一个又大又稀疏的矩阵（细胞 × 基因，整数计数，90%+ 是 0），经过一串变换，输出「每个细胞是什么类型」。

QC过滤 → 归一化 → 选高变基因 → 降维(PCA→UMAP) → 聚类(Leiden) → 注释 → 差异表达

数据装在 AnnData 对象里，整条管道就是不断往它身上加东西：

.X —— 表达矩阵（行=细胞，列=基因）
.obs —— 每个细胞的元数据（行注释：QC 指标、聚类标签、细胞类型…）
.var —— 每个基因的元数据（列注释：是否线粒体、是否高变…）
.uns —— 非结构化的全局结果（邻居图参数、配色…）

记住这个结构：后面每一步在改谁（加 obs 列？裁 var？换 X？），就一目了然。

import scanpy as sc
import numpy as np
sc.settings.verbosity = 1
sc.settings.set_figure_params(dpi=80)

第 1 步：加载数据 + 认识 AnnData

干什么：下载 PBMC 3k（外周血单个核细胞，~2700 个，Scanpy 官方入门标配），并把基因名去重。

看什么：输出 AnnData object with n_obs × n_vars = 2700 × 32738

n_obs = 细胞数（行）
n_vars = 基因数（列）

动手感受数据结构（把注释逐行打开跑一下）：

adata.X[:5,:5].toarray() —— 稀疏矩阵的一角，绝大多数是 0（这叫 dropout：不是真不表达，是没测到）
adata.obs.head() —— 此刻几乎是空的，后面每一步会往这里加列
adata.var.head() —— 每个基因的信息（现在只有 gene_ids）

对比 bulk RNA-seq：以前是把上百万细胞混在一起测平均，丢了异质性。单细胞把这个平均「拆」到每个细胞——这才是能区分组织里有哪些细胞类型的基础。

adata = sc.datasets.pbmc3k()      # 自动下载
adata.var_names_make_unique()     # 基因名去重，养成习惯
print(adata)
# print(adata.X[:5,:5].toarray())  # 稀疏矩阵的一角，大多是 0
# print(adata.obs.head())   # 每细胞元数据
print(adata.var.head())  # 每基因元数据

第 2 步：质控 QC —— 扔掉「不是真细胞」的行

液滴测序里混着空泡、濒死细胞、还有 doublet（两个细胞挤进一个液滴被当成一个）。这些是技术噪声，不清掉会污染后面所有步骤。

先算三个每细胞指标（下面这个 cell）：

n_genes_by_counts —— 这个细胞测到多少种基因。太少 = 空液滴/破细胞
total_counts —— 总分子数，即测序深度
pct_counts_mt —— 线粒体基因占比。太高 = 细胞在死（膜破了，胞质 mRNA 漏走，线粒体 mRNA 还在，占比飙高）

adata.var['mt'] = adata.var_names.str.startswith('MT-') 是在标记哪些基因是线粒体基因（人类基因名以 MT- 开头），有了这个标记才能算出 mt 占比。

# 给每个基因打一个布尔标签，名字以MT-开头 -> True(线粒体基因)，其余的是False
# 先告诉程序"哪些基因是线粒体基因"，再让 scanpy 一次性算出所有 QC 指标，为后续过滤低质量细胞做准备。
adata.var['mt'] = adata.var_names.str.startswith('MT-')   # 线粒体基因（人类是 MT- 开头）
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

看分布（先看再砍）—— 怎么读这几张图

先搞懂：小提琴图（violin）在画什么

把每张小提琴图想象成一个竖着的、对称的「烛台/陀螺」：

纵轴 = 指标数值；横轴没含义。
某高度上「胖」还是「瘦」= 有多少细胞落在那个数值。肚子最胖处 = 大多数细胞扎堆的值；越往上越细 = 越少细胞有那么高的值。本质就是直方图竖过来再左右镜像（两边一样是镜像，不是两组数据）。
黑点 = 一个个细胞；它们左右散开（jitter）只为不重叠，横向位置随机、没意义——只看每个点的高度。

一句话：胖肚子 = 正常细胞扎堆处；向上拖出的细尾巴 + 稀疏黑点 = 数值异常的细胞，QC 要砍的就是这些尾巴。

逐张看

图1 n_genes_by_counts（每细胞测到多少种基因）：胖肚子 ~500-1200（正常）；少数点拖到 2500-3500 → 可疑 doublet（两个细胞挤一起，基因种类翻倍）。→ 下一步 < 2500 剪掉这条尾巴。
图2 total_counts（每细胞总分子数 = 测序深度）：胖肚子 ~1500-3500，尾巴到 16000。同理，特别高的也常是 doublet（和图1联动）。
图3 pct_counts_mt（线粒体占比 %）：胖肚子 ~1-3%（健康）；点拖到 5/10/20% → 濒死细胞（膜破，胞质 mRNA 漏走，线粒体的还在，占比飙高）。→ 下一步 < 5 砍掉死细胞。

两张散点图（每点仍是一个细胞，看「两指标的关系」）

total_counts(x) vs pct_counts_mt(y)：看上方点（线粒体高）。典型死细胞 = 总数低 + 线粒体高。确认 5% 这条线合不合理。
total_counts(x) vs n_genes_by_counts(y)：正常应正相关（点排成往右上的曲线带）；偏离主带的点 = 可疑 doublet。确认数据行为正常、揪离群点。

看图的目的 = 用眼睛定阈值

看哪张	看到什么	定哪个阈值
图1 n_genes	主群到 ~2000，再上是稀疏尾巴	`n_genes_by_counts < 2500`
图3 pct_mt	主群在 3% 以下，5% 以上零星	`pct_counts_mt < 5`

「主群在哪结束、垃圾尾巴从哪开始」这条线是拿眼睛估的，没有标准答案——这正是清单 3.0 的 L2 判断点。将来把这步做成 Agent 的 tool，这两个数就是要暴露给用户调、或写护栏（「砍掉超过 X% 细胞就报警」）的参数。

# 看图定阈值
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')

过滤看什么：细胞数从 2700 降到 ~2638，基因数降到 ~13714。 ③ 坑/决策点：2500 和 5% 这两个阈值不是金科玉律，是看着上面的图定的。砍太狠会丢稀有真细胞，太松留垃圾。← 这就是 3.0 说的 L2 判断；将来这就是你 tool 里要暴露/加护栏的参数。

sc.pp.filter_cells(adata, min_genes=200)      # 细胞至少测到200种基因
sc.pp.filter_genes(adata, min_cells=3)        # 基因至少在3个细胞里出现
adata = adata[adata.obs.n_genes_by_counts < 2500, :]   # 砍掉疑似doublet
adata = adata[adata.obs.pct_counts_mt < 5, :]          # 砍掉濒死细胞
adata

第 3 步：归一化 —— 让细胞之间可比

问题：每个细胞捕获的总分子数天差地别（纯技术原因）。细胞 A 测到 5000 条、B 测到 500 条，直接比没意义。

两步：

normalize_total(target_sum=1e4) —— 每个细胞缩放到相同总量（1 万）
log1p —— 取 log(1+x)。基因表达跨好几个数量级，log 压缩量纲、让方差更均匀、防止高表达基因主导。等价于你做 ML 前的特征缩放。

adata.raw = adata —— 存一份当前（全基因、已归一化）状态的备份。第 4 步会把基因裁到 ~2000 个，但第 9 步画 marker 图想用全部基因，靠这份备份。别漏。

⚠️ 你这次跑出了 WARNING: adata.X seems to be already log-transformed。这是 scanpy 的保守提示——通常是这个 cell 被重复执行了（log 又套了一层）。干净做法：从第 1 步重新跑一遍（菜单 Kernel → Restart Kernel and Run All Cells），保证每个预处理步骤只执行一次。本次教学不影响，但要养成「预处理步骤不重复跑」的习惯。

sc.pp.normalize_total(adata, target_sum=1e4)   # 每细胞归一到总量1万
sc.pp.log1p(adata)                             # log(1+x)
adata.raw = adata                              # 存一份当前状态（注释画图时用全基因）

第 4 步：选高变基因 HVG —— 只留有区分力的列

~2 万基因里，大部分是「管家基因」（每个细胞都差不多，对区分类型没用）。挑出变化最大的 ~2000 个。

为什么：降维 + 去噪，只盯真正能区分细胞的基因。就是特征选择。

看什么：sc.pl.highly_variable_genes 的图里，黑点 = 高变（保留）、灰点 = 不变（丢弃）。最后裁完 n_vars 从 ~13714 降到 ~1800-2000。

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata)
adata = adata[:, adata.var.highly_variable]    # 只留高变基因
adata

第 5 步：缩放 + PCA —— 降到几十维

regress_out（可选、较慢）—— 回归掉 total_counts、线粒体比例这些技术变量的残留影响。想快可以先跳过。
scale(max_value=10) —— 每个基因做 z-score 标准化，截断到 10（防极端值）。之后 .X 变成有正有负的标准化值。
pca —— 把 ~2000 维压到几十个主成分，抓主要变异轴、去噪。

看什么：pca_variance_ratio 图找「肘部」（曲线拐弯变平处），决定下一步 n_pcs 用多少。PBMC3k 大概 30-40 个就够。

sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])  # 可选：回归掉技术变量
sc.pp.scale(adata, max_value=10)                             # z-score，截断到10
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True)

第 6 步：邻居图 + UMAP

neighbors(n_neighbors=10, n_pcs=40) —— 在 PCA 空间建 k 近邻图。这张图是后面聚类和 UMAP 共同的基础。 n_pcs 用上一步看的肘部值。
umap —— 再压到 2 维只为了画图看。

看什么：一团点散成几坨——已经有结构了，但还没分组、没颜色。

⚠️ 铁律 caveat：UMAP 只能看、不能量。两坨在图上挨得近 ≠ 生物学相近，UMAP 会扭曲全局距离。任何定量结论都不能基于 UMAP 的二维坐标。

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)   # n_pcs 用上一步看的肘部值
sc.tl.umap(adata)
sc.pl.umap(adata)

第 7 步：聚类 Leiden —— 把相似细胞分组

在上一步那张近邻图上做社区发现（Leiden 算法），把抱团的细胞分成一簇簇 = 候选细胞类型/状态。

最重要的旋钮 resolution：调高 = 簇更多更细，调低 = 簇更少更粗。PBMC3k 在 1.0 附近大约 8 个簇（但随 scanpy/leiden 版本会有出入，第 9 步贴标签时要按你实际的簇数来）。

动手体会：把 resolution 改成 0.5 和 2.0 各跑一遍，直观看到簇变少/变多。这就是「欠聚类 / 过度聚类」——没有标准答案，要靠下一步 marker 验证。

小提示：新版 scanpy 可能提示 flavor 的 FutureWarning，不影响。想消掉写 sc.tl.leiden(adata, resolution=1.0, flavor="igraph", n_iterations=2, directed=False)。

sc.tl.leiden(adata, resolution=2.0)
sc.pl.umap(adata, color=['leiden'])

第 8 步：找每簇的特征基因 —— 为贴标签做准备

rank_genes_groups 找出每个簇相比其他簇显著高表达的基因（Wilcoxon 检验）。这些就是给簇「贴标签」的线索。

看什么：每个簇排名最高的一串基因名。

⚠️ caveat：scRNA-seq 的 DE p 值普遍偏乐观（用同一份数据先聚类、又在同一份数据上检验，叫 double dipping）。当线索用，别当统计铁证。

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

第 9 步：注释 —— 给每个簇贴上「这是什么细胞」

看每个簇高表达哪些已知 marker 基因，对照贴标签。用 dotplot 看 marker 在各簇的表达（点越大 = 表达该基因的细胞比例越高，越红 = 平均表达越强）：

marker	细胞类型
IL7R	CD4 T 细胞
CD14, LYZ	CD14+ 单核细胞
MS4A1	B 细胞
CD8A	CD8 T 细胞
GNLY, NKG7	NK 细胞
FCGR3A, MS4A7	FCGR3A+ 单核细胞
FCER1A, CST3	树突状细胞 (DC)
PPBP	巨核细胞

这步最吃判断：marker 模糊、簇对不上、或某簇像 doublet 都正常——回到清单 3.0，「生物学上成不成立」该拉师兄师姐/导师把关，不归你。

marker_genes = ['IL7R','CD14','LYZ','MS4A1','CD8A','GNLY','NKG7','FCGR3A','MS4A7','FCER1A','CST3','PPBP']
sc.pl.dotplot(adata, marker_genes, groupby='leiden')

⚠️ 下面那个 cell 报错了：`new categories need to have the same number of items`

原因：new_names 写了 8 个名字，但你这次 Leiden 实际聚出的簇不是 8 个（簇数随 scanpy/leiden 版本和 resolution 变化）。盲抄别人的 8 个名字必然对不上——这正是清单 3.0 说的：不能照搬，要对应你自己的结果。

第一步，先看你到底有几个簇、各簇高表达什么基因：

print(adata.obs['leiden'].cat.categories)        # 你实际有几个簇、编号是什么
import pandas as pd
pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(10)   # 每簇 Top10 基因，对照第 9 步 marker 表

第二步，用「按簇编号映射」的稳健写法（不依赖簇的数量和顺序）：

# 按 dotplot + 上面 Top 基因的判断逐个填；键=簇编号字符串，值=细胞类型
cluster2type = {
    '0': 'CD4 T',
    '1': 'CD14+ Mono',
    # ... 把你实际的每一个簇都填上（有几个填几个）
}
adata.obs['cell_type'] = adata.obs['leiden'].map(cluster2type).astype('category')
sc.pl.umap(adata, color='cell_type', legend_loc='on data')

这样即使簇数是 7 或 9，也只是 cluster2type 多写少写几行，不会再因「数量对不上」报错。原来的 rename_categories(new_names) 要求名字数量 == 簇数量且按顺序，很脆；换成 map 字典更工程化——这也呼应你「领域是广度、工程是深度」的思路：用工程手段把脆弱的领域操作变稳。

new_names = ['CD4 T', 'CD14+ Mono', 'B', 'CD8 T', 'NK', 'FCGR3A+ Mono', 'DC', 'Megakaryocyte']
adata.rename_categories('leiden', new_names)
sc.pl.umap(adata, color='leiden', legend_loc='on data')

第 10 步：存盘

adata.write('pbmc3k_annotated.h5ad') 把带注释的结果存成 .h5ad。

这就是单细胞分析的标准产物。后面学空转（Squidpy + 10x Visium）时，数据格式完全一样是 .h5ad，80% 的步骤你已经会了——只多了「空间坐标」那一维。
✅ 跑到这里、能看到一张 UMAP 上每坨细胞标着 T/B/NK/单核/DC 等名字，你就完成了学习清单 3.4 检查点第 1 条：用 Scanpy 把 scRNA-seq 从 QC 跑到细胞类型注释。

adata.write('pbmc3k_annotated.h5ad')