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 表达定量问题：给定一组测序 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）需要：

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

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

Pseudo-alignment 的算法思路

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

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

算法设计：

索引构建：将转录组序列分解为 k-mer，建立 k-mer 到转录本列表的映射
兼容性查询：对于每个 read，提取其 k-mer，查询索引，确定兼容转录本集合
丰度估计：使用 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 算法迭代：

E-step：给定当前丰度估计 $\theta^{(t)}$ ，计算每个 read 的期望归属
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
单 read 处理	$O(	r
全样本处理	$O(n \times	r
EM 迭代	不涉及	$O(

实践中，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 更新）和合理的初始化策略来缓解这一问题，但不能完全消除局部最优的风险。