使用牛津纳米孔测序进行甲基化检测的信号处理与深度学习框架.

0. TL; DR

DeepMod2 是一个用于利用纳米孔测序的离子电流信号进行甲基化检测的综合深度学习框架。DeepMod2 实现了双向长短期记忆 (bidirectional long short-term memory, BiLSTM) 模型和 Transformer 模型，并且可以分析在 R9 和 R10 流通池上生成的 POD5 和 FAST5 信号文件。

DeepMod2 能够处理经过分型 (phased) 的 reads，从而推断出表观单倍型 (epihaplotypes) 或单倍型特异性的甲基化调用，为印记区域等研究提供了有力工具。作者证明了使用自适应采样进行的简化代表性测序与全基因组测序之间存在高度相关性 ($r=0.96$)，为大规模、经济高效的甲基化谱分析提供了一种可行的策略。

1. 背景介绍

DNA 甲基化，即在 DNA 分子的特定核苷酸上添加甲基基团的过程，是人类和其他物种基因组中一种重要的表观遗传变化。5-甲基胞嘧啶 (5mC) 是人类中最普遍的 DNA 甲基化形式，它在基因组稳定性、基因组印记、衰老和癌变等多种生物学过程中都扮演着至关重要的角色。

传统的甲基化检测方法，如甲基化微阵列 (methylation microarrays) 和基于亚硫酸氢盐转化 (bisulfite conversion) 的短读长测序，虽然为 5mC 的研究做出了巨大贡献，但它们也存在显著的局限性，例如，它们难以分析基因组中的重复区域和表观单倍型 (epihaplotypes)。

牛津纳米孔技术 (Oxford Nanopore Technologies, ONT) 测序等长读长技术为解决这些问题提供了新的途径。它通过检测 DNA 分子穿过纳米孔时产生的离子电流信号变化，来直接区分甲基化和未甲基化的胞嘧啶。近年来，已有一系列工具被开发出来用于此目的，如 Nanopolish, Megalodon, DeepSignal 等。

此外，ONT 测序平台上的自适应采样 (adaptive sampling) 技术，允许在测序过程中实时选择目标 DNA 分子。这使得对 CpG 岛或印记区域等特定区域进行富集测序成为可能，实现了类似于简化代表性亚硫酸氢盐测序 (RRBS) 的简化代表性甲基化测序 (RRMS)。

然而，近年来，公开和开源的甲基化检测工具的开发严重滞后于 ONT 测序技术的发展。例如，大多数现有工具仅支持即将被淘汰的 R9.4 流通池，而无法处理新一代 R10.4 流通池产生的、信号特征完全不同的数据。此外，它们也无法处理 ONT 用于替代 FAST5 的新型 POD5 文件格式。

为了解决这些问题，作者提出了 DeepMod2，一个全面的深度学习框架，用于从 ONT 测序中进行甲基化检测。DeepMod2 对其前身 DeepMod 进行了重大改进，旨在提供一个准确、高效且能兼容所有 ONT 流通池和数据格式的开源工具。

2. DeepMod2方法

DeepMod2 的工作流程是一个集成了信号处理、特征提取和深度学习模型推理的综合框架。它接收 POD5/FAST5 文件中的离子电流信号和 BAM 文件中的 read 序列作为输入，最终输出逐 read 和逐位点的甲基化预测结果。

2.1 特征提取 (Feature Extraction)

为了从原始的纳米孔信号中提取用于模型预测的特征，DeepMod2 需要建立碱基检出 (basecalled) 序列与原始信号之间的映射。这一步依赖于 Guppy 或 Dorado 等碱基检出软件 (basecaller) 生成的移动表 (move table)。特征提取过程如下：

  1. 信号归一化：首先，使用中位数绝对偏差 (median absolute deviation) 对原始信号进行归一化。
  2. CpG 位点定位：在每个 read 上，确定所有匹配 CG 序列基序的位点。如果提供了比对好的 BAM 文件和参考基因组，则还会确定那些在参考基因组上映射到 CG 基序的位点。
  3. 特征矩阵生成：对于每个待预测的 CpG 位点，以其为中心构建一个21 bp的窗口。对于这个窗口中的每一个碱基，提取以下19个特征，最终构成一个 $21 \times 19$ 的特征矩阵：
    • 碱基信号的长度（log10）
    • 碱基信号的均值、标准差、中位数、中位数绝对偏差
    • 碱基信号四个四分位的均值
    • 碱基检出 (basecall) 质量分数的概率
    • read 碱基的独热编码 (one-hot encoding) (4维)
    • 比对上的参考基因组碱基的独热编码 (5维，包含’N’)（可选特征）

2.2 深度学习模型推理 (Deep Learning Model Inference)

DeepMod2 提供了两种深度学习模型用于甲基化预测：BiLSTM 和 Transformer。

BiLSTM 模型包含两个隐藏单元数为128的 BiLSTM 层。BiLSTM 的输出被展平后，送入一个包含两个全连接层的网络，最终通过 sigmoid 激活函数得到甲基化概率。与只使用中心时间步输出的 DeepMod 不同，DeepMod2 使用了所有时间步的输出，从而能更好地整合上下文信息。

Transformer 模型将 BiLSTM 层替换为4个包含8个自注意力 (self-attention) 头的 Transformer 编码器层。

为了在 CPU 上实现快速推理，作者对模型进行了剪枝 (pruning)，移除了部分权重最低的连接。这使得模型在 CPU 上的推理速度提升了约5倍，而对准确率的影响可以忽略不计。

2.3 逐位点预测 (Per-site Prediction)

DeepMod2 将所有覆盖同一个 CpG 位点的 reads 的预测结果进行合并，以估计该位点的甲基化水平。

模型会分别报告参考基因组正链和负链的甲基化情况。如果输入的是一个经过分型 (phased) 的 BAM 文件（包含 HP 标签），DeepMod2 还能额外提供每个单倍型的甲基化计数。

最终的预测结果会以文本文件形式输出，并且 BAM 文件中也会添加标准的甲基化标签（MM 和 ML），方便在 IGV 等基因组浏览器中进行可视化。

3. 实验分析

作者使用多个公开可用的和新生成的数据集，在不同场景下对 DeepMod2 的性能进行了评估。

3.1 逐读长基准评估 (Per-read Benchmark Evaluation)

作者首先在 Ashkenazim 三人组（HG002, HG003, HG004）和 NIH3T3 小鼠细胞系的 R9.4.1 和 R10.4.1 流通池数据上，对 DeepMod2 和其他先进的甲基化检测工具（如 Guppy, Rockfish, f5c）进行了逐 read 的性能评估。

3.2 逐位点基准评估 (Per-site Benchmark Evaluation)

作者接着在全基因组尺度以及七个特定的基因组区域（CpG 岛、shores, shelves 等）中，对各个工具的逐位点预测性能进行了评估。

在 R9.4.1 数据集上，DeepMod2 BiLSTM 和 Rockfish 的性能始终最优，其全基因组 F1-score 均在 99.85-99.92% 的极窄范围内。在 R10.4.1 数据集上，DeepMod2 BiLSTM 和 Guppy 的性能相当，F1-score 在 99.94-99.97% 之间。

一个有趣的现象是，尽管某些工具（如 Nanopolish/f5c）的逐 read 性能较差，但它们的逐位点 F1-score 与其他工具的差距在0.5%以内，这表明通过聚合多个 reads 的信息可以有效弥补单个 read 预测的不足。

作者计算了各个工具预测的甲基化频率与短读长测序“金标准”之间的皮尔逊相关系数 (Pearson’s correlation coefficient)。结果显示，DeepMod2 BiLSTM 在 Ashkenazim 三人组上的相关性最高（95-97.35%）。热图显示，R10 流通池的数据与金标准的相关性更好，尤其是在中等甲基化水平的位点上。

3.3 简化代表性甲基化测序 (RRMS) 的评估

作者对 HG002 基因组进行了 RRMS 测序，并使用 DeepMod2 进行了分析。

RRMS 是一种可行的、经济高效的、用于对复杂基因组区域进行大规模甲基化谱分析的策略。

3.4 单倍型特异性甲基化检测

作者展示了 DeepMod2 在两个已知的印记控制区域（ICRs）中的单倍型特异性甲基化检测结果。DeepMod2 能够准确地识别出两个单倍型上相反的甲基化模式，这与已知的父源/母源印记模式完全一致。

作者进一步将这个流程应用于在 HG002, HG003, HG004 三个基因组中系统性地检测 ICRs。在21个已知的 ICRs 中，该流程成功地在三个基因组中分别找出了19, 20, 19个。将 DeepMod2 与 NanoCaller（一个 SNV 检测和分型工具）相结合，可以准确地识别印记区域。

3.5 运行时与准确性评估

作者对 DeepMod2 和 Dorado 两个框架的端到端运行时间，以及不同 basecaller 模型对准确率的影响进行了评估。

综上所述，DeepMod2 是一个功能全面、性能卓越且高效的开源工具，为基于纳米孔测序的甲基化研究提供了强大的支持。