CNN6mA: 基于位置特异性CNN与跨交互网络的可解释性6mA位点预测神经网络模型.

0. TL; DR

本文开发了一个基于 CNN 的 6mA 位点预测器CNN6mA。CNN6mA 在每个位置上应用具有不同窗口大小的位置特异性一维卷积层 (position-specific 1-D convolutional layer)提取与特定位置相关的、多层次的特征。交叉互动网络 (cross-interactive network)探索了查询序列内所有核苷酸模式之间的关系，通过将特征在两个不同维度上进行广播和求和，来捕捉特征之间的交互信息。

在多个物种的基准测试中，CNN6mA 的性能优于现有的先进模型。更重要的是，模型能够生成一个贡献度得分向量 (contribution score vector)，它明确地指出了序列中哪些位置对预测结果贡献最大，从而为模型的预测机制提供了清晰的解释。

1. 背景介绍

N6-甲基腺嘌呤 (6mA) 是一种在不改变 DNA 序列的情况下调控遗传活动的表观遗传修饰。它在 DNA 复制、修复、沉默、转录和基因表达等多种表观遗传过程中扮演着至关重要的角色，并与癌症等疾病的发生发展密切相关。

尽管已存在多种实验方法（如单分子实时 (SMRT) 测序）用于在全基因组范围内以单碱基分辨率识别 6mA 位点，但这些方法通常耗时、费力，且有时测序质量和信噪比较差。

为了补充实验方法并加速对 6mA 调控机制的理解，基于机器学习 (machine learning) 和深度学习 (deep learning) 的计算机模拟 (in silico) 方法应运而生。

然而，尽管这些方法取得了成功，但它们并未深入探讨其预测机制。例如，BERT6mA 忽略了每个核苷酸的位置信息，而这在 DNA 修饰预测中可能非常重要。

为此，作者开发了一个名为 CNN6mA 的新型 CNN 预测器。该模型通过两种新颖的架构——位置特异性一维卷积层和交叉互动网络——来分别捕捉与位置相关的特征和特征间的相互关系，旨在实现高性能的同时，提供对预测机制的清晰解释。

2. CNN6mA 方法

2.1 数据集

作者使用了从 Lv 等人的 iDNA-MS 服务器获取的、覆盖11个不同物种（如拟南芥、线虫、果蝇、人类等）的 6mA 和非 6mA 序列数据集。所有序列长度均为41 bp。

2.2 深度学习模型 (CNN6mA)

作者提出的 CNN6mA 模型是一个基于 CNN 的深度学习模型，其整体架构如图所示，主要由四个子网络组成。

2.2.1 嵌入层 (Embedding Layer)

输入一个长度为 $N$ 的 DNA 序列。将 A, T, G, C 四种核苷酸分别转换为索引 1, 2, 3, 4。使用 nn.Embedding (PyTorch) 将每个索引表示为一个 $M$ 维的特征向量（$M$ 设为64）。输出一个 $N \times M$ 的特征矩阵 $X$。

2.2.2 位置特异性一维卷积层 (Position-specific 1-D convolutional layer)

与在所有位置共享相同滤波器的传统 1-D CNN 不同，该层为每个序列位置 $i$ 都创建了一组独立的滤波器。

在每个位置 $i$ 创建 $F$ 个窗口大小为 $W$ 的滤波器 $f_i^{(W)}(k,l,j)$。将这些位置特异性的滤波器应用于输入特征矩阵 $X$，生成输出矩阵 $C$。

\[c_i^{(W)}(k) = \sum_{l=1}^{W} \sum_{j=1}^{M} f_i^{(W)}(k,l,j) x(i-W/2-1+l, j)\]

为了捕捉不同层次的特征模式，作者使用了多种不同窗口大小（$W \in {3, 5, 7}$）的滤波器，并将它们的输出在每个位置进行拼接，得到最终的特征矩阵 $B$。

\[B[i,k] = \text{concat}(d_i^{(W_1)}(k), d_i^{(W_2)}(k), \dots, d_i^{(W_S)}(k))\]

2.2.3 交叉互动网络 (Cross-interactive network)

交叉互动网络旨在探索查询序列内所有核苷酸模式之间的关系。

将上一层输出的特征矩阵 $B$ (形状为 $N \times F_{cat}$) 在两个不同的新维度上进行广播和扩展，得到两个形状分别为 $N \times 1 \times N \times F_{cat}$ 和 $1 \times N \times N \times F_{cat}$ 的张量。将这两个扩展后的张量相加，得到一个考虑了所有成对核苷酸模式之间关系的特征张量 $H$。对 $H$ 进行最大池化 (max-pooling)，得到最终的特征矩阵 $Q$。

2.2.4 输出层

将特征矩阵 $Q$ 送入一个包含三个全连接层的输出网络。最终的输出是一个贡献度得分向量 (contribution score vector)，其每个元素的值都在0到1之间，代表了序列中对应位置的核苷酸对最终预测的贡献度。通过对贡献度得分向量进行全局最大池化 (global max-pooling)，得到最终的预测得分。如果得分大于0.5，则判定为阳性样本。

3. 实验分析

3.1 CNN6mA 的优化

作者首先探究了位置特异性一维卷积层中多尺度滤波器的最佳组合。作者测试了窗口大小为{3}, {3,5}, {3,5,7}, {3,5,7,9}的四种组合。结果显示，使用窗口大小为 {3, 5, 7} 的三个滤波器组合时，模型在 SN, ACC, MCC 和 AUC 等指标上均取得了最佳性能。多尺度滤波器的组合对于提升预测性能是有效的。窗口大小为9可能过长，无法有效捕捉核苷酸模式。

3.2 CNN6mA 与现有先进方法的比较

作者将 CNN6mA 与六种现有的先进模型（iDNA-MS, iDNA-ABT, SNNRice6mA, DeepTorrent, Deep6mA, BERT6mA）在11个物种的数据集上进行了性能比较。

在11个物种中，CNN6mA 在6个物种上取得了最高的 ACC，在8个物种上取得了最高的 AUC。与基于 BERT 的 BERT6mA 和 iDNA-ABT 相比，CNN6mA 在许多物种上表现相当或更优。

CNN6mA 具有很高的学习能力和泛化能力，其性能在多个物种上均优于或可比于现有的最先进模型。这主要归功于其能够捕捉位置相关信息和局部核苷酸模式之间关系的创新架构。

3.3 位置特异性一维卷积层的有效性

为了验证位置特异性卷积层的有效性，作者将其与使用共享滤波器的普通卷积层进行了性能比较。

结果显示，在几乎所有物种和两种不同的编码方法（索引嵌入和独热编码 (one-hot encoding)）下，使用位置特异性卷积层的模型的 ACC 和 AUC 都显著高于使用普通卷积层的模型。这明确地表明，对于 6mA 的预测，位置特异性的信息是至关重要的，而位置特异性卷积层能够有效地捕捉这些信息。

3.4 交叉互动网络的有效性

作者比较了包含和不包含交叉互动网络的 CNN6mA 模型的性能。结果显示，在保持模型参数数量大致相同的情况下，加入交叉互动网络显著提升了模型的预测性能。交叉互动网络通过探索序列内不同核苷酸模式之间的关系，为模型提供了额外的重要信息，对于 6mA 位点的预测是有用的。

3.5 预测机制的可解释性

CNN6mA 的一个重要特性是它能够生成一个贡献度得分向量，用于解释其预测机制。

作者使用 UMAP 和 t-SNE 对阳性和阴性样本的平均贡献度得分向量进行了可视化。结果显示，两类样本的向量被清晰地分开了，表明这些向量是分类的关键。

热图显示，阳性样本的贡献度得分向量普遍高于阴性样本。更重要的是，在阳性样本中，得分较高的位置与使用 sequence logo 工具识别出的保守核苷酸的位置高度一致，都集中在序列的中部（18-28位）。

作者进一步分析了在特定位置得分较高的样本。结果显示，这些样本在该位置周围确实富集了特定的核苷酸模式。例如，在第20位得分高的序列富集了 GAGG 基序，这与之前的研究报道相符。

贡献度得分向量能够准确地反映出对 6mA 预测有重要贡献的核苷酸位置和模式。与 BERT6mA 的注意力权重相比，贡献度得分向量是基于每个核苷酸计算的，因此能够更直观、更明确地揭示 6mA 相关的模式，为发现新的生物学基序提供了有力工具。