单细胞测序数据分析

引言

单细胞 RNA-seq（scRNA-seq）把转录组分析的分辨率从「一群细胞」提升到「单个细胞」。传统的 bulk RNA-seq 测的是成千上万个细胞的平均值，掩盖了细胞异质性——而恰恰是这种异质性（肿瘤的亚克隆、免疫细胞的不同亚型、发育中的过渡态）往往是生物学问题的核心。单细胞的价值就在于揭示「平均」看不到的结构。

但单细胞分析的技术栈与 bulk 完全不同。它处理的对象是「数万个细胞 × 数万个基因」的超稀疏矩阵，99% 以上的值是 0（因为每个细胞捕获到的 mRNA 分子有限）。这带来了独特的挑战：如何区分「基因真的不表达」与「没测到」（dropout）？如何在极度稀疏的数据上做可靠的降维聚类？如何避免把技术噪声当成生物学异质性？

工程上的难点有三个：一是数据规模与稀疏性，几万个细胞的矩阵若按稠密存储会占用数十 GB 内存；二是分析步骤的主观性，从质控阈值到聚类分辨率，每一步都有大量「经验参数」，不同选择给出不同结果；三是可复现性差，UMAP 图每次运行都可能略有不同，细胞类型注释更是高度依赖专家判断。

本文按「技术路线 → 数据特征 → 上游处理 → 质控 → 归一化 → 降维 → 聚类注释 → 差异分析 → 轨迹与整合」的顺序展开。代码以 Scanpy（Python）为主、Seurat（R）为辅。读完你应该能独立完成一个单细胞分析，并理解每个参数背后的权衡。

目录

1. 单细胞测序的技术路线
2. 数据特征：稀疏矩阵与零膨胀
3. 上游处理：Cell Ranger 与计数矩阵
4. 质控：细胞与基因过滤
5. 归一化与高变基因选择
6. 降维：PCA 与 UMAP
7. 聚类与细胞类型注释
8. 标记基因与差异表达
9. 轨迹推断与数据整合

1. 单细胞测序的技术路线

单细胞测序的核心难题是「分离单个细胞」并「给每个细胞的转录本打上标签」。主流技术路线：

10x Chromium 是绝对主流：用微流控把单个细胞与一个凝胶珠（gel bead）包裹在油滴里。凝胶珠上有数百万条捕获探针，每条探针包含：

[测序引物][细胞条形码 Cell Barcode (16bp)][UMI (12bp)][poly-T 捕获 mRNA]

• 细胞条形码（Cell Barcode）：标记「这条读段来自哪个细胞」。同一个液滴里的所有 mRNA 共享同一 barcode。
• UMI（Unique Molecular Identifier）：标记「这条读段来自哪个原始 mRNA 分子」。用于去除 PCR 重复，准确计数。
• poly-T：捕获 mRNA 的 poly-A 尾，所以 10x 只测到 mRNA 的 3’ 端。

这个设计决定了数据的形式：每个「细胞 × 基因」的计数，是「该细胞中带特定 barcode 且特定基因的 UMI 数量」。理解 barcode 与 UMI 的作用，是理解单细胞计数矩阵的基础。

空间转录组是近年热点：它把「哪个细胞」的问题扩展为「在哪个位置」。10x Visium 用带位置条形码的载玻片捕获组织切片上的 mRNA，保留空间信息。分析时需要专门的空间方法（如 squidpy）。

2. 数据特征：稀疏矩阵与零膨胀

单细胞数据的第一个特征是极度稀疏。一个典型的 10x 数据集：5000 个细胞 × 20000 个基因 = 1 亿个矩阵元素，但非零元素可能只有 2000 万（20%），其余 80% 是 0。在更复杂的组织里，稀疏度可达 95% 以上。

import scanpy as sc
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
print(adata.X.shape)              # (5000, 20000)
print(adata.X.nnz / (5000*20000)) # 非零比例，通常 < 0.2

# 每个细胞的基因数分布
import numpy as np
genes_per_cell = np.asarray((adata.X > 0).sum(axis=1)).flatten()
print(f"中位基因数: {np.median(genes_per_cell):.0f}")
# 典型值：1000-3000 个基因/细胞

为什么是稀疏的？ 两个原因：一是每个细胞的 mRNA 捕获效率有限（通常只捕获到 10-20% 的转录本），大量低表达基因「没测到」；二是很多基因在特定细胞类型中确实不表达。前者是技术噪声（dropout），后者是生物学信号——区分两者是单细胞分析的核心难题。

零膨胀（zero inflation） 指数据中 0 的比例远超理论预期。这带来统计挑战：一个基因在细胞 A 中计数为 0，可能是「不表达」，也可能是「表达但没捕获到」。处理方式：

• 不要轻易把 0 当「不表达」：需要结合该基因在其他细胞中的表达判断；
• imputation（插补）：用统计方法填补 dropout 的 0（如 MAGIC、ALRA），但有争议——插补可能引入假信号；
• 稀疏矩阵存储：用 CSR/CSC 格式（scipy.sparse）存储，节省内存。参见 Python 内存与 GC 。

数据规模的工程影响：一个 10 万细胞的数据集，稠密矩阵需要 100000 × 20000 × 8 字节 = 16 GB，而稀疏存储只需几百 MB。所有单细胞工具都用稀疏矩阵，且很多步骤需要足够的内存（如降维、聚类）。

3. 上游处理：Cell Ranger 与计数矩阵

从 FASTQ 到计数矩阵由 Cell Ranger（10x 官方）或开源替代（STARsolo、alevin）完成：

# Cell Ranger 完整流程
cellranger count \
  --id=sample1 \
  --transcriptome=refdata-gex-GRCh38-2024-A \
  --fastqs=/data/fastq \
  --sample=sample1 \
  --expect-cells=5000 \
  --localcores=16 --localmem=64

上游处理的四个关键步骤：

1. barcode 校正：测序错误的 barcode（1 个碱基差异）会被校正到已知白名单中最接近的 barcode；
2. 比对：读段比对到参考转录组（不是基因组），因为只需要知道是哪个基因；
3. UMI 计数：同一个 barcode + 同一个基因 + 同一个 UMI 的读段合并为 1 个分子计数；
4. 细胞判定：区分「真细胞」与「空液滴」。空液滴里只有背景 mRNA，UMI 数很低。

细胞判定的核心是「barcode 排序曲线的拐点」。把 barcode 按总 UMI 数降序排列，真细胞在曲线顶部（高 UMI），空液滴在底部（低 UMI），拐点处是分界。Cell Ranger 的 --expect-cells 参数影响这个判定。

输出三个矩阵文件：

filtered_feature_bc_matrix/   ← 过滤后的真细胞（分析用这个）
raw_feature_bc_matrix/        ← 全部 barcode（含空液滴，调试用）
molecule_info.h5              ← 分子级信息（高级分析）

一个常见错误是「细胞数估计不准」：--expect-cells 设得太低会漏掉细胞，太高会混入空液滴。对于未知样本，可以先跑一次看 Cell Ranger 报告的「Estimated Number of Cells」，再调整。

4. 质控：细胞与基因过滤

单细胞质控是「过滤掉坏细胞和坏基因」。核心指标有三个：

import scanpy as sc
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()

# 标记线粒体基因（人类以 MT- 开头，小鼠以 mt- 开头）
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

# 查看分布，据此设阈值
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].hist(adata.obs["n_genes_by_counts"], bins=100); axes[0].set_title("genes/cell")
axes[1].hist(adata.obs["total_counts"], bins=100); axes[1].set_title("UMI/cell")
axes[2].hist(adata.obs["pct_counts_mt"], bins=100); axes[2].set_title("% mito")
plt.show()

线粒体比例是最重要的质控指标。活细胞的 mRNA 主要在胞质，线粒体基因占比低（人类通常 < 10-20%）。死细胞或濒死细胞的细胞膜破裂，胞质 mRNA 流失，只剩下被膜保护的线粒体 mRNA，导致线粒体比例飙升。所以高线粒体比例 = 死细胞，应过滤。

阈值设定经验（人类组织）：

# 典型过滤阈值（需按组织调整）
adata = adata[adata.obs["n_genes_by_counts"].between(200, 5000), :]
adata = adata[adata.obs["total_counts"] < 20000, :]
adata = adata[adata.obs["pct_counts_mt"] < 20, :]

双细胞（doublet） 是另一个常见问题：两个细胞被同一个液滴包裹，看起来像一个「表达两套标记基因」的异常细胞。检测工具：

# Scrublet 检测双细胞
import scrublet as scr
scrub = scr.Scrublet(adata.X)
doublet_scores, predicted_doublets = scrub.scrub_doublets()
adata.obs["doublet_score"] = doublet_scores
adata = adata[~predicted_doublets, :]

质控阈值的设定没有绝对标准，必须看分布而非套用固定值。不同组织（如心肌细胞的线粒体天然高、神经元 RNA 含量高）的合理范围差异很大。机械套用「线粒体 < 20%」会误杀某些细胞类型。

5. 归一化与高变基因选择

质控后要处理「测序深度差异」：有的细胞测到 5000 UMI，有的只有 1000，直接比较计数是不公平的。归一化方法：

# 1. 归一化到相同总计数（默认 1e4）
sc.pp.normalize_total(adata, target_sum=1e4)

# 2. 对数变换（压缩动态范围）
sc.pp.log1p(adata)

# 3. 保存原始计数（供后续差异分析）
adata.raw = adata

为什么是「归一化 + log」而不是 TPM？因为单细胞数据的稀疏性使 TPM 类方法不稳定（大量 0 导致除零或放大噪声）。normalize_total + log1p 是单细胞的标准做法：先归一化到相同总计数，再取对数压缩动态范围（避免高表达基因主导降维）。

高变基因（Highly Variable Genes, HVG） 选择是降维前的关键步骤：不是所有 20000 个基因都包含生物学信息，很多基因在所有细胞里表达量相近（无区分度）。只保留「在细胞间变异大」的基因：

# 选择高变基因
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=2000,      # 保留 2000 个
    flavor="seurat_v3",    # 用 seurat_v3 方法（对稀疏数据更稳健）
    subset=True
)
print(f"保留 {adata.n_vars} 个高变基因")

HVG 选择的价值：一是降维时聚焦有信息的基因，减少噪声；二是大幅降低计算量。2000 个 HVG 足以捕获主要细胞类型差异。

一个争议点：HVG 选择可能「偏向高表达基因」，导致稀有细胞类型的标记基因被排除。如果关注稀有细胞，可以适当增加 HVG 数量（如 3000-5000），或用「所有基因」做降维（代价是噪声增加）。

6. 降维：PCA 与 UMAP

单细胞降维分两步：PCA 做线性降维（2000 维 → 50 维），UMAP/t-SNE 做非线性可视化（50 维 → 2 维）。

# PCA：先标准化，再降维
sc.pp.scale(adata, max_value=10)       # 标准化（限制极端值）
sc.tl.pca(adata, n_comps=50, svd_solver="arpack")

# 确定用几个主成分（看方差解释率拐点）
sc.pl.pca_variance_ratio(adata, n_pcs=50, log=True)

# 近邻图（UMAP 与聚类的基础）
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)

# UMAP 可视化
sc.tl.umap(adata)
sc.pl.umap(adata, color=["leiden", "n_genes_by_counts"])

关键参数与含义：

PCA 的 n_pcs 选择：看方差解释率曲线（elbow plot），在拐点处取。太多主成分会引入噪声，太少会丢失信号。通常 20-50 之间。

UMAP 的解读陷阱：UMAP 图上「距离的远近」不等于「生物学相似度」。「两个簇靠得近」不代表它们相似，「簇的大小」也不代表细胞数量。UMAP 主要用于「看结构」（有几个簇、是否连续），而非精确的定量比较。这是初学者最容易误解的地方——不要从 UMAP 图直接读出「A 细胞比 B 细胞多 3 倍」这样的结论。

t-SNE 与 UMAP 的对比：t-SNE 保留局部结构好但全局结构差、慢；UMAP 兼顾局部与全局、快。当前 UMAP 是主流，但两者都应「只看局部关系，谨慎解读全局」。

7. 聚类与细胞类型注释

聚类把「转录组相似的细胞」分组。主流算法是 Leiden（Louvain 的改进版）：

# Leiden 聚类
sc.tl.leiden(adata, resolution=1.0, key_added="leiden")

# 分辨率影响簇的数量：分辨率高→簇多，低→簇少
for res in [0.3, 0.5, 1.0, 1.5]:
    sc.tl.leiden(adata, resolution=res, key_added=f"leiden_{res}")

# 可视化
sc.pl.umap(adata, color=["leiden", "leiden_0.3", "leiden_1.5"])

分辨率（resolution）是最重要的参数，它控制聚类的「粒度」：resolution 越大，簇越多、越细。没有「正确」的分辨率——取决于你要回答的问题。粗略的细胞类型（如 T 细胞、B 细胞）用低分辨率（0.3-0.5），细分的亚型（如 CD4+ 记忆 T、CD8+ 效应 T）用高分辨率（1.0-1.5）。建议「多分辨率聚类」，结合生物学知识选择。

细胞类型注释是单细胞分析中最需要生物学知识的步骤。三种方法：

方法一：标记基因手动注释。用已知的标记基因（marker）判断每个簇的身份：

# 已知标记基因（示例）
marker_genes = {
    "T cell": ["CD3D", "CD3E", "CD2"],
    "B cell": ["CD79A", "MS4A1", "CD19"],
    "Monocyte": ["CD14", "LYZ", "FCGR3A"],
    "NK cell": ["NKG7", "GNLY", "KLRD1"],
    "Endothelial": ["PECAM1", "VWF"],
}
sc.pl.dotplot(adata, marker_genes, groupby="leiden", standard_scale="var")

方法二：自动注释工具。用参考数据集自动匹配：

# CellTypist 自动注释
import celltypist
predictions = celltypist.annotate(adata, model="Immune_All_Low.pkl")
adata = predictions.to_adata()

方法三：与参考数据集比对。用 scmap、SingleR 等把细胞映射到已知注释的参考。

三种方法应结合使用：自动工具给初步建议，标记基因验证，专家知识定夺。自动注释工具的结果必须人工审核——它们依赖参考数据集的覆盖范围，对参考中没有的细胞类型会强行归类。

8. 标记基因与差异表达

聚类之后要找「每个簇的标记基因」，用于注释和后续验证：

# 找每个簇的标记基因
sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=10, sharey=False)

# 提取结果
import pandas as pd
result = sc.get.rank_genes_groups_df(adata, group="0")
print(result.head(20))

单细胞差异表达与 bulk 有本质区别，需要注意几个陷阱：

陷阱一：伪重复（pseudo-replication）。如果处理组有 3 个样本、对照组有 3 个样本，每个样本有几千个细胞，不能把细胞当独立样本——同一只小鼠的细胞不是独立的，它们共享个体效应。正确做法是把「样本」作为统计单位，用伪批量（pseudobulk）方法：

# 伪批量：把每个样本的细胞聚合成一个「bulk」样本，再做差异分析
import scanpy as sc
pseudobulk = sc.get.aggregate(adata, by="sample", func="sum")
# 然后用 DESeq2/edgeR 分析伪批量数据（见 RNA-seq 专题）

陷阱二：只按 p 值排序。单细胞细胞数多，p 值容易「显著」但不一定「有生物学意义」。应结合效应量（log fold change）与表达比例（在多少细胞中表达）。

陷阱三：忽略技术协变量。测序深度、细胞周期、批次都会影响表达。用 sc.pp.regress_out 或把协变量纳入模型。

单细胞的差异分析方法选择：

有生物学重复的实验，伪批量 + DESeq2 是最严谨的选择，因为它的统计单位是样本而非细胞。

9. 轨迹推断与数据整合

轨迹推断（trajectory inference） 分析细胞的分化/发育路径。细胞在分化过程中转录组连续变化，轨迹推断试图重建这条「路径」：

# PAGA：分析簇之间的连接关系
sc.tl.paga(adata, groups="leiden")
sc.pl.paga(adata)

# 扩散伪时间（diffusion pseudotime）
sc.tl.dpt(adata)
sc.pl.umap(adata, color=["dpt_pseudotime"])

轨迹推断的假设是「细胞状态连续变化」，适用于发育、分化、细胞周期等场景。不适用于离散的细胞类型（如成熟 T 细胞和 B 细胞之间没有「轨迹」）。常用工具还有 Monocle3、Slingshot。

数据整合（integration） 是处理「批次效应」的关键。当数据来自多个样本、多个实验批次时，同一细胞类型可能因为批次差异而分开。整合方法把不同批次对齐：

# Harmony 整合（快，效果好）
sc.external.pp.harmony_integrate(adata, key="batch")
sc.pp.neighbors(adata, use_rep="X_pca_harmony")

# BBKNN（基于近邻图）
sc.external.pp.bbknn(adata, batch_key="batch")

# scVI（深度生成模型，慢但强大）
import scvi
scvi.model.SCVI.setup_anndata(adata, batch_key="batch")
model = scvi.model.SCVI(adata); model.train()
adata.obsm["X_scVI"] = model.get_latent_representation()

整合的取舍：整合太弱，批次效应未消除，同一细胞类型分开；整合太强，可能把真实的生物学差异也抹掉（over-correction）。评估方法是看「同一细胞类型是否跨批次聚在一起」且「不同细胞类型是否仍分开」。

一个原则：能通过实验设计避免批次，就不要靠整合补救。但如果数据已经产生，整合是必要的。

权衡取舍

常见坑清单

1. 机械套用质控阈值：不同组织线粒体比例天然不同，误杀细胞；看分布定阈值。
2. 忽略双细胞：双细胞看起来像「表达两套标记」的异常类型，污染聚类；用 Scrublet 检测。
3. 从 UMAP 读定量结论：UMAP 距离不代表相似度、簇大小不代表细胞数；只看结构。
4. 把细胞当独立样本：伪重复导致假阳性；有重复时用伪批量。
5. 聚类分辨率固定：不同问题需要不同粒度；多分辨率对比再选。
6. 自动注释不审核：参考没有的细胞类型被强行归类；结合标记基因人工验证。
7. 过度插补：imputation 引入假信号；谨慎使用或不用。
8. 整合过度：把真实生物学差异也抹平；评估整合前后的标记基因保留情况。
9. 忽略细胞周期：增殖细胞因周期基因聚集，掩盖真实类型；必要时回归掉周期效应。
10. 稀疏矩阵稠密化：内存爆掉；全程用 scipy.sparse，避免 .toarray()。

小结

单细胞分析把转录组的分辨率推到了细胞层面，也把「统计严谨性」的要求推到了新高度。稀疏矩阵、零膨胀、伪重复、批次效应——每一个都是 bulk 分析中没有的挑战。工具（Scanpy/Seurat）让分析变得容易，但「参数选择」和「结果解读」仍高度依赖对生物学与统计的理解。

工程上，单细胞最该建立的三个习惯：一是「质控看分布」，每个数据集的组织、细胞类型不同，阈值必须自适应；二是「注释必验证」，自动工具的结果一定要用标记基因核对；三是「统计单位是样本」，有生物学重复时用伪批量而非细胞级检验。另一个实用建议是「保存中间对象」：单细胞分析步骤多、耗时长，把每个关键步骤的 adata 对象存盘（.h5ad），避免从头重跑。

下一步可以看 生信流程编排 了解如何把单细胞流程工程化，或 RNA-seq 转录组分析 对比 bulk 与单细胞的差异。如果你的聚类结果「全是同一个大簇」或「碎成几十个小簇」，先检查质控、HVG 数量、聚类分辨率与整合是否合适。