Basecalling 算法

快速概览

Basecalling 是将 Nanopore 原始离子电流信号解码为 DNA/RNA 碱基序列的核心算法。它通过建模信号与 k-mer 的对应关系，实现纳米孔测序的「读出」功能。

• Nanopore 测量 k-mer 穿过纳米孔时的离子电流变化，而非直接识别碱基
• HMM 通过状态转移和发射概率建模信号生成过程
• 神经网络（CNN/RNN/Transformer）直接学习信号到碱基的端到端映射
• 现代 basecaller 已超越传统 HMM，采用深度学习方法

前置知识

所属板块 分析方向与案例

把基础对象与算法方法重新放回真实分析任务与工作流。

阅读目标 帮助建立阅读上下文

先判断这页与你当前问题的关系，再决定是否深入展开。

建议前置 先建立相关基础对象与方法直觉

建议先建立相关基础对象与方法直觉，再进入本页。

问题引入：如何”阅读”电流信号

Nanopore 测序的物理原理与 Illumina 或 PacBio 截然不同：

Illumina：检测荧光标记的核苷酸 PacBio：检测荧光标记的核苷酸掺入事件 Nanopore：检测 DNA 穿过纳米孔时的离子电流变化

核心问题是：如何将连续的电流信号时间序列转换为离散的碱基序列？

这类似于语音识别问题——将连续的声波信号转换为离散的音素/文字序列。Basecalling 算法就是解决这一解码问题的核心。

物理基础：信号从何而来

当单链 DNA 在电场驱动下穿过蛋白质纳米孔时，孔道中的 k-mer（k 个连续碱基）会阻挡离子流，产生特征性的电流信号。

信号生成模型

电流信号可建模为：

$$I(t) = f(k(t)) + \epsilon(t)$$

其中：

• $I(t)$：时刻 $t$ 的观测电流
• $k(t)$：时刻 $t$ 位于孔内的 k-mer（通常 $k=5$）
• $f(\cdot)$：k-mer 到电流的映射函数（非线性）
• $\epsilon(t)$：噪声（热噪声、电子噪声、分子动力学波动）

关键观察：

• 不同 k-mer 产生不同的特征电流（典型范围 20-200 pA）
• 同一 k-mer 由于噪声产生电流分布（近似高斯）
• 相邻信号高度重叠（一个碱基移动只改变 k-mer 的一个位置）

数学框架：解码问题形式化

给定观测信号序列 $S = s_1, s_2, \ldots, s_T$，目标是找到最可能的碱基序列 $B = b_1, b_2, \ldots, b_L$：

$$B^* = \arg\max_B P(B | S) = \arg\max_B P(S | B) \cdot P(B)$$

其中：

• $P(S | B)$：发射概率（信号生成模型）
• $P(B)$：序列先验（如 k-mer 频率、上下文约束）

由于信号与碱基的比例关系不确定（一个碱基移动可能产生多个信号点），这是一个序列到序列的推断问题。

HMM 模型：概率图方法

模型结构

隐马尔可夫模型（HMM）将 basecalling 建模为：

• 隐状态：当前位于孔内的 k-mer（共 $4^k$ 个可能状态）
• 观测值：电流信号 $s_t$
• 状态转移：由于 DNA 是连续移动的，从 k-mer $x_1x_2\ldots x_k$ 只能转移到 $x_2\ldots x_k y$（其中 $y \in \{A,C,G,T\}$）
• 发射概率：$P(s_t | \text{k-mer})$，通常建模为高斯分布

发射概率建模

对于每个 k-mer $k_i$，观测电流服从：

$$P(s | k_i) = \mathcal{N}(s; \mu_i, \sigma_i^2) = \frac{1}{\sqrt{2\pi\sigma_i^2}} \exp\left(-\frac{(s - \mu_i)^2}{2\sigma_i^2}\right)$$

参数 $(\mu_i, \sigma_i)$ 通过训练数据估计。

Viterbi 解码

使用 Viterbi 算法找到最优状态序列：

$$\delta_t(i) = \max_j \left[ \delta_{t-1}(j) \cdot a_{ji} \right] \cdot b_i(s_t)$$

其中 $\delta_t(i)$ 是时刻 $t$ 处于状态 $i$ 的最优路径概率。

复杂度：$O(T \times 4^k \times 4) = O(T \cdot 4^k)$。对于 $k=5$（1024 状态），这在计算上是可行的。

HMM 的局限性

1. 特征独立性假设：发射概率仅依赖当前 k-mer，忽略上下文
2. 固定 k-mer 大小：难以捕捉长程依赖
3. 高斯噪声假设：实际噪声分布更复杂

这些局限性推动了神经网络方法的发展。

神经网络模型：端到端学习

架构演进

现代 basecaller 采用深度神经网络，架构经历了三代演进：

第一代：CNN（如 DeepNano）

Signal → 1D Conv → ReLU → Pool → ... → FC → Softmax

• 局部特征提取强，但长程依赖建模弱

第二代：RNN/LSTM（如 Guppy 早期版本）

Signal → LSTM → LSTM → ... → FC → Softmax

• 自然处理变长序列，建模时间依赖
• 训练较慢，并行化困难

第三代：Transformer（如 Dorado）

Signal → Embedding → Multi-Head Attention → ... → Softmax

• 全局上下文建模，高度并行化
• 计算资源需求高，但准确性最佳

训练与损失函数

使用大量（信号，参考序列）配对数据训练，优化交叉熵损失：

$$L = -\sum_t \sum_{b \in \{A,C,G,T\}} y_{t,b} \log(\hat{p}_{t,b})$$

关键挑战是建立信号与参考的对齐（信号与碱基的比例关系不确定）。

关键技术与挑战

信号预处理

归一化：消除测序批次和温度差异

• Z-score 归一化：$s'_t = (s_t - \mu) / \sigma$
• 分位数归一化：将信号映射到固定范围

上下文建模：

• 滑动窗口：处理信号与碱基的非确定性对应关系
• 注意力机制：在 Transformer 中聚合全局上下文

主要挑战

同聚物问题： 同一碱基连续重复（如 AAAA）时，信号变化缓慢，难以确定重复次数。这是 Nanopore basecalling 的主要误差来源。

速度波动： DNA 穿过纳米孔的速度不均匀，导致信号时间扭曲。

修饰碱基： 甲基化等修饰改变电流信号，需要专门训练模型或将其作为额外输出维度。

示例

简化 HMM 示例

假设 k=2，只有 4 个状态：AA, AC, AG, AT（简化示例）

信号序列：[10.5, 12.3, 11.8, 13.2]

发射分布（均值）：

• AA: μ=10.0, σ=0.5
• AC: μ=11.0, σ=0.5
• AG: μ=12.0, σ=0.5
• AT: μ=13.0, σ=0.5

转移概率：均匀（每个状态到下一个状态概率相等）

步骤 1：计算发射概率

对于信号 10.5：

• P(10.5|AA) ≈ exp(-(10.5-10)²/2) ≈ 0.88
• P(10.5|AC) ≈ exp(-(10.5-11)²/2) ≈ 0.88
• P(10.5|AG) ≈ exp(-(10.5-12)²/2) ≈ 0.32
• P(10.5|AT) ≈ exp(-(10.5-13)²/2) ≈ 0.04

类似计算其他信号点。

步骤 2：Viterbi 递推

初始化：

• δ₁(AA) = 0.88
• δ₁(AC) = 0.88
• δ₁(AG) = 0.32
• δ₁(AT) = 0.04

递推（简化）：

• δ₂(AC) = max(δ₁·A) · P(12.3|AC) ≈ 0.88 · 0.04 = 0.04
• δ₂(AG) = max(δ₁·A) · P(12.3|AG) ≈ 0.88 · 0.88 = 0.77
• …

步骤 3：回溯

得到最优状态序列，例如：AA → AG → AG → AT

步骤 4：转换为碱基

取每个 k-mer 的第一个碱基：A → A → A → A

最终序列：AAAA

复杂度分析

HMM

时间复杂度

• Viterbi 算法：O(T · N²)，T 为信号长度，N 为状态数（4ᵏ）
• 由于转移约束，实际为 O(T · N · 4) = O(T · N)

总时间复杂度：O(T · 4ᵏ)

空间复杂度

• DP 表：O(T · N)

总空间复杂度：O(T · 4ᵏ)

神经网络

时间复杂度

• CNN：O(T · d)，d 为模型复杂度（参数量 × 计算量）
• RNN：O(T · d)
• Transformer：O(T² · d)（自注意力）

总时间复杂度：O(T · d) 或 O(T² · d)

空间复杂度

• 模型参数：O(d)
• 中间激活：O(T · h)，h 为隐藏层大小

总空间复杂度：O(d + T · h)

与真实工具或流程的连接

Guppy

Oxford Nanopore 的官方 basecaller：

• 多种模型：fast, high-accuracy, super-accuracy
• 实时模式：支持 streaming basecalling
• GPU 加速：利用 GPU 加速推理
• 多种化学：适配不同测序化学版本

Bonito

开源的神经网络 basecaller：

• 模块化设计：易于自定义和扩展
• 支持多种架构：CNN, LSTM, Transformer
• 可训练：用户可以训练自己的模型

Dorado

最新的 basecaller，基于 Transformer：

• 更高准确性：利用 Transformer 的强大建模能力
• GPU 优化：高度优化的 GPU 实现
• 修饰检测：集成碱基修饰检测

常见误区

• Basecalling 只需要一次：
• 不是。通常需要多次迭代，特别是对于低质量数据。
• 信号归一化不重要：
• 重要。归一化显著影响模型性能，必须正确处理。
• HMM 已经过时：
• 不是。HMM 仍用于某些场景，特别是作为神经网络的后处理步骤。

算法优化

硬件加速

• GPU 加速：神经网络推理在 GPU 上可加速 10-100 倍
• 专用芯片：Google Coral、Edge TPU 等边缘设备
• 量化：降低数值精度（FP16, INT8）减少计算量

算法优化

• 模型剪枝：移除不重要的神经元
• 知识蒸馏：用大模型训练小模型
• 早停机制：实时 basecalling 中提前终止

流式处理

• 滑动窗口缓存：避免存储整个信号
• 增量解码：边接收信号边 basecalling
• 自适应采样：根据信号质量调整采样率

参考资料

• Timp, W., et al. (2016). Nanopore sequencing: electrical measurements of the code of life. IEEE Transactions on Nanobioscience, 15(3), 260-267.
• Boža, T. T., et al. (2020). Nanopore signal processing with hidden Markov models. Bioinformatics, 36(7), 2105-2111.
• Wick, R. R., et al. (2019). Integrating mapping, assembly and haplotype variant estimation of single-molecule sequencing data using the Medaka genome analysis toolkit. BioRxiv.

PacBio 与 Nanopore

长读长测序技术平台对比

延伸阅读

Consensus 算法

Basecalling 后的序列纠错与共识

延伸阅读

结构变异检测

长读长在 SV 检测中的应用

延伸阅读

序列深度学习

深度学习在序列分析中的更多应用

延伸阅读