论文解读：《基于注意力的多标签神经网络用于12种广泛存在的RNA修饰的综合预测和解释》

文章链接：https://www.nature.com/articles/s41467-021-24313-3
DOI：https://doi.org/10.1038/s41467-021-24313-3
期刊：Nature Communications（一区）
影响因子：14.919
发表时间：2021 年 6 月 29 日
数据：http://www.xjtlu.edu.cn/biologicalsciences/multirm
代码：https://github.com/Tsedao/MultiRM
服务器：http://www.xjtlu.edu.cn/biologicalsciences/multirm
补充资料：https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-021-24313-3/MediaObjects/41467_2021_24313_MOESM1_ESM.pdf

1.文章概述

最近的研究表明，通过转录后RNA修饰的表达转录调控对所有类型的RNA都是至关重要的。精确地确定RNA修饰位点对于理解RNA的功能和调控机制至关重要。作者提出一种方法：MultiRM，可以综合预测和解释转录后RNA修饰。该方法建立在一个基于注意力机制的多标签预测深度学习框架上，MultiRM不仅可以同时预测十二个广泛发生的RNA转录组修饰（m6A，m1A，m5C，m5U，m6Am，m7G，ψ，I，Am，Cm，Gm，Um）的位置，而且还可以返回对正面预测贡献最大的关键序列内容。作者提出的模型从相关序列上下文的角度考虑，揭示了不同类型RNA修饰之间的联系。作者还提供了检测多个RNA修饰的解决方案，使得对这些RNA修饰的综合分析成为可能，并且获得了对基于序列RNA修饰的机制更好理解。

2.背景

转录后修饰 RNA修饰增加了RNA分子的结构和功能多样性，调节RNA生命1的所有阶段。因此，准确地确定RNA修饰位点对于理解各种RNA的功能和调控机制至关重要。但是，现有的大多数研究只集中在单一的RNA修饰类型，未能通过综合预测模型支持多个同时RNA修饰预测。因此，研究不同修饰之间的相互作用是有限的。iRNA toolkit是最早也是最全面的RNA修饰预测方法，被广泛用作衡量不同RNA修饰预测方法性能的金标准。然而，iRNA toolkit是以多个独立研究的形式呈现的，每个研究都针对一个单一的修饰。iMRM web服务器是通过一个同时支持5个RNA修饰预测的方法，但是，它基于5个独立的二进制预测器，来分别对应5个RNA修饰，没有考虑不同修饰之间潜在的相互作用。鉴于不同RNA修饰固有的生化和生物物理特性，为一种修饰建立的预测框架通常可以方便地迁移到另一种修饰的预测上。因此，同时测试多个RNA修饰的计算框架是有益和有效的。
大多数现有的研究工作依赖于单一来源的有限数据（单一数据库或单一实验产生的数据集），未能充分利用现有的表观转录组信息。此外，以前已经报道过不同的表观转录组技术之间存在很大差异。因此，利用多种正交技术产生数据，尽可能减少这些数据集可用时潜在的技术偏差是至关重要的。
RNA修饰领域的大部分工作，都集中在预测精度上，但未能对其预测结果提供一个清晰和直观的解释。尽管一些现有的方法解释了他们训练有素的预测模型，但没有一个工作提供了对预测的过程有所解释。随着在可解释的RNA、DNA模型方面的进展，使得提取低水平的CNN内核并将它们可视化为位置权重矩阵成为可能。但是这些模式只能提供模糊的见解，特别是对于多层DNN，不能提供在核苷酸水平上的解释。然而，这对于正面RNA修饰预测的关键序列内容仍然具有重大意义，这将有助于理解个别RNA修饰位点的序列依赖性形成机制。
一些RNA修饰类型（如m6am）的预测框架还没有开发出来，尽管它通过表达谱分析技术检测已经开发出来，而且分析数据是公开的。由于这些原因，所以可以利用深度学习技术，开发一个统一的预测框架，通过整合多种技术生成的数据集，支持预测多个RNA修饰。
作者提出了MultiRM，一个以注意力为基础的多标签神经网络方法，可以综合预测和解释RNA修饰从原始RNA序列(或相应的DNA序列)中。作者开发的模型支持十二种RNA修饰类型，包括m6A、m1A、m5C、m5U、m6Am、m7G、Ψ、I、Am、Cm、Gm和Um。以上这些是比较广泛发生的RNA修饰，可以用现有的基础分辨率技术，对整个转录组进行分析，这是RNA修饰的一个非常理想的特性，可以用来进行可靠的大规模预测。作者提出的方法采用多标签架构，既能兼顾不同修改的共用结构，又能充分利用其独特的特性。针对多标签学习中训练数据不平衡的问题，利用在线硬实例挖掘（online hard examples mining，OHEM）和不确定性加权进行了改进。一些被广泛采用的最先进的机器学习算法：XGBoost和CatBoost也被实现作为基准测试。重要的是，作者使用积分梯度（integrated gradient，IG）和注意力加权来深入了解训练过的整体模型，并解释每个个体的预测。最后，开发了一个网络服务器，可以免费访问。

3.数据

作者从15种不同的基础分辨率技术中获得了20个表观转录组图谱，用于12种不同类型的RNA修饰（m6A、m1A、m5C、m5U、m6Am、m7G、Ψ、I、Am、Cm、Gm和Um），如下表所示。所使用的数据涵盖了所有广泛存在的RNA修饰，这些修饰可以在基本分辨率下，对整个转录组进行分析。为预测因子构建最可靠的阴性对照数据（非修饰核苷酸）。从含有阳性位点的同一转录本未经修饰的碱基中随机选择阴性位点。

肌苷（I）位点从RADAR数据库中收集，Am、Cm、Gm和Um修饰位点从Nm-seq中报道的RMBase中收集。当存在代表修饰的基序时，即m6A的DRACH基序和m6Am的BCA基序，该基序用于进一步限制相应修饰的正数据和负数据。对于m6A，由于需要检查使用各种技术的现有大量基础分辨率研究的可靠性，因此之前确定的总共87616个m6A位点仅通过将其从负m6A位点中排除而间接使用。
最终，共收集了30多万个位点。然后将每个类别中的RNA序列（对应于一种RNA修饰）分为三组，即训练集、验证集和测试集。在这里，不同类别（修饰类型）的训练集是不平衡的，即不同RNA修饰的位点数量不同，而验证集和测试集的样本大小分别为150和50。一般来说，超参数是基于验证集进行优化的，而报告的最终预测性能是在测试集上实现的。传统上，K-fold用于缓解许多机器学习问题中的过度拟合，尤其是对于那些训练数据较小的人。这是因为使用K-fold验证模型可以更好地估计模型的结果将如何推广到一个独立的数据集，尤其是在有限的数据集中，其中一个小测试无法反映数据的整个分布。然而，在作者的所搜集的数据中，总共有约30万个训练数据。因此，5%的测试集和验证集已经可以给出一个很好的估计。

4.方法

4.1 嵌入

（1）One-hot
A：（1，0，0，0）
C：（0，1，0，0）
G：（0，0，1，0）
U：（0，0，0，1）
N：（0，0，0，0）
最后每个长度为L的序列，被编码为L*4的形式。
（2）隐马尔科夫模型
将每个RNA序列x^m转换成n维向量，b：∑^L→R^d*L和g：R^d*L→Rⁿ，f(x^m)=g(b(x^m))和b(x^m)=[μ₁,μ₂,…,μ_L]。每个测试点都会对x^m中不同位置的潜在长程交互进行总结，g 则会对整个数据集的交互信息进行聚合和无固定维嵌入。在我们的工作中，我们在循环模块和多标签之前添加了 hmm 层模块，形成一个端到端的解决方案，从修饰位点的训练到预测。
（3）Word2vec
作者训练自己的RNA嵌入，把每个RNA序列作为一个句子，用k-mer作为句子中的单词。在Gene2vec中已经证明，3-mer在m6A位点具有最好的预测性能。因此，作者将输入数据改成3-mer形式。更具体地说，一个3碱基的滑动窗口在1001碱基的样本序列上移动，可以创建999个单词序列。每个单词对应于一个索引，该索引来自所有可能的3-mer组合（在作者的训练数据中有104种不同的组合）。然后，通过一个长度为五的相邻词窗口来学习相邻词之间的内在联系，生成一个300维的特征向量。最后，每个嵌入的RNA序列被转换成一个999乘以300矩阵。

4.2 模型搭建

首先根据样本的有效数量为每个任务使用一个恒定的权重。随后，作者设计了多标记模型，在训练过程中通过学习每个任务的权重来自我调节，只对损失较大（OHEM）的样本进行反向传播，以便在整个训练过程中共同对任务和样本进行优先排序。作者还测试了焦点损失的性能，它基本上降低了共同二元交叉熵损失分配给高分类样本的损失。

4.3 统计显著性。

预测概率的统计显著性是用p-value的上限来评估的，表示在同一核苷酸的所有出现中观察到的概率的极端程度。它是根据假定的RNA修饰位点的相对排序计算出来的，也就是说，如果只有1% 的核苷酸报告的概率大于某个特定位点，那么该位点的p-value的上界是0.01。这就是预测的截止点。但是，必须注意的是，即使采用相同的p-value截止值，假阳性预测的比例，在更丰富的RNA修饰和不丰富的RNA修饰之间，仍然存在本质上的差异。

4.4 解释

除了精确预测RNA修饰，作者使用注意力权重和积分梯度（IG）来直观地解释模型是如何做出决策的。具体来说，集中在所提出的模型最有价值的，同时作出不同的预测，并获得了具体核苷酸的贡献值，同时通过注意力的权重和树枝作出积极的预测。
Bahdanau注意力最初是作为处理序列到序列模型的长输入序列的解决方案而引入的。作者通过该方法将输入RNA序列映射到12个上下文向量，然后将其移植到本文所提出的方法中。由于它可以访问整个输入RNA序列，并且能够从序列中挑选出特定的元素来产生输出，因此这种机制使得模型可以根据需要，集中注意每个预测任务输入RNA序列的相关核苷酸。因此，通过可视化注意力的权重，表示每个预测任务中每个输入RNA序列核苷酸的权重，作者可以识别模型中输入序列中最关键的部分，同时做出不同的预测。
通过计算输出神经元相对于其输入的梯度，基于梯度的归因方法可以反映输入特征通过网络对特定输出的贡献程度。作者使用了一种叫做IG的归因方法。在这里，感兴趣的目标神经元是每个修饰的分类层。当输入沿从基线或参考到输入的线性路径变化时，IG计算输出神经元的平均梯度。它测量每个输入对修饰预测的贡献，并为输入序列中的重要核苷酸分配更高的分数。基于每个输入核苷酸位置的贡献分数，作者将归因图可视化为序列标志，其中高度代表该位置在预测中的重要性。正向核苷酸的大小代表了预测RNA修饰出现的一个重要水平。
对特定RNA修饰的每个输入序列的属性图进行可视化，不仅可以给出正面预测时的重要位置，而且可以揭示其相应修饰的潜在目标基序。为了计算每个RNA修饰的共识基序贡献最大，根据以前的研究，作者积累了对应于所有真正样本的每个位置的归因值，预测分数在前10% 。然后，对于每个样本，通过在所需长度的滑动窗口中取得最高的平均分，去掉它的邻居，然后在下一个主题中重复，来搜索整个属性图中的最高k主题。在多重序列比对之后，UMAP被用来嵌入最高级别的图案，DBSCAN被用来聚集这些嵌入的图案。最后通过计算它们的PWM将这些模块聚合在一起，并使用序列标志将它们可视化。

4.5 基线性能。

由于对于许多生物信息学问题，基于树的分类算法通常具有最好的现成准确性，在这项工作中，作者将所提出模型与两个梯度提升决策树（XGBoost和CatBoost）进行比较。XGBoost已经被广泛用于生物信息学预测。这里使用CatBoost是因为它具有处理范畴变量的内置技术，从而避免了在这个问题中将核苷酸转换为one-hot的额外步骤。梯度-升压决策树作为多级分类器用于所有13类，包括12个广泛发生的RNA修饰和不发生修饰核苷酸类。为了得到最优结果，作者搜索每种方法的超参数，并在补充表2中提供了最优选择。

4.6 DNN模型的注意力

作者分别对以注意力为基础的DNN模型进行了101碱基，51碱基，或21碱基 窗口的短RNA序列训练。对于word2vec，作者预先训练RNA序列的3-mer达到1001碱基，然后从中提取相应长度的短RNA子序列。然后将它们设计成端到端的方式，用one-hot编码输入序列。在训练过程中，使用了128的小批量作为输入，并训练了超过100次。另外，在训练中使用了Adam优化和128的小批量尺寸。同时，对合适的模型进行了学习率衰减，并在五个连续时间段内泛化损失增加时引入提前停止，以防止对训练数据的过度拟合。最后，使用验证集搜索单个模型的最佳超参数，并使用测试集从各种具有最佳性能的模型中选择最佳模型。

5.结果

5.1 参数优化

作者以21碱基，51碱基，和101碱基的RNA序列作为输入评估了多标记模型，如表1所示，51碱基序列的输入在所有修饰中获得了最佳的平均性能，并且在测试的十二个RNA修饰中，这种设置也返回了最佳性能。

如补充表1所示，输入序列的51碱基对于XGBoost方法也是最佳的。

然后，针对训练数据不平衡的问题，在优化的多标记模型上实现了OHEM、不确定性加权（UW）和焦点损失，并用51碱基的输入序列测试了它们的性能。然后根据优化的模型得到改进的，如表2所示。OHEM和UW都有利于改进。所以作者将两者进行了结合。

在优化设置(51碱基输入，UW、OHEM)的基础上，比较了新开发的方法MultiRM与基线方法和其他嵌入技术。每个模型的优化超参数可以在补充表2中找到。

如表3所示，新提出的方法MultiRM获得了最好的平均性能和中位性能，并且在十二个RNA修饰方面获得了最好的平均性能。采用的XGBoost算法在4个改进方案上取得了最佳性能，在平均性能和中位性能略有不足。

根据各自的ROC曲线选择最大G-Mean值的最佳阈值。表4计算并提供了相应的性能评估指标，包括每种修饰的SN、SP、ACC、MCC。



多元回归法的精度和回收率曲线和 ROC曲线曲线见补充文件图1和图2。关于多标签分类方案下的MultiRM性能指标，请参阅补充表3。

5.2 解释性

注意力范围内的序列内容对 rna 修饰有较大的影响，这些区域内的突变更有可能导致RNA修饰位点的增加或减少，如补充文件图3所示。

为了进一步量化MultiRM和DREME、STREME获得的序列之间的相似性，应用序列比较工具TOMTOM生成p-value。如图2所示，充分的p-value表示在一定程度上具有一致性。

如图3所示，RNA修饰显示出彼此之间具有强烈而显著的正相关性，包括那些来自不同核苷酸的修饰。

为了进一步验证上述发现，作者计算了两个任意RNA修饰之间的成对距离，并与随机值进行了比较。虽然不能完全排除实验偏差的可能性，例如，多聚体选择，但观察到在这个分析中考虑的大多数RNA修饰中有很强的聚集效应。很明显，两个任意RNA修饰之间的距离可能比随机的更近，如补充图4所示。

5.3 Web

6.讨论

在补充图5中提供了关于不平衡样本大小的建议的MultiRM的性能指标，这反映了目前对修饰在现实世界中的分布的知识，但是，必须注意的是，收集到的RNA修饰位点的数量受到生物技术的检测灵敏度和现有实验数据的强烈影响，而不是它们的真实数量。当这些RNA修饰的整体丰度更容易获得时，将需要更可靠的错误发现率控制。

7.想法

1. 单修饰不同物种间的多标签预测是否有意义？
2. CatBoost可以应用一下；
3. 正负样本不均衡时，可以试一下OHEM方法对样本进行平衡；
4. 此文为RNA修饰，是否可以做一下DNA和RNA的之间同类型碱基的修饰，例如DNA中的6mA和RNA中的m6A，同样都是腺嘌呤发生修饰。