单细胞转录组分析 · 从入门到精通

本教材基于「生信学习库」知识库内容深度整理而成，涵盖：

• 📖 Single-cell Best Practices 全书 47 章详解
• 💻 Starlitnightly 18 篇实战教程完整代码
• 🔧 Scanpy / Bioconductor / omicverse 工具链全参考
• 🧬 多模态（ATAC/CITE-seq/空间组/TCR）扩展分析
• ⚙️ 工程化（Snakemake/Biopython/DeepChem/Linux/算法）

全书 33 章，约 13,000 行，可直接作为教科书使用。

第一篇：入门基础

1. 单细胞RNA测序技术原理

1.1 测序技术演进：一代→二代→三代

1.1.1 第一代测序——Sanger测序

第一代测序技术，又称Sanger测序（桑格测序），由Frederick Sanger于1977年发明，是基于双脱氧核苷酸末端终止法的测序技术。其基本原理是：在DNA聚合酶催化的测序反应中，掺入带有荧光标记的ddNTP（双脱氧核苷酸），当ddNTP被掺入到延伸的DNA链中时，由于ddNTP缺乏3’-OH基团，链的延伸被终止。通过毛细管电泳分离不同长度的DNA片段，即可读出DNA序列。

第一代测序的特点：

• 读长较长（700-1000 bp）
• 准确率极高（>99.99%）
• 通量极低，每次只能测定一条DNA序列
• 成本高昂（人类基因组计划花费约30亿美元）
• 适用于小规模、高精度需求的验证性测序

在单细胞领域的局限： Sanger测序通量太低，完全无法满足单细胞水平上对数千至数万个细胞进行转录组分析的需求。

1.1.2 第二代测序——高通量测序（NGS）

第二代测序技术（Next-Generation Sequencing, NGS），又称高通量测序，约在2005-2007年间商业化。代表性平台包括Illumina（Solexa）、Roche 454、SOLiD等，其中Illumina平台占据主导地位。

Illumina测序的基本原理（边合成边测序，Sequencing by Synthesis, SBS）：

1. 文库构建：将RNA逆转录为cDNA，片段化后加上接头（adaptor）
2. 桥式PCR扩增：在flow cell表面进行簇生成（cluster generation）
3. 测序：每次加入一种带荧光标记的dNTP，合成一个碱基后拍摄荧光信号，确定碱基类型
4. 去荧光：切除荧光基团，进入下一个循环

第二代测序的特点：

• 通量极高（一次运行可产生数十亿条reads）
• 单碱基成本极低（全基因组测序成本降至数百美元）
• 读长较短（Illumina: 75-300 bp; 已废弃的454: ~700 bp）
• 需要PCR扩增，可能引入扩增偏差
• 适用于大规模转录组测序，是单细胞RNA-seq的基础技术平台

在单细胞领域的应用： 当前几乎所有的scRNA-seq协议均基于Illumina平台进行测序。10x Genomics Chromium系统的测序深度通常为每个细胞50,000-100,000 reads。

1.1.3 第三代测序——长读长测序

第三代测序技术（Third-Generation Sequencing, TGS），又称长读长测序或单分子测序，不需要PCR扩增，直接对单条DNA/RNA分子进行测序。代表性平台包括PacBio（SMRT测序）和Oxford Nanopore（纳米孔测序）。

PacBio SMRT测序原理：

• 使用零模波导孔（Zero-Mode Waveguide, ZMW）技术
• 单个DNA聚合酶固定在ZMW底部
• 荧光标记的核苷酸掺入时产生荧光脉冲信号
• 平均读长10-20 kb，超长模式可达>50 kb
• HiFi模式准确率>99.9%

Oxford Nanopore测序原理：

• DNA/RNA分子穿过生物纳米孔（蛋白质通道）
• 不同碱基通过时引起电流变化
• 通过检测电流信号变化来识别碱基
• 读长可达数百万碱基
• 设备便携（MinION仅U盘大小）

第三代测序在单细胞领域的应用：

• 单细胞全长转录本测序（Iso-Seq），可检测可变剪接和转录本异构体
• 结合10x Genomics barcode的长读长测序
• 纳米孔直接RNA测序（无需逆转录，保留RNA修饰信息）

三代测序技术对比：

1.2 scRNA-seq协议分类：全长式 vs 标签式

scRNA-seq协议根据cDNA扩增和测序策略的不同，可分为两大类：全长式（Full-length）和标签式（Tag-based/UMI-based）。这一分类是最核心的技术区分，直接决定了下游数据处理方法的选择。

1.2.1 全长式协议（Full-length Protocols）

全长式协议对每个细胞中所有转录本的完整长度进行测序，不使用UMI（Unique Molecular Identifier，唯一分子标识符）。

Smart-seq2

Smart-seq2由Picelli等人于2013年开发，是目前最成熟的全长式scRNA-seq协议。

原理：

1. 单细胞通过FACS（荧光激活细胞分选）分选到96孔或384孔板中
2. 裂解细胞，释放RNA
3. 使用oligo-dT引物进行逆转录，合成第一链cDNA
4. 利用模板切换（Template Switching）机制：当逆转录酶到达mRNA的5’端时，会添加数个非模板的胞嘧啶（C），然后加入含有核糖鸟嘌呤（rG）的模板切换寡核苷酸（TSO），逆转录酶切换到TSO模板继续延伸
5. 通过PCR扩增全长的cDNA
6. 构建测序文库，进行Illumina测序

优势：

• 覆盖转录本全长，可检测可变剪接（Alternative Splicing）
• 可检测SNP和等位基因特异性表达
• 灵敏度高，对低表达基因检测能力更强
• 适用于转录本异构体分析

劣势：

• 不使用UMI，无法消除PCR扩增偏差
• 通量低（通常96-384个细胞/批次）
• 成本高（每个细胞成本约为10x的10-50倍）
• 手工操作步骤多，实验周期长

Smart-seq3

Smart-seq3由Hagemann-Jensen等人于2020年开发，是Smart-seq2的升级版本。

关键改进：

• 引入了5’端的UMI标记
• 结合了全长式和标签式的优势
• 通过5’端的UMI实现分子计数，同时保留全长转录本信息
• 灵敏度进一步提升

原理：

1. 与Smart-seq2类似的模板切换过程
2. 但TSO中包含了UMI序列
3. 因此每个cDNA分子的5’端带有唯一标识
4. 测序时先读取5’端的UMI，再读取转录本序列

CEL-seq / CEL-seq2

CEL-seq（Cell Expression by Linear amplification and Sequencing）由Hashimshony等人于2012年开发，是最早的标签式协议之一。

原理：

1. 使用FACS分选单细胞到孔板中
2. 使用带有细胞barcode和poly-dT的引物进行逆转录
3. 使用体外转录（IVT）进行线性扩增（而非PCR指数扩增）
4. 测序3’端

MARS-seq

MARS-seq（Massively Parallel RNA Single-cell sequencing）由Jaitin等人于2014年开发，同样基于IVT线性扩增和3’端测序。

1.2.2 标签式协议（Tag-based Protocols）

标签式协议仅对转录本的3’端或5’端进行测序，使用UMI和细胞barcode进行分子计数。

10x Genomics Chromium

10x Genomics Chromium是目前最广泛使用的scRNA-seq平台，于2016年推出。

原理：

1. 细胞悬液与含有barcode的凝胶珠（Gel Bead）在微流控芯片中混合
2. 每个凝胶珠表面带有约750,000条寡核苷酸引物，每条引物包含：
   • Illumina P5接头序列
   • 16 nt的细胞barcode（同一凝胶珠上的所有引物共享相同的barcode）
   • 10 nt的UMI（同一凝胶珠上的不同引物有不同的UMI）
   • 30 nt的poly-dT序列（用于捕获mRNA的poly-A尾）
3. 在微流控芯片中，单个细胞与单个凝胶珠被包裹在同一个液滴（GEM, Gel Bead-in-Emulsion）中
4. 液滴内的细胞被裂解，mRNA被带有barcode和UMI的引物捕获
5. 逆转录产生带有barcode和UMI的cDNA
6. 打破液滴，收集所有cDNA，混合进行PCR扩增和文库构建
7. 测序时读取3’端序列

核心优势：

• 超高通量：每次可捕获500-10,000个细胞（Chromium X可达数万）
• 使用UMI消除PCR扩增偏差
• 自动化程度高，实验流程标准化
• 可重复性好
• 成本相对较低

核心劣势：

• 仅测3’端，无法检测可变剪接和转录本异构体
• 对低表达基因的灵敏度低于Smart-seq2
• Dropout率较高（部分基因在某些细胞中未被检测到）
• 双细胞（Doublet）问题：两个细胞被包裹在同一个液滴中

Drop-seq

Drop-seq由Macosko等人于2015年开发，是10x Chromium的前驱技术。

原理：
与10x Chromium类似，也是基于液滴微流控技术：

1. 细胞与带有barcode的微珠在微流控装置中共包裹在液滴中
2. 每个微珠上带有约10^8条引物，包含细胞barcode和UMI
3. 液滴内裂解细胞、捕获mRNA、逆转录

与10x的差异：

• 使用普通微珠（非凝胶珠），细胞裂解和逆转录在液滴内完成
• 细胞捕获率较低
• 成本更低但数据质量略逊
• 是学术界的开源方案

inDrop

inDrop由Klein等人于2015年开发，同样基于液滴微流控技术。

原理与特点：

• 使用可降解的水凝胶微珠（degradable hydrogel beads）
• 微珠上的barcode引物可以在液滴内释放
• 使用IVT（体外转录）进行线性扩增，而非PCR
• 线性扩增可避免PCR偏差，但扩增效率较低

1.2.3 全长式 vs 标签式详细对比表

1.3 UMI原理与作用

1.3.1 什么是UMI

UMI（Unique Molecular Identifier，唯一分子标识符），又称Molecular Barcode，是在逆转录之前附加到每个原始RNA分子上的一段随机核苷酸序列（通常8-12 nt）。它的核心作用是区分来自同一原始分子的多个PCR拷贝和来自不同原始分子的cDNA。

1.3.2 为什么需要UMI

在scRNA-seq实验中，原始RNA分子需要经过以下步骤：

1. 逆转录（RT）：mRNA → cDNA
2. PCR扩增：将cDNA扩增到足够的量以进行文库构建和测序

PCR扩增是指数级的，不同的cDNA分子可能被扩增不同次数。假设原始分子A被扩增了100次，原始分子B被扩增了5次——测序reads中A的计数远高于B，但这并不反映原始mRNA的表达水平差异，而是PCR偏差。

没有UMI的情况：

• 测序reads数 = 原始分子数 × PCR扩增倍数
• 不同分子的扩增倍数不同，reads数无法直接反映原始分子数
• 需要依赖统计方法（如移位对数归一化）来校正

有UMI的情况：

• 每个原始分子在逆转录时被赋予唯一的UMI
• PCR扩增后，来自同一原始分子的所有reads共享相同的UMI
• 去重（deduplication）：将相同UMI的reads合并为一个计数
• UMI计数 = 原始分子数

1.3.3 UMI去重的数学原理

对于每个基因$g$在细胞$c$中的表达量，UMI计数过程可以建模为：

$$X_{gc} = \text{number of unique UMIs for gene } g \text{ in cell } c$$

而不是简单的reads计数。UMI去重的基本步骤：

1. 对每个细胞中每个基因的所有reads，按照UMI序列分组
2. 相同UMI的reads视为来自同一原始分子，只计为1
3. 不同UMI的reads视为来自不同原始分子

UMI碰撞（Collision）问题： 由于UMI长度有限（如10 nt = 4^10 = 1,048,576种可能），当同一基因在同一细胞中的原始分子数很多时，不同的原始分子可能被分配到相同的UMI，导致计数低估。这被称为UMI碰撞或”生日悖论”问题。

校正公式（基于生日悖论）：

$$\hat{N} = -M \ln\left(1 - \frac{n}{M}\right)$$

其中$M$是UMI空间大小（$4^L$，$L$为UMI长度），$n$是观测到的不同UMI数，$\hat{N}$是估计的原始分子数。

1.3.4 UMI的优势总结

1.4 液滴式 vs 板式完整对比

1.4.1 液滴式（Droplet-based）

代表平台： 10x Genomics Chromium, Drop-seq, inDrop

原理： 利用微流控技术将单个细胞与单个barcode微珠包裹在水-油乳液的液滴中。每个液滴作为一个独立的反应体系，其中的细胞被裂解，mRNA被带有特异性barcode的引物捕获。

流程（以10x Chromium为例）：

1. 制备单细胞悬液（要求细胞活性>80%，无团块）
2. 在Chromium Controller上将细胞悬液与凝胶珠和试剂混合
3. 微流控芯片生成GEM液滴（约100,000个液滴/运行）
4. GEM内裂解→逆转录→cDNA带barcode+UMI
5. 打破GEM，收集cDNA，PCR扩增
6. 构建Illumina测序文库

优势：

• 通量极高（每运行可处理500-10,000+细胞）
• 实验流程高度自动化和标准化
• 单细胞成本低
• 适合发现稀有细胞类型
• 适合构建大规模细胞图谱

劣势：

• 仅测3’端或5’端
• Dropout率较高
• 双细胞问题（约0.8%-8%的液滴包含两个细胞）
• 环境RNA（Ambient RNA/Soup）污染
• 对细胞悬液质量要求高

1.4.2 板式（Plate-based）

代表平台： Smart-seq2/3, CEL-seq2, MARS-seq

原理： 通过FACS（荧光激活细胞分选）或手动操作将单个细胞分选到96孔或384孔板中，每个孔作为独立的反应体系。

流程（以Smart-seq2为例）：

1. 制备单细胞悬液
2. 使用FACS将单个细胞分选到含裂解液的96/384孔板中
3. 在孔中进行裂解→逆转录（含模板切换）→PCR扩增
4. 每孔的cDNA单独进行文库构建（或使用indexing PCR合并）
5. Illumina测序

优势：

• 全长转录本覆盖（Smart-seq2/3）
• 可检测可变剪接和转录本异构体
• 灵敏度高，对低表达基因检测更好
• Dropout率相对较低
• 无双细胞问题（除非FACS分选错误）
• 可结合FACS表面标记信息进行预筛选

劣势：

• 通量极低（通常<1000细胞/实验）
• 实验操作繁琐，手工步骤多
• 单细胞成本高
• 不使用UMI（Smart-seq2），PCR偏差无法校正
• 难以发现稀有细胞类型（通量限制）

1.4.3 液滴式 vs 板式详细对比表

1.5 每种技术的适用场景

1.5.1 场景推荐表

1.5.2 选择决策树

1
2
3
4
5
6
7
8
9
10

你的研究目标是什么？
├── 需要全长转录本信息？
│   ├── 是 → Smart-seq2（无UMI需求）/ Smart-seq3（需UMI）
│   └── 否 → 继续判断
│       ├── 需要高通量（>1000细胞）？
│       │   ├── 是 → 10x Chromium
│       │   └── 否 → CEL-seq2 或 MARS-seq
│       └── 预算有限？
│           ├── 是 → Drop-seq（开源、低成本）
│           └── 否 → 10x Chromium（标准化、可重复性好）

1.6 关键参考文献

1. Sanger测序： Sanger, F., Nicklen, S. & Coulson, A.R. (1977). DNA sequencing with chain-terminating inhibitors. Proceedings of the National Academy of Sciences, 74(12), 5463-5467.
2. Illumina测序： Bentley, D.R. et al. (2008). Accurate whole human genome sequencing using reversible terminator chemistry. Nature, 456(7218), 53-59.
3. Smart-seq2： Picelli, S. et al. (2013). Smart-seq2 for sensitive full-length transcriptome profiling in single cells. Nature Methods, 10(11), 1096-1098.
4. Smart-seq3： Hagemann-Jensen, M. et al. (2020). Single-cell RNA counting at allele and isoform resolution using Smart-seq3. Nature Biotechnology, 38(6), 708-714.
5. 10x Genomics Chromium： Zheng, G.X.Y. et al. (2017). Massively parallel digital transcriptional profiling of single cells. Nature Communications, 8, 14049.
6. Drop-seq： Macosko, E.Z. et al. (2015). Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell, 161(5), 1202-1214.
7. inDrop： Klein, A.M. et al. (2015). Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell, 161(5), 1187-1201.
8. CEL-seq： Hashimshony, T. et al. (2012). CEL-seq: single-cell RNA-seq by multiplexed linear amplification. Cell Reports, 2(3), 666-673.
9. MARS-seq： Jaitin, D.A. et al. (2014). Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types. Science, 343(6172), 776-779.
10. UMI原理： Kivioja, T. et al. (2011). Counting absolute numbers of molecules using unique molecular identifiers. Nature Methods, 9(1), 72-74.
11. scRNA-seq综述： Aldridge, S. & Teichmann, S.A. (2020). Single cell transcriptomics comes of age. Nature Communications, 11, 4307.
12. scRNA-seq技术演进： Wen, L. & Tang, F. (2022). The evolution of single-cell RNA sequencing technology and application. Science China Life Sciences, 65, 838-847.

2. 原始数据处理

2.1 Cell Ranger完整流程

Cell Ranger是10x Genomics官方提供的配套分析软件，用于处理Chromium单细胞RNA-seq数据。它是一个高度集成的分析管线，可以从原始测序数据（BCL格式）直接生成基因表达矩阵和初步的分析结果。

Cell Ranger的核心流程包括三个主要步骤：mkfastq → count → aggr。

2.1.1 cellranger mkfastq：BCL→FASTQ转换

功能： 将Illumina测序仪产生的原始BCL（Base Call）文件转换为FASTQ格式的测序数据。该步骤借用了Illumina的bcl2fastq工具，但进行了优化以适配10x Genomics的测序结构。

原理：

1. Illumina测序仪输出的原始数据是BCL格式（每个循环每个cluster的碱基调用）
2. mkfastq将BCL文件按照样本index（I7/I5）拆分为不同样本的FASTQ文件
3. 对于10x Genomics的测序文库，通常一个样本会产生3个FASTQ文件：
   • I1（Index read）：包含样本index序列
   • R1（Read 1）：包含26 bp的cell barcode + UMI序列
   • R2（Read 2）：包含98 bp的cDNA插入片段序列

命令示例：

1
2
3
4
5
6
7
8
9
10
11
12

# 基本用法
cellranger mkfastq \
  --run=/path/to/bcl_folder \
  --output-dir=/path/to/fastq_output \
  --csv=sample_sheet.csv \
  --localcores=16 \
  --localmem=64

# sample_sheet.csv 格式示例：
# Lane,Sample,Index
# 1,Sample1,SI-GA-A1
# 2,Sample2,SI-GA-B1

参数说明：

2.1.2 cellranger count：FASTQ→计数矩阵

功能： 对单个样本的FASTQ文件进行比对、过滤、UMI计数和基因表达定量，生成gene-barcode计数矩阵。

完整流程：

1. 比对（Alignment）

   • 使用STAR（Spliced Transcripts Alignment to a Reference）将R2 reads比对到参考基因组
   • Cell Ranger内置了修改版的STAR，针对10x数据进行了优化
   • 比对时考虑了剪接位点，可以将跨越内含子的reads正确比对

2. Barcode校正（Barcode Correction）

   • 对R1中读取的26 bp序列进行解析：
     • 前16 bp为cell barcode（与已知的737,000个10x barcode whitelist比对）
     • 后10 bp为UMI
   • 对cell barcode进行容错校正：如果读到的barcode与whitelist中的某个barcode只有1个碱基不同，则校正为whitelist中的barcode

3. UMI计数（UMI Counting）

   • 对于比对到同一基因、具有相同cell barcode的reads
   • 按UMI分组，相同UMI视为来自同一原始分子
   • 对UMI也进行容错校正（编辑距离≤1的UMI合并）
   • 最终得到每个基因在每个细胞中的UMI计数

4. 细胞呼叫（Cell Calling）

   • 区分真实细胞和空液滴（不含细胞但可能含有环境RNA）
   • Cell Ranger使用空液滴模型：先假设所有barcode都是空液滴，拟合背景分布
   • 然后识别总UMI数显著高于背景的barcode为真实细胞
   • 默认方法：对每个barcode的总UMI数取log10排序，寻找”拐点”（knee point）

5. 生成输出文件

   • filtered_feature_bc_matrix/：过滤后的基因表达矩阵（仅包含真实细胞）
   • raw_feature_bc_matrix/：未过滤的基因表达矩阵（包含所有barcode）
   • metrics_summary.csv：质控指标汇总
   • web_summary.html：交互式HTML报告
   • molecule_info.h5：分子级别信息（用于后续aggr）

命令示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 基本用法
cellranger count \
  --id=sample1 \
  --transcriptome=/path/to/refdata-gex-GRCh38-2020-A \
  --fastqs=/path/to/fastq_output/sample1 \
  --sample=sample1 \
  --localcores=16 \
  --localmem=64 \
  --expect-cells=5000

# 参数说明：
# --id: 输出目录名称
# --transcriptome: 参考基因组路径（需预先用cellranger mkref构建）
# --fastqs: FASTQ文件所在目录
# --sample: 样本名称（与mkfastq输出一致）
# --expect-cells: 预期细胞数（帮助细胞呼叫算法）
# --localcores: CPU核心数
# --localmem: 内存（GB）

构建参考基因组：

1
2
3
4
5
6
7
8
9

# 下载10x参考基因组
# 人类：refdata-gex-GRCh38-2020-A
# 小鼠：refdata-gex-mm10-2020-A

# 或自定义参考基因组
cellranger mkref \
  --genome=custom_ref \
  --fasta=genome.fa \
  --genes=genes.gtf

2.1.3 cellranger aggr：多样本聚合

功能： 将多个cellranger count的输出结果进行整合，生成合并的基因表达矩阵，并进行测序深度的归一化。

为什么需要aggr：

• 同一样本可能在不同lane上测序
• 同一组的不同样本需要合并分析
• 需要消除不同文库之间测序深度的差异

归一化策略：

• 对每个文库计算一个归一化因子
• 默认基于映射到转录组的reads数进行归一化
• 公式：$f_i = \frac{M}{N_i}$，其中$M$是所有文库的中位有效reads数，$N_i$是第$i$个文库的有效reads数
• 归一化后的UMI计数：$x’{ij} = \lfloor f_i \times x{ij} \rfloor$

命令示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 创建aggregation CSV文件 (aggregation.csv)
# library_id,chemical,lanes,molecule_h5
# sample1,three_prime_v3,1,/path/to/sample1/outs/molecule_info.h5
# sample2,three_prime_v3,1,/path/to/sample2/outs/molecule_info.h5

# 运行aggr
cellranger aggr \
  --id=pooled_analysis \
  --csv=aggregation.csv \
  --normalize=map \
  --localcores=32 \
  --localmem=128

# 参数说明：
# --id: 输出目录名称
# --csv: 聚合信息表
# --normalize: 归一化方法（map=基于映射reads数；none=不归一化）

注意事项：

• 所有文库必须使用相同版本的参考基因组
• chemical列需要注明化学版本（three_prime_v2, three_prime_v3等）
• 聚合后应重新进行下游分析（归一化、聚类等）

2.2 alevin-fry流程（作为替代方案）

2.2.1 alevin-fry简介

alevin-fry是一个高效、灵活的单细胞测序数据处理工具，由Srivastava等人于2022年发表在Nature Methods上。它是salmon软件生态的一部分，可以作为Cell Ranger的替代方案，具有速度更快、内存占用更少的优势。

核心优势：

• 处理速度比Cell Ranger快5-10倍
• 内存占用极低（可处理百万级细胞数据）
• 支持多种定量模式（pseudoalignment-based）
• 开源免费，不依赖10x Genomics许可
• 输出格式与Cell Ranger兼容

2.2.2 alevin-fry完整流程

alevin-fry的处理流程分为三步：salmon alevin → alevin-fry generate-permit-list → alevin-fry quantify

步骤1：使用salmon alevin进行伪比对

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# 构建salmon索引
salmon index \
  -t transcriptome.fa \
  -i salmon_index \
  -k 31

# 运行salmon alevin（伪比对+初始定量）
salmon alevin \
  -l ISR \
  -i salmon_index \
  -1 R1.fastq.gz \
  -2 R2.fastq.gz \
  -o alevin_output \
  -p 16 \
  --chromiumV3

# 参数说明：
# -l ISR: 文库类型（I=反向，S=有义链，R=配对）
# -i: salmon索引路径
# -1: R1 FASTQ（含barcode+UMI）
# -2: R2 FASTQ（含cDNA序列）
# -o: 输出路径
# -p: 线程数
# --chromiumV3: 指定10x Chromium v3化学方案

步骤2：生成允许列表（permit list）

1
2
3
4
5
6
7
8
9
10
11
12

# 识别真实细胞的barcode
alevin-fry generate-permit-list \
  -i alevin_output \
  -o alevin_fry_output \
  --expected-ori 5000 \
  --knee-dist

# 参数说明：
# -i: salmon alevin的输出路径
# -o: alevin-fry的输出路径
# --expected-ori: 预期细胞数
# --knee-dist: 使用knee-point方法识别细胞

步骤3：定量

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# UMI去重和定量
alevin-fry quantify \
  -i alevin_fry_output \
  -o alevin_fry_output/quant \
  -m salmon_index \
  --use-mtx \
  --rad

# 参数说明：
# -i: alevin-fry输出路径
# -o: 定量结果输出路径
# -m: 映射文件（来自salmon索引）
# --use-mtx: 以MTX格式输出计数矩阵
# --rad: 使用resolution-aware deduplication

输出文件：

• alevin_fry_output/quant/alevin/featureDump.txt：质控指标
• alevin_fry_output/quant/alevin/quants_mat.mtx：计数矩阵
• alevin_fry_output/quant/alevin/quants_mat_rows.txt：细胞barcode
• alevin_fry_output/quant/alevin/quants_mat_cols.txt：基因列表

2.2.3 将alevin-fry结果导入Python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

import scanpy as sc
import pandas as pd
from scipy.io import mmread
from scipy.sparse import csr_matrix

# 读取计数矩阵
matrix = csr_matrix(mmread("alevin_fry_output/quant/alevin/quants_mat.mtx"))
barcodes = pd.read_csv("alevin_fry_output/quant/alevin/quants_mat_rows.txt", header=None)[0].values
genes = pd.read_csv("alevin_fry_output/quant/alevin/quants_mat_cols.txt", header=None)[0].values

# 创建AnnData对象
adata = sc.AnnData(
    X=matrix.T,  # 转置：行为细胞，列为基因
    obs=pd.DataFrame(index=barcodes),
    var=pd.DataFrame(index=genes)
)

print(adata)
# AnnData object with n_obs × n_vars = 5000 × 33538

2.3 FASTQ→计数矩阵的完整转化过程

2.3.1 数据流转全景

1
2
3
4
5
6
7
8
9
10
11

原始测序数据 (BCL)
    ↓ [cellranger mkfastq / bcl2fastq]
FASTQ文件 (R1: barcode+UMI, R2: cDNA, I1: sample index)
    ↓ [cellranger count / salmon alevin]
比对到参考基因组 (SAM/BAM)
    ↓ [Barcode校正 + UMI去重]
基因-细胞计数矩阵 (Gene × Cell count matrix)
    ↓ [细胞呼叫/过滤]
过滤后的计数矩阵 (filtered count matrix)
    ↓ [导入scanpy/Seurat]
AnnData / SingleCellExperiment对象

2.3.2 每一步的详细说明

Step 1: BCL → FASTQ

• 将Illumina原始信号转换为碱基序列
• 按样本index拆分
• 输出paired-end FASTQ文件

Step 2: FASTQ → Alignment

• 将R2 reads比对到参考基因组
• Cell Ranger使用STAR（splicing-aware比对）
• alevin-fry使用伪比对（pseudoalignment，更快但不提供精确比对坐标）

Step 3: Alignment → Gene Assignment

• 将比对结果分配到基因
• 基于GTF注释文件确定read属于哪个基因
• 对于比对到多个基因的reads，分配给唯一基因或丢弃

Step 4: Barcode Processing + UMI Deduplication

• 解析R1中的cell barcode和UMI
• 校正barcode中的测序错误
• 对同一基因、同一细胞、同一UMI的reads去重

Step 5: Cell Calling

• 区分真实细胞和空液滴
• 基于总UMI数的分布，识别显著高于背景的barcode

Step 6: Output Count Matrix

• 输出稀疏矩阵格式（MTX）
• 三个文件：矩阵、行索引、列索引

2.4 比对策略：STAR vs Cell Ranger

推荐策略：

• 一般分析：alevin-fry（速度快、资源消耗少）
• 需要BAM文件进行下游分析（如变异检测）：Cell Ranger
• 正式发表：Cell Ranger（标准化程度高，审稿人更熟悉）

2.5 多文库聚合策略

2.5.1 为什么需要聚合

• 一个实验可能包含多个样本
• 同一样本可能在不同lane上测序
• 不同样本可能使用不同的化学方案（v2/v3）
• 需要消除批次效应和测序深度差异

2.5.2 聚合方法

方法1：Cell Ranger aggr（推荐）

• 在原始UMI计数层面进行聚合
• 基于映射reads数进行归一化
• 优点：保留了UMI计数的整数性
• 缺点：需要molecule_info.h5文件

方法2：下游分析时合并

• 分别对每个样本运行cellranger count
• 在Python中使用sc.concat或anndata.concat合并
• 需要后续进行批次校正（如Harmony、scVI等）

1
2
3
4
5
6
7
8
9
10
11
12
13

import scanpy as sc

# 方法2：在Python中合并多个样本
adata1 = sc.read_10x_mtx('sample1/filtered_feature_bc_matrix/')
adata1.obs['sample'] = 'sample1'

adata2 = sc.read_10x_mtx('sample2/filtered_feature_bc_matrix/')
adata2.obs['sample'] = 'sample2'

# 合并
adata = sc.concat([adata1, adata2], join='inner')
print(adata)
# AnnData object with n_obs × n_vars = 10000 × 33538

2.6 将Cell Ranger结果导入Python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import scanpy as sc

# 方法1：读取过滤后的计数矩阵
adata = sc.read_10x_mtx(
    'path/to/filtered_feature_bc_matrix/',
    var_names='gene_symbols',
    cache=True
)
print(adata)
# AnnData object with n_obs × n_vars = 5000 × 33538

# 方法2：读取h5格式
adata = sc.read_10x_h5(
    'path/to/filtered_feature_bc_matrix.h5'
)

# 方法3：读取未过滤的矩阵（用于SoupX环境RNA校正）
adata_raw = sc.read_10x_mtx(
    'path/to/raw_feature_bc_matrix/',
    var_names='gene_symbols'
)

# 保存为h5ad格式（后续分析使用）
adata.write_h5ad('raw_counts.h5ad')

2.7 注意事项

1. Cell Ranger版本一致性：同一项目的所有样本应使用相同版本的Cell Ranger，不同版本的结果可能存在细微差异
2. 参考基因组版本：确保使用与实验设计一致的参考基因组（如GRCh38 vs GRCh37）
3. 内存需求：Cell Ranger count通常需要32-64 GB内存，建议在服务器上运行
4. FASTQ文件完整性：运行前检查FASTQ文件大小和read数量是否正常
5. 细胞悬液质量：Cell Ranger的细胞呼叫算法依赖于细胞barcode的UMI分布，低质量悬液可能导致错误呼叫
6. 双细胞率：10x Chromium的双细胞率约为每1000个细胞0.8%，即8000个细胞中约有64个双细胞
7. 环境RNA污染：死细胞释放的mRNA会污染液滴，导致计数偏差，可通过SoupX校正
8. alevin-fry的定量模式：cr-like模式输出与Cell Ranger最接近的结果；full模式提供更精细的UMI去重
9. 多映射reads处理：Cell Ranger将多映射reads分配给最可能的基因，alevin-fry默认丢弃多映射reads
10. 数据备份：始终保留原始FASTQ文件和Cell Ranger完整输出

3. 数据结构

3.1 AnnData完整结构

AnnData（Annotated Data）是Python scverse生态系统的核心数据结构，专为存储和操作带注释的多维数据而设计。它是scanpy、scvi-tools、squidpy等单细胞分析工具的基础。

3.1.1 AnnData结构概览

1
2
3
4
5
6
7
8
9
10

AnnData Object
├── X              # 核心计数矩阵 (n_obs × n_vars, 通常为稀疏矩阵)
├── obs            # 观测（细胞）的一维元数据 (DataFrame, n_obs × k)
├── var            # 变量（基因）的一维元数据 (DataFrame, n_vars × k)
├── obsm           # 观测的多维注释 (Dict[str, ndarray], n_obs × m)
├── varm           # 变量的多维注释 (Dict[str, ndarray], n_vars × m)
├── uns            # 非结构化元数据 (Dict, 任意类型)
├── layers         # 与X对齐的额外矩阵 (Dict[str, ndarray])
├── obsp           # 观测的成对矩阵 (Dict[str, ndarray], n_obs × n_obs)
└── varp           # 变量的成对矩阵 (Dict[str, ndarray], n_vars × n_vars)

3.1.2 各属性详解

X — 核心表达矩阵

• 类型：scipy.sparse.csr_matrix 或 numpy.ndarray
• 维度：n_obs × n_vars（细胞 × 基因）
• 通常存储当前分析步骤使用的表达矩阵
• 可以是原始UMI计数、归一化后的值、缩放后的值等
• 建议始终在layers中保留原始计数

1
2
3
4

# 访问X
adata.X  # 完整矩阵
adata.X[0:5, 0:5]  # 前5个细胞和5个基因的子矩阵
adata.X.sum(axis=1)  # 每个细胞的总counts

obs — 细胞元数据

• 类型：pandas.DataFrame
• 行索引：细胞barcode
• 每列是一个细胞级别的注释信息

1
2
3
4
5
6
7
8
9
10
11
12

# 常见obs列
# - 'n_genes_by_counts': 每个细胞检测到的基因数
# - 'total_counts': 每个细胞的总UMI数
# - 'pct_counts_mt': 线粒体基因比例
# - 'cell_type': 细胞类型注释
# - 'batch': 批次信息
# - 'leiden': Leiden聚类结果
# - 'sample': 样本来源

# 添加细胞元数据
adata.obs['cell_type'] = ['CD4 T cell', 'B cell', 'Monocyte', ...]
adata.obs['condition'] = ['control', 'treatment', ...]

var — 基因元数据

• 类型：pandas.DataFrame
• 行索引：基因名/Ensembl ID
• 每列是一个基因级别的注释信息

1
2
3
4
5
6
7
8
9
10
11
12

# 常见var列
# - 'gene_ids': Ensembl基因ID
# - 'mt': 是否为线粒体基因
# - 'ribo': 是否为核糖体基因
# - 'n_cells_by_counts': 表达该基因的细胞数
# - 'mean_counts': 基因的平均counts
# - 'highly_variable': 是否为高可变基因
# - 'dispersions': 基因离散度

# 添加基因元数据
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))

obsm — 细胞的多维注释

• 类型：Dict[str, numpy.ndarray]
• 每个值的维度：n_obs × m（m为嵌入维度）
• 通常存储降维结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 常见obsm键
# - 'X_pca': PCA降维结果 (n_obs × n_pcs)
# - 'X_umap': UMAP嵌入 (n_obs × 2)
# - 'X_tsne': t-SNE嵌入 (n_obs × 2)
# - 'X_mde': MDE嵌入 (n_obs × 2)
# - 'X_diffmap': 扩散映射 (n_obs × n_components)

# 访问obsm
adata.obsm['X_pca']  # PCA结果矩阵
adata.obsm['X_umap']  # UMAP坐标

# omicverse特殊格式
adata.obsm['scaled|original|X_pca']  # 缩放后的PCA
adata.obsm['counts|original|X_pca']  # 基于原始计数的PCA

varm — 基因的多维注释

• 类型：Dict[str, numpy.ndarray]
• 每个值的维度：n_vars × m
• 较少使用，可存储基因的PCA载荷（loadings）等

1
2

# 示例：存储PCA载荷
adata.varm['PCs'] = pca_model.components_.T  # n_vars × n_pcs

uns — 非结构化元数据

• 类型：Dict
• 可存储任意Python对象
• 通常存储分析参数、颜色映射、邻域图信息等

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 常见uns键
# - 'hvg': 高可变基因选择参数和结果
# - 'neighbors': KNN邻域图参数
# - 'leiden': Leiden聚类参数
# - 'pca': PCA参数
# - 'rank_genes_groups': 差异表达基因结果

# 添加自定义元数据
adata.uns['experiment_info'] = {
    'organism': 'Homo sapiens',
    'tissue': 'PBMC',
    'platform': '10x Chromium v3',
    'date': '2024-01-15'
}

layers — 额外表达矩阵

• 类型：Dict[str, numpy.ndarray] 或 Dict[str, scipy.sparse.matrix]
• 每个值与X维度相同：n_obs × n_vars
• 存储不同处理阶段的表达矩阵

1
2
3
4
5
6
7
8

# 常见layers用法
adata.layers['counts'] = adata.X.copy()  # 保存原始UMI计数
adata.layers['log1p_norm'] = log_normalized  # 移位对数归一化
adata.layers['scran_norm'] = scran_normalized  # scran归一化
adata.layers['pearson_residuals'] = pearson  # Pearson残差

# 在不同layers之间切换
adata.X = adata.layers['counts'].copy()  # 恢复原始计数

obsp — 细胞成对矩阵

• 类型：Dict[str, numpy.ndarray] 或 Dict[str, scipy.sparse.matrix]
• 维度：n_obs × n_obs
• 存储细胞间的距离矩阵、连接性矩阵等

1
2
3
4
5
6
7

# 常见obsp键
# - 'distances': KNN距离矩阵
# - 'connectivities': KNN连接性矩阵（用于Leiden聚类）

# 访问
adata.obsp['distances']  # 稀疏的距离矩阵
adata.obsp['connectivities']  # 稀疏的连接性矩阵

varp — 基因成对矩阵

• 类型：Dict[str, numpy.ndarray]
• 维度：n_vars × n_vars
• 存储基因间的相关矩阵、调控网络等

1
2

# 示例：基因相关矩阵
adata.varp['gene_correlation'] = correlation_matrix

3.1.3 AnnData创建、切片、合并的完整代码

创建AnnData对象

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

import numpy as np
import pandas as pd
import anndata as ad
from scipy.sparse import csr_matrix

# ========== 方法1：从计数矩阵创建 ==========
# 创建稀疏计数矩阵（100个细胞 × 2000个基因）
counts = csr_matrix(np.random.poisson(1, size=(100, 2000)), dtype=np.float32)
adata = ad.AnnData(counts)
print(adata)
# AnnData object with n_obs × n_vars = 100 × 2000

# ========== 方法2：带完整注释创建 ==========
# 计数矩阵
n_obs = 1000  # 细胞数
n_vars = 5000  # 基因数
counts = csr_matrix(np.random.poisson(1, size=(n_obs, n_vars)), dtype=np.float32)

# 细胞元数据
obs = pd.DataFrame({
    'cell_type': np.random.choice(['T cell', 'B cell', 'Monocyte'], n_obs),
    'batch': np.random.choice(['batch1', 'batch2'], n_obs),
    'n_genes': np.random.randint(500, 3000, n_obs),
}, index=[f'cell_{i}' for i in range(n_obs)])

# 基因元数据
var = pd.DataFrame({
    'gene_id': [f'ENSG{i:09d}' for i in range(n_vars)],
    'mt': [gene.startswith('MT-') for gene in [f'gene_{i}' for i in range(n_vars)]],
}, index=[f'gene_{i}' for i in range(n_vars)])

# 降维结果
obsm = {
    'X_pca': np.random.randn(n_obs, 50).astype(np.float32),
    'X_umap': np.random.randn(n_obs, 2).astype(np.float32),
}

# 非结构化元数据
uns = {
    'experiment': 'PBMC_analysis',
    'organism': 'Homo sapiens',
}

# 额外表达矩阵层
layers = {
    'counts': counts.copy(),
    'normalized': csr_matrix(np.log1p(counts.toarray()), dtype=np.float32),
}

# 创建AnnData对象
adata = ad.AnnData(
    X=counts,
    obs=obs,
    var=var,
    obsm=obsm,
    uns=uns,
    layers=layers,
)

print(adata)
# AnnData object with n_obs × n_vars = 1000 × 5000
#     obs: 'cell_type', 'batch', 'n_genes'
#     var: 'gene_id', 'mt'
#     obsm: 'X_pca', 'X_umap'
#     uns: 'experiment', 'organism'
#     layers: 'counts', 'normalized'

切片（Slicing）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# ========== 细胞切片 ==========
# 通过索引切片
adata_subset = adata[:500]  # 前500个细胞
print(adata_subset.shape)  # (500, 5000)

# 通过布尔掩码切片
mask = adata.obs['cell_type'] == 'T cell'
adata_t_cells = adata[mask]
print(adata_t_cells.shape)  # (~333, 5000)

# 通过细胞名切片
adata_selected = adata[['cell_0', 'cell_1', 'cell_2']]

# ========== 基因切片 ==========
# 通过基因名切片
adata_genes = adata[:, ['gene_0', 'gene_1', 'gene_2']]
print(adata_genes.shape)  # (1000, 3)

# 通过布尔掩码切片（如选择高可变基因）
mask_hvg = adata.var['highly_variable'] == True  # 假设已计算
adata_hvg = adata[:, mask_hvg]

# ========== 同时切片 ==========
adata_sub = adata[:500, :1000]  # 前500个细胞，前1000个基因

# ========== 注意：切片返回视图（View） ==========
print(adata_subset.is_view)  # True
# 修改视图会影响原始对象
# 若需独立副本：
adata_copy = adata_subset.copy()
print(adata_copy.is_view)  # False

合并（Merging/Concatenating）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

import anndata as ad

# ========== 沿观测（细胞）维度合并 ==========
adata1 = ad.AnnData(np.random.randn(100, 2000))
adata1.obs['batch'] = 'batch1'
adata1.obs_names = [f'cell_b1_{i}' for i in range(100)]
adata1.var_names = [f'gene_{i}' for i in range(2000)]

adata2 = ad.AnnData(np.random.randn(150, 2000))
adata2.obs['batch'] = 'batch2'
adata2.obs_names = [f'cell_b2_{i}' for i in range(150)]
adata2.var_names = [f'gene_{i}' for i in range(2000)]

# 内连接（仅保留共同基因）
adata_concat = ad.concat([adata1, adata2], join='inner')
print(adata_concat.shape)  # (250, 2000)

# 外连接（保留所有基因，缺失值填充为0）
adata_concat_outer = ad.concat([adata1, adata2], join='outer')

# 保留批次标签
adata_concat = ad.concat(
    [adata1, adata2],
    join='inner',
    label='batch',
    keys=['batch1', 'batch2'],
    merge='same'  # 合并策略：相同值保留，不同值设为缺失
)

# ========== 沿变量（基因）维度合并 ==========
adata_genes1 = ad.AnnData(np.random.randn(100, 1000))
adata_genes1.var_names = [f'gene_{i}' for i in range(1000)]

adata_genes2 = ad.AnnData(np.random.randn(100, 500))
adata_genes2.var_names = [f'gene_{1000+i}' for i in range(500)]

adata_var_concat = ad.concat(
    [adata_genes1, adata_genes2],
    axis=1,  # 沿变量维度合并
    join='inner'
)
print(adata_var_concat.shape)  # (100, 1500)

# ========== 合并时的merge参数 ==========
# merge='same': 相同值保留，不同值设为缺失（默认）
# merge='unique': 仅保留唯一值
# merge='first': 保留第一个对象的值
# merge=None: 不合并元数据

读写AnnData对象

1
2
3
4
5
6
7
8
9
10
11
12
13

# 保存为h5ad文件
adata.write_h5ad('my_data.h5ad', compression='gzip')

# 读取h5ad文件
adata_loaded = ad.read_h5ad('my_data.h5ad')

# 大文件可使用backed模式（按需读取，不全部载入内存）
adata_backed = ad.read_h5ad('my_data.h5ad', backed='r')
# 注意：backed模式下只能读取，不能修改

# 从10x Genomics输出读取
adata = sc.read_10x_mtx('path/to/filtered_feature_bc_matrix/')
adata = sc.read_10x_h5('path/to/filtered_feature_bc_matrix.h5')

3.2 MuData多模态数据结构

MuData是scverse生态中用于存储多模态单细胞数据的数据结构，建立在AnnData之上。它允许在同一个对象中存储来自不同组学层次的数据（如RNA + ATAC + 蛋白质）。

3.2.1 MuData结构

1
2
3
4
5
6
7
8
9
10
11
12

MuData Object
├── mod            # 模态字典 (Dict[str, AnnData])
│   ├── 'rna'      # scRNA-seq数据 (AnnData)
│   ├── 'atac'     # scATAC-seq数据 (AnnData)
│   └── 'prot'     # CITE-seq蛋白质数据 (AnnData)
├── obs            # 全局细胞元数据 (DataFrame)
├── var            # 全局变量元数据 (DataFrame)
├── obsm           # 全局细胞嵌入 (Dict)
├── varm           # 全局变量注释 (Dict)
├── uns            # 全局非结构化元数据 (Dict)
├── obsp           # 全局细胞成对矩阵 (Dict)
└── varp           # 全局变量成对矩阵 (Dict)

3.2.2 MuData代码示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

from mudata import MuData
import anndata as ad
import numpy as np
import pandas as pd
from scipy.sparse import csr_matrix

# 创建RNA模态数据
n_obs = 1000
rna_data = ad.AnnData(
    X=csr_matrix(np.random.poisson(1, size=(n_obs, 20000)), dtype=np.float32),
    obs=pd.DataFrame(index=[f'cell_{i}' for i in range(n_obs)]),
    var=pd.DataFrame(index=[f'gene_{i}' for i in range(20000)]),
)
rna_data.obs['n_genes'] = np.random.randint(500, 3000, n_obs)

# 创建蛋白质模态数据（CITE-seq）
prot_data = ad.AnnData(
    X=csr_matrix(np.random.poisson(5, size=(n_obs, 100)), dtype=np.float32),
    obs=pd.DataFrame(index=[f'cell_{i}' for i in range(n_obs)]),
    var=pd.DataFrame(index=[f'protein_{i}' for i in range(100)]),
)

# 创建MuData对象
mdata = MuData({'rna': rna_data, 'prot': prot_data})
print(mdata)
# MuData object with n_obs × n_vars = 1000 × 20100
#   2 modalities
#     rna:  1000 x 20000
#     prot: 1000 x 100

# 访问单个模态
print(mdata['rna'])  # 返回RNA的AnnData对象

# 更新MuData（在修改模态后）
mdata.update()

# 保存和读取
mdata.write('my_mudata.h5mu')
mdata_loaded = MuData.read('my_mudata.h5mu')

# backed模式读取
mdata_backed = MuData.read('my_mudata.h5mu', backed=True)

3.3 SingleCellExperiment（R）数据结构

SingleCellExperiment（SCE）是R/Bioconductor生态系统中存储单细胞数据的核心数据结构，类似于Python中的AnnData。

3.3.1 SCE结构

1
2
3
4
5
6
7
8
9
10
11
12
13

SingleCellExperiment Object
├── assays           # 表达矩阵列表 (List[matrix])
│   ├── 'counts'     # 原始计数
│   └── 'logcounts'  # 归一化后的值
├── colData          # 细胞元数据 (DataFrame, n_cells × k)
├── rowData          # 基因元数据 (DataFrame, n_genes × k)
├── reducedDims      # 降维结果 (List[matrix])
│   ├── 'PCA'        # PCA结果
│   └── 'UMAP'       # UMAP结果
├── colPairs         # 细胞成对数据 (List)
├── rowPairs         # 基因成对数据 (List)
├── metadata         # 非结构化元数据 (List)
└── alternativeExperiments  # 替代实验数据

3.3.2 SCE代码示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

# 创建SCE对象
library(SingleCellExperiment)

# 计数矩阵
counts <- matrix(rpois(10000, lambda = 1), nrow = 200, ncol = 50)
rownames(counts) <- paste0("gene_", 1:200)
colnames(counts) <- paste0("cell_", 1:50)

# 细胞元数据
cell_meta <- DataFrame(
    cell_type = sample(c("T cell", "B cell", "Monocyte"), 50, replace = TRUE),
    batch = sample(c("batch1", "batch2"), 50, replace = TRUE)
)

# 基因元数据
gene_meta <- DataFrame(
    gene_id = paste0("ENSG", sprintf("%09d", 1:200)),
    mt = grepl("^MT-", rownames(counts))
)

# 创建SCE对象
sce <- SingleCellExperiment(
    assays = list(counts = counts),
    colData = cell_meta,
    rowData = gene_meta
)

# 添加归一化矩阵
logcounts(sce) <- log1p(counts)

# 添加降维结果
reducedDim(sce, "PCA") <- prcomp(t(logcounts(sce)))$x[, 1:30]

# 访问数据
counts(sce)          # 原始计数
logcounts(sce)       # 归一化值
colData(sce)         # 细胞元数据
rowData(sce)         # 基因元数据
reducedDim(sce, "PCA")  # PCA结果

3.4 Python/R互操作

在单细胞分析中，Python（scanpy）和R（Seurat/SingleCellExperiment）是最常用的两大分析生态。互操作性允许研究者在两种语言之间灵活切换，利用各自的优势工具。

3.4.1 数据格式映射表

3.4.2 转换工具

方法1：zellkonverter（R → Python）

1
2
3
4
5
6
7
8
9

# 在R中安装
BiocManager::install("zellkonverter")

# 将SCE保存为h5ad格式
library(zellkonverter)
writeH5AD(sce, "my_data.h5ad")

# 从h5ad格式读取
sce <- readH5AD("my_data.h5ad")

方法2：anndata2ri（Python → R，通过rpy2）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# 在Python中安装
# pip install anndata2ri

import anndata2ri
from rpy2.robjects import globalenv
import anndata as ad

# 激活anndata2ri转换
anndata2ri.activate()

# 读取AnnData对象
adata = ad.read_h5ad('my_data.h5ad')

# 将AnnData传递到R环境
globalenv['adata'] = adata

# 在R中（通过rpy2的%%R magic）：
# library(SingleCellExperiment)
# sce <- adata  # 自动转换为SCE对象

方法3：MuDataSeurat（MuData ↔ Seurat）

1
2
3
4
5
6
7
8
9

# 安装
remotes::install_github("zilongstudio/MuDataSeurat")

# 从h5ad/h5mu读取到Seurat
library(MuDataSeurat)
seu <- ReadH5AD("my_data.h5ad")

# 从h5mu读取多模态数据到Seurat
seu <- ReadH5MU("my_data.h5mu")

方法4：sceasy

1
2
3
4
5
6
7
8
9
10
11

# pip install sceasy

import sceasy
import anndata as ad

# AnnData → Seurat (保存为RDS)
adata = ad.read_h5ad('my_data.h5ad')
sceasy.convert(adata, 'my_data.rds', output_file_type='anndata2seurat')

# Seurat → AnnData
sceasy.convert('my_data.rds', 'my_data_from_r.h5ad', output_file_type='seurat2anndata')

方法5：直接在Python中调用R函数（rpy2）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
from rpy2.robjects.conversion import localconverter

# 在Python中调用R的scran归一化
import anndata2ri
anndata2ri.activate()

# 准备数据
adata_pp = adata.copy()
sc.pp.normalize_total(adata_pp)
sc.pp.log1p(adata_pp)
sc.pp.pca(adata_pp)
sc.pp.neighbors(adata_pp)

# 传递数据到R
data_mat = adata_pp.X.T
if hasattr(data_mat, 'toarray'):
    data_mat = data_mat.tocsc().toarray()

with localconverter(ro.default_converter + pandas2ri.converter):
    ro.globalenv["data_mat"] = data_mat

# 在R中运行scran
ro.r('''
library(scran)
size_factors <- sizeFactors(
    computeSumFactors(
        SingleCellExperiment(list(counts=data_mat)),
        min.mean = 0.1
    )
)
''')
size_factors = ro.r('size_factors')

3.4.3 注意事项

1. 数据类型一致性：R和Python对稀疏矩阵的支持不同，转换时注意矩阵格式
2. 行/列方向：R中表达矩阵通常是基因×细胞，Python中是细胞×基因，需要转置
3. 基因名规范：确保基因名在Python和R中一致（如使用Ensembl ID而非symbol）
4. 版本兼容：zellkonverter和anndata2ri的版本需要与anndata包版本匹配
5. 内存管理：大对象在Python和R之间传递时需要双倍内存
6. 大文件处理：对于超大数据集，建议通过磁盘文件（h5ad/h5mu）进行转换，而非内存传递

3.5 关键参考文献

1. AnnData： Virshup, I. et al. (2023). The scverse project provides a computational ecosystem for single-cell omics data analysis. Nature Biotechnology, 41(5), 604-606.
2. Scanpy： Wolf, F.A., Angerer, P. & Theis, F.J. (2018). SCANPY: large-scale single-cell gene expression data analysis. Genome Biology, 19, 15.
3. MuData： Danilenko, A. et al. (2022). MuData: multimodal data object for Python. Journal of Open Source Software, 7(73), 4176.
4. SingleCellExperiment： Amezquita, R.A. et al. (2020). Orchestrating single-cell analysis with Bioconductor. Nature Methods, 17(2), 137-145.
5. Cell Ranger： 10x Genomics. Cell Ranger - Official Documentation. https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger
6. alevin-fry： Srivastava, A. et al. (2022). Alevin-fry unlocks rapid, accurate and memory-frugal quantification of single-cell RNA-seq data. Nature Methods, 19, 316-322.
7. STAR： Dobin, A. et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29(1), 15-21.
8. zellkonverter： Zappia, L. & Lun, A. Converting single-cell data structures between Bioconductor and Python. Bioconductor Vignette.

第二篇：预处理全流程

4. 质量控制

4.1 三类数据缺陷

单细胞RNA测序数据由于技术限制，不可避免地存在三类主要缺陷，质量控制的目的就是识别和去除这些有问题的数据。

4.1.1 Drop-out（基因丢失）

定义： Drop-out是指基因在某些细胞中表达量不为零，但由于技术原因（捕获效率低、测序深度不足等）未被检测到，导致计数矩阵中出现大量零值。

成因：

• mRNA捕获效率有限（10x Chromium约10-15%）
• 逆转录效率有限（约30-50%）
• 测序深度不足
• PCR扩增偏差

影响：

• 计数矩阵稀疏（通常50-80%为零）
• 下游分析中假阴性率高
• 差异表达分析中统计功效降低
• 低表达基因可能被完全遗漏

特征： 在UMI计数数据中，Drop-out率与基因表达水平负相关——高表达基因的Drop-out率低，低表达基因的Drop-out率高。

4.1.2 死亡细胞/低质量细胞

定义： 在样本制备过程中，部分细胞可能已经死亡或正在死亡，这些细胞的转录组特征与正常活细胞不同，会引入噪声。

特征：

• 细胞膜破裂，胞质mRNA泄漏到环境中
• 线粒体mRNA相对富集（因为线粒体较大，不易泄漏）
• 总UMI数低
• 检测到的基因数少

为什么需要去除：

• 死亡细胞的转录组不代表真实的生物学状态
• 环境RNA（来自死亡细胞的mRNA）会污染其他细胞的计数
• 死亡细胞引入假变异，影响聚类和差异表达分析

识别指标：

• 高线粒体基因比例（>15-20%通常指示死亡细胞）
• 低总UMI数
• 低检测基因数

4.1.3 双细胞（Doublets）

定义： 在液滴式scRNA-seq中，两个或多个细胞可能被包裹在同一个液滴中，被赋予同一个cell barcode，导致测序数据实际上是多个细胞的混合信号。

产生原因：

• 细胞悬液浓度过高
• 液滴包裹的随机性
• 细胞团块未完全解离

双细胞率估算：

• 10x Chromium约0.8% per 1000 cells
• 即8000个细胞中约有64个双细胞
• 公式：$\text{doublet rate} \approx 0.001 \times \frac{N_{\text{cells}}}{1000}$

双细胞类型：

• 同源双细胞（Homotypic doublet）： 两个相同类型的细胞。难以检测，但影响较小（看起来像一个正常的该类型细胞）
• 异源双细胞（Heterotypic doublet）： 两个不同类型的细胞。更容易检测，影响更大（看起来像一个中间态细胞，可能被误认为新的细胞类型）

影响：

• 可能被误认为新的细胞类型或中间态
• 影响差异表达分析
• 影响细胞类型注释

4.2 QC指标详解

质量控制通过计算每个细胞（barcode）的一系列QC指标来评估数据质量，常用的核心指标如下：

4.2.1 nUMIs（total_counts）

定义： 每个细胞中检测到的UMI总数，即该细胞的文库大小（library size）。

计算方法：

$$\text{nUMIs}c = \sum{g=1}^{G} x_{gc}$$

其中$x_{gc}$是细胞$c$中基因$g$的UMI计数。

生物学意义：

• 反映细胞的总RNA含量
• 低nUMIs可能指示低质量细胞、死亡细胞或空液滴
• 高nUMIs可能指示双细胞
• 不同细胞类型的nUMIs可能自然不同（如大细胞通常有更多RNA）

典型值：

• 10x Chromium v3：5,000-50,000 UMIs/细胞
• 低于500-1000 UMIs通常被认为是低质量细胞

代码示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import scanpy as sc

# 计算QC指标
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))
adata.var['hb'] = adata.var_names.str.contains('^HB[^(P)]')

sc.pp.calculate_qc_metrics(
    adata,
    qc_vars=['mt', 'ribo', 'hb'],
    percent_top=None,
    log1p=False,
    inplace=True
)

# nUMIs存储在 adata.obs['total_counts']
print(adata.obs['total_counts'].describe())
# count     14814.000000
# mean      12000.523456
# std        8000.123456
# min          50.000000
# 25%        7000.000000
# 50%       11000.000000
# 75%       16000.000000
# max      150000.000000

4.2.2 检测基因数（n_genes_by_counts / detected_genes）

定义： 每个细胞中表达量非零的基因数量。

计算方法：

$$\text{n_genes}c = \sum{g=1}^{G} \mathbb{1}(x_{gc} > 0)$$

生物学意义：

• 反映细胞的转录活性
• 低n_genes可能指示低质量细胞
• 高n_genes可能指示双细胞
• 不同细胞类型的基因表达谱宽度不同

典型值：

• 10x Chromium v3：1,000-5,000 genes/细胞
• 低于200-250 genes通常被认为是低质量细胞

4.2.3 线粒体基因比例（pct_counts_mt / mito_perc）

定义： 来自线粒体基因的UMI计数占细胞总UMI计数的比例。

计算方法：

$$\text{pct_counts_mt}c = \frac{\sum{g \in \text{MT}} x_{gc}}{\sum_{g=1}^{G} x_{gc}} \times 100$$

线粒体基因识别（人类）： 基因名以”MT-“开头（如MT-ND1, MT-CO1, MT-ATP6等，共13个蛋白编码基因）

线粒体基因识别（小鼠）： 基因名以”mt-“开头

生物学原理：

• 线粒体DNA编码的转录本位于线粒体内部
• 当细胞死亡时，细胞膜破裂，胞质mRNA泄漏
• 但线粒体较大且有自己的膜结构，mRNA不易泄漏
• 因此死亡细胞中线粒体mRNA的相对比例升高

典型阈值：

• 通常使用5-20%作为过滤阈值
• 具体阈值取决于组织类型（某些组织如心脏、肌肉天然具有高线粒体基因比例）

注意： 高线粒体基因比例不一定意味着细胞质量差。某些细胞类型（如心肌细胞、运动神经元）由于高能量需求，天然具有高线粒体基因表达。因此，线粒体基因比例应与其他QC指标联合使用。

4.2.4 核糖体基因比例（pct_counts_ribo）

定义： 来自核糖体蛋白基因的UMI计数占细胞总UMI计数的比例。

核糖体基因识别： 基因名以”RPS”（核糖体蛋白小亚基）或”RPL”（核糖体蛋白大亚基）开头。

生物学意义：

• 核糖体基因是管家基因，在所有活细胞中高度表达
• 通常占总counts的15-40%
• 极端高值可能指示细胞应激或技术偏差
• 可作为辅助QC指标

4.2.5 血红蛋白基因比例（pct_counts_hb）

定义： 来自血红蛋白基因的UMI计数占比。

血红蛋白基因识别： 基因名匹配^HB[^(P)]模式（如HBA1, HBA2, HBB等，排除HBP基因）

适用场景： 主要用于血液相关样本（如外周血单核细胞PBMC），检测红细胞污染。

4.3 MAD（中位数绝对偏差）过滤方法

4.3.1 原理

MAD（Median Absolute Deviation，中位数绝对偏差）是一种基于数据分布的鲁棒统计量，用于衡量数据的离散程度。与标准差相比，MAD对异常值（outlier）更加鲁棒，因为它是基于中位数而非均值计算的。

4.3.2 数学公式

MAD的定义：

$$\text{MAD} = \text{median}(|x_i - \text{median}(x)|)$$

即：先计算所有数据点与中位数的绝对偏差，再取这些绝对偏差的中位数。

基于MAD的异常值检测：

对于下限过滤（如nUMIs过低）：

$$\text{lower bound} = \text{median}(x) - n \times \text{MAD}(x)$$

对于上限过滤（如线粒体比例过高）：

$$\text{upper bound} = \text{median}(x) + n \times \text{MAD}(x)$$

其中$n$是MAD的倍数，通常取3或5。$n$值越大，过滤越宽松。

MAD与标准差的关系： 对于正态分布，MAD约等于0.6745倍的标准差：

$$\sigma \approx 1.4826 \times \text{MAD}$$

因此，3个MAD约等于2个标准差，5个MAD约等于3.4个标准差。

4.3.3 MAD过滤的优势

4.3.4 MAD过滤代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

import numpy as np
import pandas as pd

def mad_filter(adata, metrics, nmads=5, direction='both'):
    """
    基于MAD的自动质量控制过滤
    
    参数：
        adata: AnnData对象
        metrics: dict, 格式为 {指标名: 过滤方向}
            'both': 双侧过滤
            'upper': 仅上限过滤
            'lower': 仅下限过滤
        nmads: float, MAD的倍数（默认5）
        direction: 保留未使用的向后兼容参数
    
    返回：
        过滤后的AnnData对象
    """
    mask = pd.Series(True, index=adata.obs_names)
    
    for metric, filt_dir in metrics.items():
        values = adata.obs[metric]
        median = np.median(values)
        mad = np.median(np.abs(values - median))
        
        lower = median - nmads * mad
        upper = median + nmads * mad
        
        if filt_dir == 'lower' or filt_dir == 'both':
            mask = mask & (values >= lower)
        if filt_dir == 'upper' or filt_dir == 'both':
            mask = mask & (values <= upper)
        
        print(f"{metric}: median={median:.2f}, MAD={mad:.2f}, "
              f"range=[{lower:.2f}, {upper:.2f}], "
              f"filtered={(~(values >= lower) | ~(values <= upper)).sum()}")
    
    return adata[mask].copy()

# 使用MAD过滤
metrics = {
    'total_counts': 'lower',       # 低nUMIs过滤（去除低质量细胞）
    'n_genes_by_counts': 'lower',   # 低基因数过滤
    'pct_counts_mt': 'upper',       # 高线粒体比例过滤
}

adata_filtered = mad_filter(adata, metrics, nmads=5)

4.4 Seurat式手动阈值 vs MAD自动阈值对比

Seurat式手动阈值代码：

1
2
3
4
5
6
7

# Seurat风格手动阈值过滤
adata = adata[
    (adata.obs['total_counts'] > 500) &
    (adata.obs['n_genes_by_counts'] > 250) &
    (adata.obs['pct_counts_mt'] < 15),
    :
].copy()

4.5 SoupX环境RNA校正

4.5.1 原理

在液滴式scRNA-seq中，细胞被裂解后释放的mRNA会溶解在缓冲液中，形成所谓的”汤”（Soup）。这些游离的mRNA（环境RNA，ambient RNA）可能被错误地分配到同一液滴中的细胞，导致背景污染。

SoupX的工作原理：

1. 使用raw_feature_bc_matrix（未过滤的矩阵，包含所有barcode）估计环境RNA的组成
2. 假设空液滴（不含细胞的barcode）中的RNA完全来自环境RNA
3. 通过低表达空液滴的基因表达谱构建环境RNA的轮廓（soup profile）
4. 估计每个细胞中环境RNA的污染比例（soup fraction，$\rho$）
5. 对每个基因在每个细胞中的计数进行校正：

$$x_{gc}^{\text{corrected}} = x_{gc} - \rho_c \cdot s_g \cdot N_c$$

其中$\rho_c$是细胞$c$的污染比例，$s_g$是基因$g$在环境RNA中的比例，$N_c$是细胞$c$的总counts。

4.5.2 SoupX代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# SoupX通常在R中运行，可通过rpy2在Python中调用
# 以下是R代码示例

# 安装SoupX
# install.packages("SoupX")

library(SoupX)

# 方法1：从10x Cell Ranger输出加载
sc <- load10X(data_dir = "path/to/cellranger_output/")

# 方法2：手动加载
sc <- SoupX::load10X(
    data_dir = "path/to/cellranger_output/",
    keep_droplets = FALSE  # TRUE则保留空液滴
)

# 自动估计污染比例
sc <- autoEstCont(sc)

# 手动指定污染比例
# sc <- setContFraction(sc, 0.05)  # 设置5%的污染比例

# 校正计数矩阵
sc <- adjustCounts(sc)

# 导出校正后的矩阵
counts_corrected <- getCounts(sc)

# 保存为10x格式
save10X(sc, "corrected_output/")

在Python中使用SoupX（通过rpy2）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import rpy2.robjects as ro
from rpy2.robjects import numpy2ri, pandas2ri
from rpy2.robjects.conversion import localconverter

# 激活转换器
numpy2ri.activate()
pandas2ri.activate()

# 在R中运行SoupX
ro.r('''
library(SoupX)

# 加载数据
sc <- load10X("path/to/cellranger_output/")

# 自动估计污染比例
sc <- autoEstCont(sc)

# 校正
sc <- adjustCounts(sc)

# 获取校正后的计数矩阵
counts_corrected <- getCounts(sc)
''')

# 将结果转回Python
counts_corrected = ro.r('counts_corrected')

4.5.3 注意事项

1. 需要raw_feature_bc_matrix：SoupX需要未过滤的计数矩阵来估计环境RNA，但该文件通常不公开提供
2. 可能产生非整数计数：校正后的矩阵可能包含小数值，某些下游工具可能不支持
3. 仅在污染严重时使用：如果环境RNA污染不严重，不建议使用SoupX，因为校正可能引入新的偏差
4. 污染比例估计：自动估计不一定准确，建议手动检查并调整

4.6 双细胞检测

4.6.1 Scrublet原理与代码

原理： Scrublet（Single-cell Remover of Doublets）由Wolock等人于2019年开发，是一种基于模拟双细胞的Python双细胞检测工具。

核心思路：

1. 模拟双细胞： 从实际数据中随机选择两个细胞的表达谱，将它们相加，创建模拟双细胞（synthetic doublet）
2. 构建KNN图： 在PCA空间中，对真实细胞和模拟双细胞构建KNN图
3. 计算双细胞评分： 对每个真实细胞，计算其K个最近邻中模拟双细胞的比例，作为双细胞评分（doublet score）
4. 阈值判定： 基于双细胞评分的分布，设定阈值将细胞分类为singlet或doublet

数学细节：

对于细胞$c$，其双细胞评分定义为：

$$\text{score}_c = \frac{\text{number of synthetic doublets among } k \text{ nearest neighbors of } c}{k}$$

预期双细胞率：

$$\text{expected doublet rate} \approx 0.01 \times \frac{N_{\text{cells}}}{1000}$$

代码示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

import scrublet as scr
import scipy.io

# 加载计数矩阵
counts_matrix = scipy.io.mmread('counts.mtx').T.tocsr()

# 初始化Scrublet
scrub = scr.Scrublet(counts_matrix, expected_doublet_rate=0.06)

# 运行双细胞检测
doublet_scores, predicted_doublets = scrub.scrub_doublets(
    min_counts=2,
    min_cells=3,
    min_gene_variability_pctl=85,
    n_prin_comps=30,
)

# 可视化双细胞评分分布
scrub.plot_histogram()

# 将结果添加到AnnData
adata.obs['doublet_score'] = doublet_scores
adata.obs['predicted_doublet'] = predicted_doublets

# 过滤双细胞
adata = adata[~adata.obs['predicted_doublet']].copy()
print(f"Removed {(predicted_doublets).sum()} doublets")
print(f"Remaining cells: {adata.n_obs}")

参数说明：

4.6.2 scDblFinder原理与代码

原理： scDblFinder由Germain等人于2022年开发，是基于R/Bioconductor的双细胞检测工具。与Scrublet不同，scDblFinder使用机器学习分类器（梯度提升树）而非KNN方法。

核心思路：

1. 模拟双细胞： 类似Scrublet，随机选择两个细胞创建模拟双细胞
2. 特征工程： 为每个细胞（真实+模拟）计算一系列特征：
   • KNN中模拟双细胞的比例
   • 表达矩阵的PCs
   • 聚类信息
   • nUMIs、n_genes等QC指标
3. 训练分类器： 使用模拟数据训练梯度提升分类器
4. 预测： 对真实细胞进行双细胞概率预测

优势（相比Scrublet）：

• 在基准测试中准确率更高
• 对异源双细胞检测能力更强
• 可同时报告双细胞评分和二元分类结果
• 可检测多细胞（multiplet）

代码示例（R）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# 安装scDblFinder
BiocManager::install("scDblFinder")

library(scDblFinder)
library(SingleCellExperiment)

# 创建SCE对象（或从文件读取）
sce <- SingleCellExperiment(
    assays = list(counts = counts_matrix)
)

# 运行scDblFinder
sce <- scDblFinder(sce)

# 查看结果
table(sce$scDblFinder.class)
# singlet doublet 
#   7500    500 

# 双细胞评分
head(sce$scDblFinder.score)

# 过滤双细胞
sce_filtered <- sce[, sce$scDblFinder.class == "singlet"]

在Python中调用scDblFinder（通过rpy2）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import anndata2ri
anndata2ri.activate()

# 将AnnData传递到R
from rpy2.robjects import globalenv
globalenv['adata'] = adata

# 在R中运行scDblFinder
%%R
library(scDblFinder)
library(SingleCellExperiment)

# 从AnnData创建SCE
sce <- SingleCellExperiment(assays = list(counts = t(adata$X)))

# 运行scDblFinder
sce <- scDblFinder(sce, dbr = 0.06)  # dbr: expected doublet rate

# 将结果传回Python
doublet_class = %R sce$scDblFinder.class
doublet_score = %R sce$scDblFinder.score

adata.obs['scDblFinder_class'] = doublet_class
adata.obs['scDblFinder_score'] = doublet_score

4.7 omicverse一站式QC代码

omicverse提供了ov.pp.qc函数，集成了双细胞检测和MAD过滤于一体，是一站式的质量控制工具。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

import omicverse as ov
import scanpy as sc

# 读取数据
adata = sc.read_10x_mtx('filtered_feature_bc_matrix/')

# 标记线粒体/核糖体/血红蛋白基因
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))
adata.var['hb'] = adata.var_names.str.contains('^HB[^(P)]')

# 计算QC指标
sc.pp.calculate_qc_metrics(
    adata,
    qc_vars=['mt', 'ribo', 'hb'],
    percent_top=None,
    log1p=False,
    inplace=True
)

# 一站式QC（Seurat模式：手动阈值）
adata = ov.pp.qc(
    adata,
    mode='seurat',
    tresh={
        'mito_perc': 0.05,      # 线粒体基因比例上限5%
        'nUMIs': 500,            # nUMIs下限500
        'detected_genes': 250,   # 检测基因数下限250
    },
    min_cells=3,                 # 基因至少在3个细胞中表达
    min_genes=200,               # 细胞至少表达200个基因
)

# 一站式QC（MAD模式：自动阈值）
adata = ov.pp.qc(
    adata,
    mode='mads',
    nmads=5,                     # 5倍MAD
    min_cells=3,
    min_genes=200,
)

# 保存QC后的数据
adata.write_h5ad('quality_control.h5ad')

4.8 完整QC流程代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

import omicverse as ov
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

ov.ov_plot_set()

# ============================================
# Step 1: 读取数据
# ============================================
adata = sc.read_10x_mtx('filtered_feature_bc_matrix/', var_names='gene_symbols')
print(f"原始数据: {adata.shape}")

# ============================================
# Step 2: 标记特殊基因
# ============================================
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))
adata.var['hb'] = adata.var_names.str.contains('^HB[^(P)]')

# ============================================
# Step 3: 计算QC指标
# ============================================
sc.pp.calculate_qc_metrics(
    adata,
    qc_vars=['mt', 'ribo', 'hb'],
    percent_top=None,
    log1p=False,
    inplace=True
)

print(f"QC指标计算完成:")
print(f"  total_counts: {adata.obs['total_counts'].median():.0f} (median)")
print(f"  n_genes_by_counts: {adata.obs['n_genes_by_counts'].median():.0f} (median)")
print(f"  pct_counts_mt: {adata.obs['pct_counts_mt'].median():.1f}% (median)")

# ============================================
# Step 4: 可视化QC指标分布
# ============================================
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

sns.histplot(adata.obs['total_counts'], bins=100, kde=False, ax=axes[0])
axes[0].set_title('Total Counts (nUMIs)')
axes[0].set_xlabel('Total Counts')

sns.histplot(adata.obs['n_genes_by_counts'], bins=100, kde=False, ax=axes[1])
axes[1].set_title('Detected Genes')

sns.histplot(adata.obs['pct_counts_mt'], bins=100, kde=False, ax=axes[2])
axes[2].set_title('Mitochondrial Gene %')

plt.tight_layout()
plt.savefig('qc_metrics_distribution.png', dpi=150)
plt.show()

# ============================================
# Step 5: 过滤低质量基因
# ============================================
# 基因至少在3个细胞中表达
sc.pp.filter_genes(adata, min_cells=3)
print(f"过滤低质量基因后: {adata.shape}")

# ============================================
# Step 6: 双细胞检测（Scrublet）
# ============================================
import scrublet as scr

scrub = scr.Scrublet(adata.X, expected_doublet_rate=0.06)
doublet_scores, predicted_doublets = scrub.scrub_doublets(
    min_counts=2,
    min_cells=3,
    min_gene_variability_pctl=85,
    n_prin_comps=30,
)

adata.obs['doublet_score'] = doublet_scores
adata.obs['predicted_doublet'] = predicted_doublets

print(f"检测到双细胞: {predicted_doublets.sum()} ({predicted_doublets.mean()*100:.1f}%)")

# 过滤双细胞
adata = adata[~adata.obs['predicted_doublet']].copy()
print(f"去除双细胞后: {adata.shape}")

# ============================================
# Step 7: 基于MAD的质量过滤
# ============================================
def compute_mad_threshold(values, nmads=5, direction='lower'):
    """计算MAD阈值"""
    median = np.median(values)
    mad = np.median(np.abs(values - median))
    if direction == 'lower':
        threshold = median - nmads * mad
    elif direction == 'upper':
        threshold = median + nmads * mad
    return threshold, median, mad

# 计算各指标的MAD阈值
nmads = 5

lower_nUMIs, med, mad = compute_mad_threshold(
    adata.obs['total_counts'], nmads, 'lower'
)
print(f"nUMIs: median={med:.0f}, MAD={mad:.0f}, threshold={lower_nUMIs:.0f}")

lower_genes, med, mad = compute_mad_threshold(
    adata.obs['n_genes_by_counts'], nmads, 'lower'
)
print(f"n_genes: median={med:.0f}, MAD={mad:.0f}, threshold={lower_genes:.0f}")

upper_mt, med, mad = compute_mad_threshold(
    adata.obs['pct_counts_mt'], nmads, 'upper'
)
print(f"pct_mt: median={med:.1f}%, MAD={mad:.1f}%, threshold={upper_mt:.1f}%")

# 应用过滤
mask = (
    (adata.obs['total_counts'] >= lower_nUMIs) &
    (adata.obs['n_genes_by_counts'] >= lower_genes) &
    (adata.obs['pct_counts_mt'] <= upper_mt)
)
adata = adata[mask].copy()
print(f"MAD过滤后: {adata.shape}")

# ============================================
# Step 8: 保存QC后的数据
# ============================================
adata.write_h5ad('quality_control.h5ad')
print(f"QC完成! 最终数据: {adata.shape}")

4.9 注意事项和常见陷阱

1. 不要过度过滤：过滤太严格可能丢失生物学有意义的细胞。例如，某些小型细胞（如血小板）天然具有低nUMIs和低基因数。
2. 线粒体基因比例的陷阱：高线粒体比例不一定意味着低质量细胞。某些细胞类型（心肌细胞、运动神经元、肾小管细胞）天然具有高线粒体含量。建议结合其他指标联合判断。
3. MAD倍数选择：nmads=5是较宽松的阈值（约3.4σ），nmads=3更严格（约2σ）。建议先用5，再根据数据特点调整。
4. 双细胞检测的时机：应在归一化之前进行双细胞检测，因为归一化可能改变双细胞的特征。
5. SoupX的使用条件：需要raw_feature_bc_matrix文件，且仅在环境RNA污染严重时使用。轻度污染不建议使用。
6. QC指标的相关性：nUMIs和n_genes高度相关，不要对两者都进行严格的上下限过滤，否则可能过度去除细胞。
7. 批次差异：不同批次的QC指标分布可能不同，建议每个批次独立设定阈值，或使用MAD方法自动适应。
8. 细胞类型特异性：过滤前考虑是否了解数据中可能存在的细胞类型及其预期QC特征。
9. 空液滴检测：对于10x数据，Cell Ranger已经进行了细胞呼叫，但有时会遗漏一些真实细胞。可以检查raw_feature_bc_matrix中的低UMI barcode是否包含被遗漏的细胞。
10. 过滤后重新计算QC指标：如果对数据进行了大幅过滤，建议重新计算QC指标并检查分布是否合理。

5. 归一化

5.1 为什么需要归一化

5.1.1 技术变异的来源

单细胞RNA测序数据中，细胞间的计数差异并不完全反映真实的生物学差异。由于技术原因，即使同一类型的细胞，其测序计数深度也可能存在巨大差异。这些技术变异的来源包括：

1. 捕获效率差异：不同细胞的mRNA被捕获的效率不同
2. 逆转录效率差异：不同细胞的mRNA被逆转录为cDNA的效率不同
3. PCR扩增偏差：不同的cDNA分子被扩增的倍数不同
4. 测序深度差异：不同细胞的reads数量不同

5.1.2 归一化的目标

归一化的核心目标是：消除技术变异（特别是细胞间测序深度的差异），使得细胞间的基因表达差异反映真实的生物学差异，而非技术偏差。

数学上，归一化旨在将”UMI的方差”缩放到指定范围，使得数据满足下游统计模型的假设（如同方差性）。

5.2 伽玛-泊松分布（Gamma-Poisson Distribution）

5.2.1 原理

从理论上和经验上建立的UMI数据模型是伽玛-泊松分布（Gamma-Poisson Distribution），又称负二项分布（Negative Binomial Distribution）。

对于基因$g$在细胞$c$中的UMI计数$Y_{gc}$，假设：

$$Y_{gc} | \lambda_{gc} \sim \text{Poisson}(\lambda_{gc})$$

$$\lambda_{gc} \sim \text{Gamma}(\mu_{gc}/\theta, 1/\theta)$$

合并后，$Y_{gc}$服从负二项分布：

$$Y_{gc} \sim \text{NB}(\mu_{gc}, \theta)$$

其方差为：

$$\text{Var}[Y] = \mu + \alpha \mu^2$$

其中：

• $\mu$代表UMI平均值（期望表达量）
• $\alpha = 1/\theta$是过度离散参数（overdispersion parameter）
• $\theta$是形状参数（shape parameter）

5.2.2 参数解释

$\theta$（形状参数）：

• $\theta$越大，过度离散越小，数据越接近泊松分布
• $\theta$越小，过度离散越大，数据方差远大于均值

$\alpha$（过度离散参数）：

• $\alpha = 0$时，$\text{Var}[Y] = \mu$，数据服从泊松分布
• $\alpha > 0$时，数据存在过度离散（overdispersion）
• 单细胞数据通常具有显著的过度离散（$\alpha > 0$）

归一化与伽玛-泊松分布的关系： 不同的归一化方法本质上是对伽玛-泊松分布中不同参数的处理策略，目标是将数据转换到方差稳定的空间。

5.3 移位对数（Shifted Logarithm）

5.3.1 原理

移位对数是最常用的归一化方法，基于delta方法（Delta Method）实现方差稳定化。

Delta方法： 对于非线性函数$f(Y)$，如果$Y$的方差为$\sigma^2$，则$f(Y)$的方差近似为：

$$\text{Var}[f(Y)] \approx [f’(E[Y])]^2 \cdot \text{Var}[Y]$$

对于伽玛-泊松分布，$\text{Var}[Y] = \mu + \alpha\mu^2$。如果我们选择$f(y) = \log(y)$，则：

$$\text{Var}[\log(Y)] \approx \frac{\mu + \alpha\mu^2}{\mu^2} = \frac{1}{\mu} + \alpha$$

当$\mu$足够大时，$\text{Var}[\log(Y)] \approx \alpha$，方差近似恒定。这就是对数变换稳定方差的数学基础。

5.3.2 公式

移位对数的变换公式：

$$f(y) = \log\left(\frac{y}{s} + y_0\right)$$

其中：

• $y$是原始UMI计数
• $s$是尺寸因子（size factor）
• $y_0$是伪计数（pseudo-count），通常为1

尺寸因子$s_c$的计算：

$$s_c = \frac{\sum_g y_{gc}}{L}$$

其中：

• $g$遍历所有基因
• $L$是目标计数深度（target sum）

不同工具的$L$值选择：

• Scanpy默认：$L = \text{median}(\sum_g y_{gc})$，即所有细胞总counts的中位数
• Seurat默认：$L = 10^4$
• omicverse默认：$L = 50 \times 10^4$

不同$L$值会影响过度离散值$\alpha$，但不影响下游分析结果的本质。

5.3.3 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

import scanpy as sc
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

# 方法1：使用scanpy的标准流程
# 先归一化到中位计数深度
scales_counts = sc.pp.normalize_total(adata, target_sum=None, inplace=False)

# 再进行log1p变换
adata.layers["log1p_norm"] = sc.pp.log1p(scales_counts["X"], copy=True)

# 方法2：直接修改adata.X（推荐用于下游分析）
sc.pp.normalize_total(adata, target_sum=None)  # 默认使用中位数
sc.pp.log1p(adata)

# 可视化归一化前后的分布
fig, axes = plt.subplots(1, 2, figsize=(10, 5))

p1 = sns.histplot(adata.obs["total_counts"], bins=100, kde=False, ax=axes[0])
axes[0].set_title("Total counts (before normalization)")

p2 = sns.histplot(adata.layers["log1p_norm"].sum(1), bins=100, kde=False, ax=axes[1])
axes[1].set_title("Shifted logarithm (after normalization)")

plt.tight_layout()
plt.show()

5.3.4 适用场景

• ✅ 大多数scRNA-seq数据分析任务
• ✅ 降维（PCA）和聚类
• ✅ 差异表达基因识别
• ✅ 2023年Nature Methods基准测试中表现最佳的方法之一
• ⚠️ 对稀有细胞类型的检测不如Pearson残差

5.4 scran反卷积归一化

5.4.1 原理

scran归一化由Lun等人于2016年开发，其核心思想是通过**反卷积（Deconvolution）**方法来估计更准确的尺寸因子。

为什么需要更好的尺寸因子：
简单地将各细胞的总counts归一化到相同深度，假设了所有细胞的RNA含量相同。但实际上，不同细胞类型可能有不同的RNA含量（如大细胞比小细胞有更多RNA）。

scran的方法：

1. 对细胞进行初步聚类（粗粒度），将相似细胞分组
2. 在每组内，将多个细胞的counts求和（pooling），形成”伪细胞”
3. 对伪细胞计算尺寸因子，减少零膨胀的影响
4. 通过线性回归从伪细胞的尺寸因子中解卷积出单个细胞的尺寸因子

数学公式：

假设将$P$个细胞池化，池化后的counts为：

$$Y_{gp}^{\text{pool}} = \sum_{c \in p} Y_{gc}$$

池化后的尺寸因子为：

$$s_p^{\text{pool}} = \sum_{c \in p} s_c$$

通过构建线性方程组，从池化尺寸因子中求解单个细胞的尺寸因子$s_c$。

5.4.2 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

import scanpy as sc
import numpy as np
from scipy.sparse import csr_matrix

# scran需要R环境，以下是通过rpy2调用的完整代码

# Step 1: 准备数据 - 初步聚类
adata_pp = adata.copy()
sc.pp.normalize_total(adata_pp)
sc.pp.log1p(adata_pp)
sc.pp.pca(adata_pp, n_comps=50)
sc.pp.neighbors(adata_pp, n_neighbors=10, n_pcs=50)
sc.tl.leiden(adata_pp, resolution=0.5, key_added='groups')

# Step 2: 将数据传递到R环境
import rpy2.robjects as ro
from rpy2.robjects import numpy2ri, pandas2ri
from rpy2.robjects.conversion import localconverter

data_mat = adata.X.T
if hasattr(data_mat, 'toarray'):
    data_mat = data_mat.tocsc().toarray()

with localconverter(ro.default_converter + numpy2ri.converter):
    ro.globalenv["data_mat"] = data_mat

with localconverter(ro.default_converter + pandas2ri.converter):
    ro.globalenv["input_groups"] = adata_pp.obs["groups"]

# Step 3: 在R中计算尺寸因子
ro.r('''
library(scran)
library(BiocParallel)

size_factors <- sizeFactors(
    computeSumFactors(
        SingleCellExperiment(list(counts=data_mat)),
        clusters = input_groups,
        min.mean = 0.1,
        BPPARAM = MulticoreParam()
    )
)
''')

# Step 4: 获取尺寸因子并应用归一化
size_factors = ro.r('size_factors')
adata.obs["size_factors"] = size_factors

# 手动应用scran归一化
scran_normalized = adata.X / adata.obs["size_factors"].values[:, None]
scran_logged = np.log1p(scran_normalized)
adata.layers["scran_normalization"] = csr_matrix(scran_logged)

5.4.3 适用场景

• ✅ 批次校正任务（scran归一化考虑了细胞类型特异性的尺寸因子）
• ✅ 需要更准确尺寸因子的分析
• ⚠️ 需要R环境
• ⚠️ 计算速度较慢
• ⚠️ 对初步聚类的结果有一定依赖

5.5 Pearson残差归一化

5.5.1 原理

Pearson残差归一化（Analytic Pearson Residuals）由Lause等人于2021年提出，基于正则化负二项回归（Regularized Negative Binomial Regression）。

核心思想：
不再通过尺寸因子和log变换来归一化，而是直接从统计模型中提取残差，将技术变异和生物变异分离。

数学推导：

对于基因$g$在细胞$c$中的计数$Y_{gc}$，假设其服从负二项分布：

$$Y_{gc} \sim \text{NB}(\mu_{gc}, \theta_g)$$

其中均值$\mu_{gc}$由细胞$c$的测序深度决定：

$$\log(\mu_{gc}) = \beta_{g0} + \beta_{g1} \log(s_c)$$

• $\beta_{g0}$：基因$g$的截距（baseline表达）
• $\beta_{g1}$：测序深度效应（通常$\beta_{g1} \approx 1$）
• $s_c$：细胞$c$的测序深度（total counts的某种度量）

Pearson残差定义：

$$Z_{gc} = \frac{Y_{gc} - \hat{\mu}{gc}}{\sqrt{\hat{\mu}{gc} + \hat{\mu}_{gc}^2/\theta_g}}$$

其中：

• $\hat{\mu}_{gc}$：模型预测的均值
• $\theta_g$：基因$g$的过度离散参数
• 分母是负二项分布的标准差

Pearson残差的性质：

• 正值：该基因在该细胞中的表达高于预期
• 负值：该基因在该细胞中的表达低于预期
• 接近零：该基因在该细胞中的表达与预期一致
• 残差的方差近似恒定（方差稳定化）

正则化： 为了避免对低表达基因的$\theta$估计不准确，使用全局的$\theta$-均值关系进行正则化（即从高表达基因估计$\theta$-均值关系，然后外推到低表达基因）。

5.5.2 代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import scanpy as sc
from scipy.sparse import csr_matrix
import seaborn as sns
import matplotlib.pyplot as plt

# 方法1：使用scanpy的experimental模块
analytic_pearson = sc.experimental.pp.normalize_pearson_residuals(adata, inplace=False)
adata.layers["analytic_pearson_residuals"] = csr_matrix(analytic_pearson["X"])

# 方法2：使用omicverse
import omicverse as ov
adata_pearson = ov.pp.preprocess(adata, mode='pearson|pearson', n_HVGs=2000)

# 可视化Pearson残差的分布
fig, axes = plt.subplots(1, 2, figsize=(10, 5))

sns.histplot(adata.obs["total_counts"], bins=100, kde=False, ax=axes[0])
axes[0].set_title("Total counts")

sns.histplot(adata.layers["analytic_pearson_residuals"].sum(1), bins=100, kde=False, ax=axes[1])
axes[1].set_title("Analytic Pearson residuals")

plt.tight_layout()
plt.show()

5.5.3 适用场景

• ✅ 识别稀有细胞类型（Pearson残差保留了生物学异质性）
• ✅ 选择高可变基因（直接基于残差方差选择HVG）
• ✅ 不需要伪计数和log变换
• ⚠️ 输出值可正可负，某些工具可能不兼容
• ⚠️ 不适合直接用于差异表达分析（差异表达通常需要计数矩阵）

5.6 三种方法完整对比