Pseudo-alignment 与表达定量

快速概览

当 RNA-seq 的分析目标是表达定量时，算法策略从精确比对转向兼容性判断与概率丰度估计。本节介绍基于 k-mer 索引的 lightweight mapping 算法及其复杂度优势。

• 问题核心：read 归属的不确定性（multi-mapping）与计算效率的权衡
• 算法策略：k-mer 索引 + 等价类（equivalence class）+ EM 推断
• 复杂度分析：从 O(reads × transcript_length) 降至 O(reads × k-mer_matches)

前置知识

• RNA-seq 工作流概览
• 序列比对基础
• TPM、FPKM、CPM 与有效长度

问题定义

RNA-seq 表达定量问题：给定一组测序 reads 和参考转录组序列，估计每个转录本（或基因）的表达丰度。

• 输入：

  • 测序 reads 集合 $R = \{r_1, r_2, ..., r_n\}$，其中 $n$ 通常为 $10^6$ 到 $10^8$ 量级
  • 参考转录组序列集合 $T = \{t_1, t_2, ..., t_m\}$，包含 $m$ 个已知转录本

• 输出：

  • 每个转录本 $t_j$ 的丰度估计值 $\hat{\theta}_j$（可以是分子计数或相对丰度）
  • 可选：read-to-transcript 归属概率矩阵

• 核心难点：

  • 由于序列相似性，一个 read 可能兼容（compatible with）多个转录本
  • 需要推断这种归属不确定性（multi-mapping uncertainty）下的最优丰度估计

从精确比对到兼容性判断

传统方法的时间复杂度

经典 spliced alignment 算法（如 STAR、HISAT2）需要：

1. 为每个 read 找到精确的碱基级比对位置
2. 处理剪接边界、错配、插入缺失
3. 计算复杂度：$O(|r| \times |t|)$ 对于每个 read-transcript 对

对于 $n$ 个 reads 和 $m$ 个转录本，总复杂度可达 $O(n \times m \times L)$，其中 $L$ 为序列长度。

Pseudo-alignment 的算法思路

关键观察：表达定量问题不需要精确的碱基级比对，只需要知道：

1. 该 read 可能来源于哪些转录本（兼容性集合）
2. 这些转录本的相对丰度如何解释观察到的 reads 分布

算法设计：

1. 索引构建：将转录组序列分解为 k-mer，建立 k-mer 到转录本列表的映射
2. 兼容性查询：对于每个 read，提取其 k-mer，查询索引，确定兼容转录本集合
3. 丰度估计：使用 EM 算法或变分推断，在多重归属不确定性下估计丰度

时间复杂度改善：

• 索引查询：$O(|r|)$ 对于每个 read（假设 k-mer 查询为 $O(1)$）
• 总复杂度：$O(n \times |r|)$，相比精确比对降低 1-2 个数量级

核心算法概念

等价类（Equivalence Classes）

定义：将具有相同兼容性集合的 reads 归入同一等价类。

形式化地，对于转录本子集 $S \subseteq T$，定义等价类：

$C_S = \{r \in R : \text{Compatible}(r) = S\}$

其中 $\text{Compatible}(r)$ 表示与 read $r$ 兼容的所有转录本集合。

算法优势：

• 不需要逐个处理每个 read
• 可以按等价类聚合计数，显著降低 EM 算法的计算负担
• 等价类数量通常远小于 reads 数量

丰度估计的统计模型

给定丰度参数 $\theta = (\theta_1, ..., \theta_m)$，read $r$ 来源于转录本 $t_j$ 的概率为：

$P(t_j | r, \theta) = \frac{\theta_j \cdot P(r | t_j)}{\sum_{k \in \text{Compatible}(r)} \theta_k \cdot P(r | t_k)}$

其中 $P(r|t_j)$ 可由序列相似性（如 k-mer 匹配率）估计。

EM 算法迭代：

1. E-step：给定当前丰度估计 $\theta^{(t)}$，计算每个 read 的期望归属
2. M-step：更新丰度参数 $\theta^{(t+1)}$ 以最大化期望似然

收敛后得到最大似然估计（MLE）或最大后验估计（MAP）。

与传统 Spliced Alignment 的算法比较

| 算法特性 | 传统比对（STAR/HISAT2） | Pseudo-alignment（Salmon/Kallisto） | |---------|------------------------|-----------------------------------| | 核心操作 | 精确序列比对，剪接处理 | k-mer 匹配，兼容性查询 | | 时间复杂度 | $O(n \cdot m \cdot L)$ | $O(n \cdot k)$，$k$ 为常数 | | 空间复杂度 | 需要存储完整比对结果 | 仅需存储等价类计数 | | 信息粒度 | 碱基级精确位置 | 转录本级归属概率 | | 适用分析 | 剪接分析、新转录本发现、变异检测 | 已知转录本定量、差异表达 | | 算法假设 | 需要完整的参考基因组 | 依赖已知转录本注释 |

算法实例

等价类构建示例

假设有三个转录本 $t_1, t_2, t_3$，五个 reads 的兼容性如下：

| Read | 兼容转录本 | |------|-----------| | $r_1$ | $\{t_1\}$ | | $r_2$ | $\{t_1, t_2\}$ | | $r_3$ | $\{t_1, t_2\}$ | | $r_4$ | $\{t_2, t_3\}$ | | $r_5$ | $\{t_3\}$ |

等价类划分：

• $C_{\{t_1\}} = \{r_1\}$，计数 $c_1 = 1$
• $C_{\{t_1, t_2\}} = \{r_2, r_3\}$，计数 $c_2 = 2$
• $C_{\{t_2, t_3\}} = \{r_4\}$，计数 $c_3 = 1$
• $C_{\{t_3\}} = \{r_5\}$，计数 $c_4 = 1$

EM 算法只需在这些等价类上迭代，而非逐个处理 5 个 reads。

与真实工具或流程的连接

这类方法常用于：

• 快速转录本层表达定量；
• 大样本批量处理；
• 把计算重点放在丰度估计而不是定位细节上。

算法复杂度与适用条件分析

维度	传统比对（STAR/HISAT2）	Pseudo-alignment (Salmon/Kallisto)
核心操作	寻找 read 在基因组上的最佳位置	判断 read 与哪个转录本兼容
比对精度	碱基级精确坐标	转录本级兼容性判断
速度	较慢（需要完整 DP 比对）	快 10-50 倍（k-mer 查询）
输出格式	BAM/SAM（完整比对信息）	定量矩阵（TPM/Counts）
下游差异表达	需额外定量步骤（featureCounts）	直接输出定量结果

时间复杂度总结

| 算法步骤 | 传统比对 | Pseudo-alignment | |---------|---------|-----------------| | 索引构建 | $O(\sum |t_j|)$，构建 FM-index 或后缀数组 | $O(\sum |t_j|)$，提取所有 k-mer | | 单 read 处理 | $O(|r| \times \log |T|)$，比对搜索 | $O(|r| - k + 1)$，k-mer 查询 | | 全样本处理 | $O(n \times |r| \times \log |T|)$ | $O(n \times |r|)$ | | EM 迭代 | 不涉及 | $O(|E| \times m)$，$|E|$ 为等价类数 |

实践中，pseudo-alignment 的总运行时间可比传统比对快 10-100 倍。

算法适用条件

选择 Pseudo-alignment 的条件：

• 分析目标为已知转录本的表达定量
• 参考注释完整且可信（如 GENCODE、Ensembl）
• 样本量大，计算效率是关键瓶颈
• 不需要碱基级变异信息或新转录本发现

必须选择传统 Spliced Alignment 的条件：

• 研究剪接模式或选择性剪接事件的精确边界
• 需要发现新转录本或融合基因
• 分析基因组变异（SNP、Indel）对表达的影响
• 数据质量较差，需要精确的错配/插入缺失处理

算法局限性与注意事项

依赖注释的完备性

Pseudo-alignment 算法假设参考转录组注释是完备的。如果：

• 样本中存在大量未注释的转录本
• 参考注释版本与实验物种/品系差异较大

则定量结果会有系统偏差，低估这些基因的真实表达。

Multi-mapping 的处理边界

虽然 EM 算法可以处理多重归属，但在极端情况下仍有限制：

• 高度重复序列：如核糖体 RNA、转座子，可能匹配数百个位点
• 同源基因家族：序列相似度 >95% 的旁系同源基因难以区分

此时需要额外的生物学假设或过滤策略。

常见误区

• Pseudo-alignment 精度不如传统比对，因此结果不可靠：
• Pseudo-alignment 在表达定量任务上的精度并不低于传统比对，甚至在某些场景下更优。因为传统比对在面对多映射 read 时也需要做出任意选择或丢弃，而 pseudo-alignment 的 EM 算法通过概率模型更合理地处理了归属不确定性。精度差异主要取决于参考注释的完备性，而非比对方法本身。
• Kallisto/Salmon 的速度优势来自"跳过了比对"：
• 这种说法不够准确。Pseudo-alignment 并非跳过比对，而是用一种更轻量的兼容性判断替代了碱基级精确比对。核心加速来自两方面：k-mer 索引的 O(1) 查询（而非 FM-index 的对数级搜索）和等价类聚合（将数百万 reads 压缩为少量等价类再运行 EM）。这是一种算法策略的转变，而非简单的"偷工减料"。
• EM 算法一定会收敛到全局最优：
• EM 算法只保证收敛到似然函数的局部极大值，而非全局最大值。在丰度参数的初始值选择不当时，EM 可能收敛到较差的解。实践中，Salmon 和 Kallisto 通过在线学习（逐 read 更新）和合理的初始化策略来缓解这一问题，但不能完全消除局部最优的风险。

RNA-seq 工作流概览

延伸阅读

TPM、FPKM、CPM 与有效长度

延伸阅读

序列比对

延伸阅读

基因层与转录本层表达定量

延伸阅读