Deep6mA：探索不同物种DNA N6甲基化腺嘌呤位点相似性模式的深度学习框架.
0. TL; DR
本文开发了一个名为 Deep6mA 的深度学习框架，用于识别 DNA 6mA 位点，且无需任何关于 6mA 的先验知识或手动设计的序列特征。
使用水稻数据训练的模型在预测其他三个物种（拟南芥、草莓、月季）的 6mA 位点时，预测准确率均超过90%。作者发现 6mA 倾向于出现在 GAGG 基序中，这意味着 6mA 位点附近的序列可能是保守的。研究还发现 6mA 在启动子的 TATA 盒中富集，这可能是其调控下游基因表达的主要机制。
1. 背景介绍
DNA 甲基化修饰，如 N4-甲基胞嘧啶 (4mC)、N6-甲基腺嘌呤 (6mA) 和 5-甲基胞嘧啶 (5mC)，在不改变序列的情况下，对基因表达的表观遗传调控起着重要作用。其中，6mA，即腺嘌呤第6位氮原子的甲基化，近年来被发现是真核生物 DNA 中一种重要的表观遗传修饰。它在 DNA 修复、DNA 复制、基因转录调控和基因表达调控等过程中扮演着重要角色。
尽管 6mA 位点在基因组中并非均匀分布，并可能受环境因素影响，但它在原核生物和真核生物中表现出相似的特征。随着高通量测序技术的发展，多种实验技术（如 MeDIP-seq, HPLC-MS/MS, SMRT）已被用于研究 6mA 的分布及其潜在功能。例如，研究发现水稻基因组中约0.2%的腺嘌呤被 6mA 甲基化，并且富含 GAGG 的序列是其最显著的富集区域。
为了辅助实验研究，一系列计算工具被开发出来用于预测 6mA 位点。早期的模型，如 i6mA-Pred（基于支持向量机）和 MM-6mAPred（基于马尔可夫模型），在利用核苷酸化学性质和邻近核苷酸转移概率方面进行了探索。后续的模型，如 iDNA6mA-rice（基于随机森林）和 SNNRice6mA（基于卷积神经网络），则引入了更复杂的机器学习和深度学习方法。
然而，这些方法仍存在一些局限性。一个关键的难点在于如何有效地捕捉 DNA 序列中核苷酸之间的长程依赖结构。CNN 虽然在图像识别等领域表现出色，但在学习长距离依赖信息方面能力有限。循环神经网络，特别是长短期记忆网络 (LSTM)，则更擅长处理序列数据中的时序信息。
因此，作者提出了一个名为 Deep6mA 的新型深度学习框架，它通过整合 CNN 和双向LSTM (BLSTM)，旨在更有效地捕捉 DNA 序列的特征和依赖关系，从而实现对 6mA 位点的高精度预测。
2. Deep6mA方法
Deep6mA 的核心是一个结合了卷积神经网络 (CNN) 和双向长短期记忆网络 (BLSTM) 的深度学习框架，旨在充分捕捉 DNA 序列中的高级特征和长程依赖关系。
2.1 基准数据集
作者使用了两个主要的水稻 6mA 数据集进行模型的训练和评估：
- 6mA-rice-Lv 数据集：包含154,000个阳性样本（含 6mA 位点）和154,000个阴性样本。每个样本都是一个长度为41 nt的序列，6mA 位点位于中心。
- 6mA-rice-Chen 数据集：一个较小的数据集，包含880个阳性样本和880个阴性样本。
此外，为了验证模型的跨物种泛化能力，作者还收集了拟南芥 (Arabidopsis thaliana)、草莓 (Fragaria vesca) 和月季 (Rosa chinensis) 的 6mA 数据集。
2.2 序列表示
作者没有使用手动设计的 DNA 序列特征，而是采用了简单而有效的独热编码 (one-hot encoding) 方法将序列转换为编码张量。A, C, G, T 分别被编码为 $(1,0,0,0), (0,1,0,0), (0,0,1,0), (0,0,0,1)$。代表不确定核苷酸的 N 被编码为 $(0,0,0,0)$。
通过这种方式，一个长度为41 nt的输入 DNA 序列被表示为一个 $4 \times 41$ 的编码矩阵。这个矩阵可以被看作是一个图像，从而启发了作者使用深度学习框架进行处理。
2.3 Deep6mA 模型架构
Deep6mA 的模型架构如图所示，它由五个 CNN 层、一个 BLSTM 层和一个全连接层组成。
CNN 层被用来提取 DNA 序列中的高级特征。第一层卷积层扮演着 6mA 位点周围基序检测器的角色，而后续的卷积层则能捕捉到更高层次的抽象特征。
模型包含五个卷积层，每个层都配置了256个大小为10的滤波器。使用 ReLU 作为激活函数。每个卷积层之后都跟一个最大池化层，用于优化特征的冗余并防止过拟合。
在 CNN 提取了高级特征之后，这些特征被送入一个 BLSTM 层。BLSTM 能够从前向和后向两个方向上学习序列中的上下文依赖关系，这对于捕捉 DNA 序列中核苷酸之间的长程依赖至关重要。BLSTM 层的隐藏单元大小为32，激活函数为 tanh。
BLSTM 层的输出被传递到一个包含32个隐藏单元的全连接层，用于整合所有学习到的特征。最后，一个 sigmoid 激活函数被用来结合 FC 层的输出，做出最终的二元分类决策（即该位点是否为 6mA 位点）。
2.4 模型训练
使用二元交叉熵 (binary cross-entropy) 作为损失函数，并加入了一个 L2 正则化项来避免过拟合。
\[L(w) = -\sum_{i=1}^{N} [y_i \log(y'_i) + (1-y_i)\log(1-y'_i)] + \alpha \|w\|^2\]其中，$y_i$ 是真实标签，$y’_i$ 是模型的预测值。
使用 Adam 优化器进行训练。在每个卷积层后都应用了批归一化和 dropout（比率为0.5）。学习率设为0.01，批大小为256。训练过程中使用了early stopping 策略，如果验证集上的性能连续5个周期没有提升，则停止训练。
3. 实验分析
作者通过一系列详尽的实验，评估了 Deep6mA 的性能，并与现有的先进方法进行了比较。
3.1 CNN 与 CNN+LSTM 的性能比较
作者首先比较了单独使用 CNN 和结合 CNN 与 LSTM 的性能。在不同的 CNN 层数和卷积核大小设置下，CNN+LSTM 框架的性能始终优于单独的 CNN。这表明，LSTM 能够有效捕捉序列的依赖结构，从而提升了模型的预测能力。
3.2 Deep6mA 模型参数选择
作者通过五折交叉验证，在30种不同的参数设置下对 CNN+LSTM 框架进行了评估，以选择最优的模型参数。最终，一个包含5个卷积层、256个卷积核、卷积核大小为10、LSTM 隐藏单元为32的配置被选为最终的 Deep6mA 模型。
3.3 6mA 位点的定位特征
作者分析了水稻12条染色体上相邻 6mA 位点之间的距离分布。结果发现，不同染色体上的分布模式相似，并且相邻 6mA 位点之间的平均距离大于64 nt。这表明，6mA 位点很少像 5mC 那样在连续区域内聚集形成差异甲基化区域 (DMR)。
然而，当作者深入观察那些确实聚集在一起的 6mA 位点时（即在30 nt的序列中有超过5个 6mA 位点），他们发现这些 6mA 位点几乎都位于启动子的 TATA 盒区域。TATA 盒是基因转录起始的关键功能元件。6mA 在此处的富集可能直接影响下游基因的表达，这可能是 6mA 甲基化修饰的一个重要调控功能。
3.4 与其他先进方法的比较
在大规模的 6mA-rice-Lv 数据集上，作者将 Deep6mA 与 MM-6mAPred 和 SNNRice6mA 进行了比较。
Deep6mA 在所有指标（SN, SP, ACC, MCC, AUC）上均表现最佳。其 AUC 达到了0.98，显著高于其他两种方法。作者认为，Deep6mA 的优越性能主要归功于 BLSTM 能够学习到核苷酸之间的长程依赖结构。
在较小的 6mA-rice-Chen 数据集上，Deep6mA 的性能同样优于其他两种方法。2
3.5 在其他植物物种上的验证
基于后续的基序分析发现 6mA 在不同物种间具有保守性，作者进一步验证了使用水稻数据训练的 Deep6mA 模型在预测其他植物物种 6mA 位点上的能力。
将在水稻数据上训练好的 Deep6mA 模型直接应用于拟南芥、草莓和月季的数据集进行预测。结果显示，Deep6mA 在这三个物种上都取得了非常高的预测准确率，并且其性能全面优于 SNNRice6mA 和 MM-6mAPred。
这有力地证明了 6mA 位点序列模式在不同物一间是保守的，并且 Deep6mA 学习到的模式具有很强的泛化能力。
3.6 跨物种的基序比较
为了探究 6mA 序列的保守性，作者使用 MEME 算法分析了四个物种中 6mA 位点周围的序列基序。结果清晰地显示，GAGG 是这四个物种中最显著的共同基序。这与之前的研究发现一致，并表明 DNA 6mA 甲基化在不同物种中很可能频繁地发生在 GAGG 基序上。
作者还通过 TOMTOM 工具，将 Deep6mA 模型第一层卷积核学习到的基序与已知的基序进行比对。结果发现，模型学习到的最显著的基序中也包含了 AGG，这再次证实了模型确实捕捉到了在不同物种间保守的序列模式。
