一种在水稻基因组中进行DNA N6腺嘌呤位点识别与基序提取的卷积计算方法.

0. TL; DR

i6mA-CNN 是基于卷积神经网络的工具，能够识别水稻基因组中的 6mA 位点。采用了一个包含四个并行分支的 CNN 架构。每个分支独立处理一种不同类型的特征编码，使得模型能够从多个角度学习和协调异构特征，避免了手动特征拼接和优化的局限性。

在基准数据集上，使用五折交叉验证，i6mA-CNN 取得了0.98的 auROC（ROC 曲线下面积）和93.97%的总体准确率，其性能优于所有现有的先进工具（如 SpineNet-6mA, iDNA6mA-Rice 等）。在另外三种植物基因组（草莓、月季、拟南芥）的独立测试集上，i6mA-CNN 同样表现出色，证明了其在植物基因组中识别 6mA 位点的泛化能力。

1. 背景介绍

DNA 表观遗传修饰，如 N6-甲基腺嘌呤 (6mA)，是一种在 DNA 复制后发生的修饰，它通过将甲基基团转移到腺嘌呤核苷酸的第六个位置来实现。这种修饰对于区分原始 DNA 链与新合成的 DNA 链、基因转录调控、复制和修复等多种生物学功能至关重要。

尽管存在多种实验方法（如 SMRT-seq）来识别全基因组的 6mA 位点，但这些方法通常耗时且成本高昂。幸运的是，计算生物学为大规模、高效地识别 6mA 位点提供了可能。

近年来，已有一系列针对水稻基因组的 6mA 位点预测工具被开发出来。

然而，尽管这些方法取得了良好的性能，但它们大多依赖于单一类型的特征表示（单一特征向量或单一特征矩阵），缺乏学习和协调多种异构特征之间关系的能力。这限制了模型从数据中提取更全面信息的能力。

为此，作者提出了 i6mA-CNN，一个基于多分支 CNN 的深度学习模型。该模型能够并行处理多种不同类型的特征表示，并自动学习它们之间的协调关系，从而实现更准确的 6mA 位点预测。

2. i6mA-CNN方法

作者提出的 i6mA-CNN 是一个基于 1D CNN 的多分支分类器，能够识别植物基因组中的 6mA 位点。其核心在于并行处理四种异构的特征表示。

2.1 基准数据集构建

从 NCBI 的 GEO 数据库获取了265,290条包含 6mA 位点的序列。使用 CD-HIT (80%) 去除冗余序列，最终得到154,000条阳性样本。从 NCBI 中构建阴性样本，并确保其 GAGG 基序的比例与阳性样本相似，以避免模型学习到这种简单的偏置。最终，构建了一个包含154,000个阳性样本和154,000个阴性样本的平衡数据集。

2.2 序列表示 (特征编码)

作者采用了四种不同的方法来表示 DNA 序列，每种方法都旨在捕捉序列的不同方面信息。

  1. PseAAC 启发的特征表示：受伪氨基酸组成 (PseAAC) 的启发，作者为每个核苷酸定义了四个维度的特征：环结构 (ring structure)、官能团 (functional group)、氢键 (hydrogen bonding) 和滞留密度 (lingering density)。一个41 bp的序列最终被表示为一个 $41 \times 4$ 的特征矩阵。
  2. 单体独热表示 (Monomer one hot representation)：将 A, T, C, G 四种核苷酸分别编码为4维的独热 (one-hot) 向量。一个41 bp的序列被表示为一个 $41 \times 4$ 的矩阵。
  3. 二聚体独热表示 (Dimer one hot representation)：将所有16种可能的二核苷酸 (dimer)（如 AA, AT…）分别编码为16维的独热向量。一个41 bp的序列包含40个重叠的二核苷酸，因此被表示为一个 $40 \times 16$ 的矩阵。
  4. 二聚体理化性质 (DPP) 表示：作者从已有的90种二核苷酸理化性质 (dimer physicochemical properties, DPP) 中，通过一种基于相关性的特征选择方法，筛选出了9种最重要且冗余度最低的性质。因此，一个41 bp的序列被表示为一个 $40 \times 9$ 的矩阵。

2.3 模型架构

i6mA-CNN 的模型架构如图所示，它由四个并行的分支组成，每个分支对应一种上述的特征表示。

每个分支都是一个包含三个顺序连接的一维卷积层 (1D convolution layers) 的 CNN。每个分支独立地从其对应的输入特征矩阵中学习类别区分特征。

每个分支的最后，CNN 的输出被展平并通过一个带 Sigmoid 激活函数的全连接 (dense) 层，输出一个介于0到1之间的值，作为该分支的特征权重代表 (feature weight representative)。

将来自四个分支的四个特征权重代表拼接成一个长度为4的向量。这个4维向量随后被送入一个最终的全连接层（带 Sigmoid 激活函数），输出最终的分类概率。

2.4 潜在基序发现 (Potential Motif Finding)

为了使模型具有可解释性，作者还设计了一个基于注意力 (attention) 的扩展模型，用于发现潜在的基序。

该模型由两部分组成：i6mA-CNN 的特征学习器部分，和一个双向 LSTM (bidirectional LSTM) 子模型。LSTM 的每个时间步接收一个核苷酸的独热表示。其输出矩阵 $O$ (维度为 $41 \times 256$) 与 CNN 部分的输出向量 $C$ (维度为 256) 进行点积，得到一个长度为41的注意力分数向量。

这个向量经过 Softmax 层归一化后，得到一个概率向量 $S$。$S$ 中的第 $i$ 个值代表了输入序列中第 $i$ 个核苷酸对于 6mA 识别的“重要性”或“激活值”。

当一个阳性样本输入到训练好的注意力模型中时，如果其输出的概率向量 $S$ 中，某个位置的概率值大于平均值（$1/41$），则该位置的核苷酸被认为是“激活的”。连续的“激活”核苷酸构成一个“激活区域”。如果在大量的阳性样本中，某个特定的核苷酸子序列频繁地作为“激活区域”出现，那么这个子序列就被识别为一个潜在的基序。

3. 实验分析

作者对 i6mA-CNN 的性能进行了全面的评估，并与多个现有的先进工具进行了比较。

3.1 与现有先进工具在基准数据集上的比较

作者在两个基准数据集上，将 i6mA-CNN 与11个现有的先进工具进行了性能比较。

在 Dataset 1 上 (880样本)使用留一交叉验证 (Jackknife testing)，i6mA-CNN 的 MCC (0.85) 和 auROC (0.98) 显著优于 i6mA-Pred (0.67, 0.91)。使用10折交叉验证，其性能也全面优于或等于其他工具。

在 Dataset 2 上 (154,000样本)使用5折交叉验证，i6mA-CNN 的 MCC (0.88) 和 auROC (0.98) 达到了最高水平，优于 iDNA6mA-Rice (0.83, 0.96) 和 SNNRice6mA (0.84, 0.97)。

展示了 i6mA-CNN 在基准数据集上进行5折交叉验证的 ROC 曲线。五条曲线高度重叠，且都远高于随机预测线，表明模型性能优异且稳定。

3.2 在独立测试集上的性能

为了检验模型的泛化能力，作者在另外三种植物（草莓、拟南芥、月季）的独立测试集上评估了 i6mA-CNN。在草莓和拟南芥数据集上，模型的成功率（准确率）分别达到了97.78%和95.78%。尽管在样本量较小的月季数据集上性能有所下降，但仍然达到了85.79%。i6mA-CNN 不仅在水稻上表现出色，还对其他植物物种具有良好的泛化能力。

3.3 在不平衡训练数据集上的模型验证

为了模拟真实世界中数据不平衡的情况，作者通过对阳性样本进行欠采样 (undersampling)，构建了不同不平衡比例（1:1 到 1:20）的训练集，并评估了模型的性能。

结果显示，即使在正负样本比例为1:20的极端不平衡情况下，模型的 auROC 仍然能达到0.91，MCC 也能达到0.73。i6mA-CNN 模型对于数据不平衡问题具有相当的鲁棒性。

3.4 特征重要性分析

作者通过分析特征学习器部分的中间输出，探究了四种不同特征对分类的贡献。

作者还进行了消融实验，即只使用四种特征中的一种来训练模型。结果显示，任何单一特征都无法达到使用全部四种特征时的性能。不同的特征对于识别不同类别的样本具有不同的重要性。i6mA-CNN 的多分支架构能够有效地学习和协调这些异构特征，从而取得比任何单一特征或简单集成学习更好的性能。

3.5 识别出的潜在基序

作者使用其提出的基序发现算法，识别出了一系列潜在的基序。排名前五的基序包括 ATAT, AGGC, GAGG, AAAA, AAAT。值得注意的是，已被实验验证的最常见的 6mA 基序 GAGG，在作者的计算方法中被识别为第三重要的基序。这证明了作者提出的基于注意力的基序发现算法的有效性，并为未来的实验验证提供了有价值的线索。