CN115394376B - 一种预测环状rna ires的方法 - Google Patents

Abstract

本发明公开了一种基于DeepCPI计算框架将特征嵌入与深度学习方法相结合实现预测环状RNA IRES的方法。在预测过程中，该方法使用S‑LSTM学习全局句子级节点来表征整个RNA序列，从而为分类任务建模上下文信息。此外，该方法通过一个加权RNA图来模拟RNA的二级结构，并基于GCN提取所需结构特征。与此同时，该方法采用下采样和模型集成解决数据不平衡带来的问题。验证实验结果表明，与目前各种深度学习方法相比，本发明所述方法不仅在序列和结构中提取特征方面具有最佳的性能，而且能有效预测不同长度的人和病毒环状RNA IRES。由此可见，本发明所提供的方法可以应用于环状RNA编码能力研究以及随后的环状RNA设计和应用，如环状RNA疫苗的开发。

Description

一种预测环状RNA IRES的方法

技术领域

发明涉及计算机辅助生物分子设计领域，具体涉及一种基于DeepCIP计算框架将特征嵌入与深度学习方法相结合实现预测环状RNA IRES的方法。

背景技术

借助于深度测序和计算分析，研究人员发现环状RNA（环状RNAs）是一大类具有多种功能的RNA，不仅存在于动植物中，也存在于病毒中（Nature 2013， 495:333–338;Front. Plant Sci. 2019,10:379; Nat. Commun. 2019; 10:2300）。虽然环状RNAs通常被认为是一类非编码RNA，具有共价键的环状结构，没有5'帽和3'端PolyA尾，但越来越多的证据表明，环状RNAs具有蛋白质编码功能(JNCI J. Natl. Cancer Inst. 2018, 110:304–315; Oncogene 2018, 37:1805–1814; Mol. Cell 2017, 66:22-37.e9 )。环状RNA编码的蛋白质可能普遍存在于多种疾病中，特别是在人类癌症中发挥抗肿瘤或促进肿瘤的作用(Mol. Cancer 2021, 20:119; Mol. Cancer 2020, 19:71; Mol. Cancer 2021; 20:66)。这一发现表明了环状RNAs翻译的潜在发展和临床应用。与线性mRNA相比，环状RNA具有更好的热稳定性、更长的有效时间和更特异的组织表达，被认为是开发下一代mRNA药物的理想平台(Mol. Ther. - Nucleic Acids 2021,25:355–371; J. Controlled Release 2022,348:84–94)。近来，有研究团队开发了编码SARS-CoV-2及其变体的刺突蛋白受体结合域（RBD）的新冠环状RNA疫苗( Cell 2022, 185:1728-1744.e16)。此外，环状RNAs在预防或治疗传染病和自身免疫病方面也有广泛的应用。

环状RNAs的翻译需要内部核糖体进入位点（IRESs），IRESs是一种RNA顺式作用的调节元件，可以在没有5'帽的情况下将小核糖体亚单位招募到翻译起始位点(Nat.Commun. 2018, 9:4475; Mol. Cell 2017, 66:9-21.e7)。IRES首次在RNA病毒基因组中发现，例如脊髓灰质炎病毒（PV）和脑心肌炎病毒（EMCV）的基因组(Nature 1988, 334:320–325; J. Virol. 1988; 62:2636–2643)。借助双顺反子分析，IRES已在病毒和细胞mRNA中广泛发现(J. Biol. Chem. 2005, 280:23425–23428)。与cap依赖性翻译相比，IRES介导的翻译机制相对未知。然而，研究人员认为，根据结构RNA对核糖体的直接招募，或在典型起始因子以及额外的IRES反式作用因子（ITAFs）的协助下与核糖体部分间接的相互作用，一级序列和RNA结构对IRES活性具有重要的功能(Oncogene 2004, 23:3200–3207; Front.Oncol. 2015, 5)。此外，与线性RNA IRES相比，环状RNA IRES含有更高的GC含量、更低的最小自由能（MFE），并且通常结构更为合理(Mol. Cell 2021; 81:4300-4318.e13)。不同的特性导致线性RNA系统和循环RNA系统中IRES活性不同。

据估计，大约10%的mRNA可能通过IRES介导的cap独立机制启动翻译，一些研究报告称IRES在多种细胞过程和多种疾病中发挥重要作用(Oncogene 2004, 23:3200–3207;Front. Oncol. 2015, 5)。然而，识别IRES元素的传统实验方法费时费力，导致确定的IRES数量仍然非常有限，这阻碍了对IRES共同特征和功能的研究。数据驱动的计算方法，如机器学习，正越来越多地应用于生物数据。从生物信息学的角度来看，IRES的预测是困难的，因为所有IRES都没有普遍的保守特征，但仍有一些预测工具被开发出来。Kolekar等人(Sci.Rep. 2016, 6:27436) 开发了一种基于支持向量机（SVM）的预测器，称为IRESPred，使用35个特征进行病毒和细胞IRES预测。然而，IRESPred中使用的训练数据集的阴性样本不是经实验验证的IRES阴性序列。2016年，Weingarten Gabbay等人(Science 2016; 351:aad4939)设计了高通量双顺反子分析，并鉴定了数千个新的人类和病毒IRES序列。新的可用IRES序列数量的大幅增加使得机器学习方法能够更好地应用于IRES识别。Gritsenko等人(Comput. Biol. 2017; 13:e1005734)基于Weingarten Gabbay等人的高通量分析数据集开发了IRESpredictor，一种随机梯度增强随机森林回归模型，用于使用6120个全局和局部序列k-mer特征预测IRES活动。但大量特征可能会导致模型过拟合和训练时间慢的问题。随后，Zhao等人(J. Genet. Genomics 2018; 45:403–406)提出了IRESfinder模型，该模型具有精心选择的19个k-mer序列特征，仅使用Weingarten Gabbay等人数据集的人类IRES来预测细胞IRES。Wang等人(BMC Bioinformatics 2019; 20:409)开发了一个名为IRESpy的XGBoost模型，该模型是最新的IRES预测工具，具有更好的性能和更短的训练时间，基于Weingarten Gabbay等人的340个全球k-mer序列特征。IRESpy中的研究还表明，基于序列/结构混合特征的模型的性能比基于序列的模型略有提高。

然而，当前报道的研究仅限于使用手工制作的特征，这些特征几乎都不是最优的，结构在IRES预测中的作用也需要进一步探讨。与传统的机器学习算法相比，深度神经网络框架尚处于探索阶段，但在IRES预测方面可能很有前景。此外，目前的IRES预测方法是针对线性mRNA设计的。由于线性RNA和环状RNA的IRES活性存在差异，因此迫切需要发展环状RNAIRES预测法。

发明内容

针对现有技术中存在的缺陷，本发明提供了基于DeepCIP计算框架将特征嵌入与深度学习方法相结合实现预测环状RNA IRES的方法。本发明首次通过序列处理、RNA图构建、特征提取和模型集成实现了对环状RNA IRES的预测。在集成模块中，通过本发明所述方法开发的DeepIRESc工具采用软投票策略来集成由不同数据集训练的三个融合模型。每个融合模型包括四个模块，包括RNA处理模块、S-LSTM模块、GCN模块和特征融合模块。首先，RNA处理模块用于对输入的RNA序列进行预处理，以便进行序列编码、结构预测和RNA图构建。其次，利用S-LSTM模块和GCN模块分别提取RNA序列和RNA二级结构的特征。最后，在特征融合模块中，将从RNA序列和结构中提取的特征进行融合。融合特征通过完全连接层和乙状函数，最终输出一个概率，表示输入的RNA序列是一个环状RNAIRES的可能性。具体地，序列处理包括高置信度数据筛选、独立测试数据集构建，三个训练子集构建和RNA序列的表示方法选择四个步骤。

为了获取序列信息，本发明首次采用句子状态LSTM（S-LSTM）自动学习序列特征；并通过构建加权RNA图，利用图卷积网络（GCN）从RNA图中提取二级结构信息。通过特征融合模块，整合序列和结构特征，识别环状RNA IRES。为了证明通过本发明方法开发的DeepIRESc工具的有效性，本发明在定义的独立测试集上对DeepIRESc和XGBoost模型进行了基准测试。基准测试结果表明，DeepIRESc在环状RNA-IRES预测方面的性能优于其他比较方法。

本发明是这样实现的。

第一方面，本发明提供了序列处理的步骤。相应步骤如下：

第一步，筛选出用于训练模型的高置信度数据；第二步，将第一步中筛选出的数据集分为训练数据集和测试数据集；第三步，通过结合下采样和模型集成相结合的方法将第二步中的训练和测试集分割成三个训练子集。

由于Weingarten-Gabbay等人通过高通量双顺反子分析，从55000个寡核苷酸中鉴定了数千个具有IRES活性的序列<Do Zotero Refresh: [@weingarten-gabbaySystematic2016]>。在此基础上，Chen等人通过构建oligo-split-eGFP 环状RNA报告鉴定了17201个eGFP（+）寡核苷酸和23654个eGFP（-）寡核苷酸。其中， 1639个寡核苷酸的IRES具有线性特异性活性（即线性IRES），而4582个寡核苷酸的IRES具有环状特异性活性（即环状IRES）。由于高质量数据集对模型的预测性能至关重要，因此，在本发明中，预测环状RNA IRES方法中的序列处理的第一步就是要筛选出用于训练模型的高置信度数据。

进一步地，本发明从Weingarten-Gabbay等人的研究中选择环状IRES作为阳性样本，同时将排除了线性序列和IRES活性大于背景活性的序列（不包括启动子活性大于0.2和剪接活性小于-2.5的序列）后的eGFP（-）样本作为阴性样本。在筛选过程中，一共获得24525个RNA序列，其中4582个阳性样本和19943个阴性样本。由于这些数据包含了一部分用于测试特定突变对IRES活性影响的合成序列，本发明采用Weingarten-Gabbay等人的方法从数据集中仅保留了标记为“CDS_screen”、“Genome_Wide_Sceen_Elements”、“High_Priority_Genes_Blocks”、“Hight_Priurity_Virus_Block”、“Human_5UTR_screen”，“IRESite_blocks”、“Viral_5UTR_Screen”和“rRNA_Matching_5UTRs”的天然序列。

优选地，本发明筛选出了一个包含4531个阳性样本和9616个阴性样本的数据集。与此同时，基于用于IRES活性测定的寡核苷酸文库是人工构建的，故本发明中RNA序列的长度选为174nt。

由于缺环状IRES的独立测试数据集，本发明将上述所获得的数据集分为训练数据集和测试数据集用于模型的训练和评估。

进一步地，本发明从CircAtlas＜Do Zotero Refresh:[@wuCircAtlas2020]>收集人类环状RNA序列，并从Viruscircbase收集病毒环状RNA序列，然后使用Blastn将数据集中的RNA序列映射到数据库中的环状RNA。结果，获得与上述所获得的数据集具有100%同一性的582个阳性样本。

更进一步地，本发明随机选择582个阴性样本，与之前描述的阳性样本一起构建独立的测试集，用于本发明所构建模型的后续评估。

更进一步地，本发明将上述所获得的数据集中除开用于独立测试集的样本后，其中剩余的3949个阳性样本和9034个阴性样本将用作训练集。为了解决由于阳性样本和阴性样本数量不平衡导致模型偏向于具有更多样本的类从而降低模型泛化能力的问题，本发明采用一种结合下采样和模型集成的方法，通过阴性样本进行随机抽样，来获得与阳性样本数量相同的阴性样本的三个子集。随机抽样的要求是不浪费数据且数据在子集之间没有重复。然后，将所有阴性子集与阳性数据集组合，以获得具有3949个阳性样本和3949个阴性样本的三个训练子集。最后将使用三个训练子集训练的模型集成到一个最终模型中。

进一步地，本发明从序列信息和结构信息两个不同的角度来表示RNA序列。

优选地，本发明采用热编码将碱基a、T/U、C和G分别表示为[1,0,0,0], [0,1,0,0], [0,0,1,0]和 [0,0,0,1]。每个序列表示为维度为[L×4]的特征矩阵，其中L表示序列的长度，在本发明的数据集中L=174。

优选地，本发明采用ViennaRNA包（版本2.5.1）中的RNApfold来捕获动态RNA二级结构,同时通过McCaskill算法计算RNA的局部稳定二级结构，并输出RNA碱基配对的概率。在RNApfold运行过程中的参数设置为W=150，c=1e-3，并且不允许产生具有孤立对（--noLP）的结构。其他参数采用默认设置。

第二方面，本发明将第一方面预测获得的碱基配对概率用于构建RNA加权图G=（V，E，W）。

进一步地，所构建RNA加权图中每个基表示为图中的节点V，其中每个节点的特征类似地使用一个热向量进行编码。

进一步地，所构建RNA加权图中边缘E包含两种不同类型的化学键信息：沿着RNA骨架连接连续核苷酸的共价键和连接成对碱基的氢键，其中共价键的重量W定义为1，氢键的重量W等于其所连接的两个碱基配对的概率。

第三方面，本发明提供了RNA序列特征提取的方法。

进一步地，用于RNA序列特征提取的模型结构是句子状态LSTM（S-LSTM）。S-LSTM是一种替代的递归神经网络，旨在解决BiLSTM的局限性。BiLSTM是递归神经网络的一种变体，由前向和后向LSTM组成，广泛应用于自然语言处理。S-LSTM使用与BiLSTM类似的门机制来控制信息流。它们的不同之处在于，S-LSTM将每个单词视为单词级节点w，并添加句子级节点g以表示完整句子。尽管S-LSTM作用于序列，但节点聚合和消息传递的方式类似于GNN的操作，可以将其归类为文本GNN。即在每个时间步骤t，可以在句子级节点和每个词级节点之间交换信息，并且每个词节点也将与其上下文节点交换信息。在每个时间步t的S-LSTM状态可以定义为公式-1：，其中 是每个字的子状态，是句子的子状态。

随着时间步长t的增加，每个h_i捕获越来越大的n元语法，同时与g交换信息，使得和g学习的上下文信息越来越丰富。最终g可用于任务分类。在默认情况下，每个词节点仅与相邻词节点交换信息，其窗口大小为1。但增加窗口大小可以允许更多信息交流。本发明中，窗口大小为3，时间步长为7。S-LSTM中词状态和句子状态的更新过程由以下公式-2到公式-7给出：

公式- 2；

公式- 3 ；

公式-4；

公式-5；

公式-6；

公式-7；

其中，是上下文窗口的串联向量，表示不同的门。和控制信息分别来自输入, 左侧上下文单元格,右侧上下文单元格, 和句子上下文单元格。为输出门。为实际输入。 W, , U, V和  b为可训练参数。σ是sigmoid函数。

与此同时，模型中S-LSTM句子状态的更新过程可以通过以下公式-8到公式-14来描述：

公式- 8；

公式- 9；

公式- 10 ；

公式- 11；

公式- 12；

公式- 13；

公式- 14 ；

其中， 和 为标准化的门用于控制来自和。 为输出门。 W,U, 和  b为可训练参数。

进一步地，每个RNA二级结构采用加权RNA图表示。

优选地，本发明通过堆叠多个GCN层进行核苷酸节点的特征向量学习。为了便于对RNA序列进行分类，本发明将每个RNA图中的核苷酸节点特征聚集到图嵌入中。全局和、最大和平均池是聚合节点功能的最常用策略。即：节点特征可以用矩阵表示，节点之间的连通性可以用具有权重的邻接矩阵表示，其中N是节点数，D是每个节点特征向量的维数。矩阵X和A作为模型输入，而Kipf和Welling中的GCN分层传播规则由公式-15给出（Semi-Supervised Classification with Graph Convolutional Networks.2017）：

公式- 15其中，

 和是单位矩阵。是的对角度矩阵。是激活矩阵， 是层中 的可训练矩阵， 是非线性激活函数。

因此，在几个GCN层之后，矩阵X被转换为公式- 16：（F是滤波器的数量），每个节点 可以用公式-17表示。

公式- 16；

公式- 17；

这里， 表示滤波器参数矩阵，其中，, 是从节点j到节点i的边缘权重。

本发明选择通过三个GCN层堆叠来学习节点级嵌入，因为GCN层太少会导致接收场较小，而太多的GCN层可能导致过度平滑。此外，在本发明中，在最后的GCN层之后添加了一个全局和池层提取图级以表示任务分类，因为这样的处理具有相对更好的性能。

为了结合S-LSTM和GCN模块中的信息序列特征和结构特征，本发明在最后一个维度中将它们连接起来。然后，使用带有sigmoid函数的分类器输出预测概率。

所有训练子集分别用于训练，得到三个模型。本发明采用软投票策略来整合这三个模型。投票机制是最典型的模型集成方法，其基本思想是输出所有分类器中投票最多的类。分类器直接给出称为硬投票的最终预测标签，并给出称为软投票的输出标签的预测概率。本发明使用软投票来计算三个模型概率的加权和，然后根据默认阈值0.5将预测标签确定为0或1。优选地，本发明将三个模型的权重设置为相等。

第四方面，本发明使用嵌套交叉验证（CV）进行超参数调整和模型选择。Errica等人使用嵌套交叉验证（CV）比较了不同的GNN模型。简单地说，嵌套CV包含外部CV和内部CV，其中CV可以选择k折叠或保持技术。K倍CV表示将数据集随机划分为  非重叠子集，其中每个子集轮流作为测试集，其余子集用于训练。保持CV意味着将数据集直接划分为两个相互排斥的子集，一个用于训练，另一个用于验证。在本发明中，所述验证采用前面所述三个训练子集中的任一个来执行。

进一步地，本发明对外部CV使用k倍技术且，对内部CV使用90%的训练分割和10%的验证分割。即，本发明训练每个外部训练折叠，并将10%随机数据作为验证集，如果在执行n个周期后没有改善性能，执行就提前停止，然后在测试折叠中进行测试。最终超参数评估分数为所有测试折叠分数的平均值。

进一步地，在完成超参数选择后，本发明使用10倍交叉验证在所有数据上重新训练和验证模型，最后在本发明所构建的独立测试集上评估模型性能。为了评估所提出模型及其各个子模块的有效性，本发明将其与使用嵌套CV的其他基线方法进行了比较。不同基线方法的描述如下：

序列的TextCNN：TextCNN是用于文本分类的卷积神经网络，它由嵌入层、卷积层、最大池层和完全连接层组成。其中序列被编码为嵌入层中的一个热向量，并且在卷积层中设置n个filters=64和filters_size=[2,3,4]。

BiLSTM：BiLSTM用于从序列中学习有效特征，并且可以多层叠加。单层、双层和三层BiLSTM用于模型比较，其中hidden_size设置为64。

TextCNN+注释二级结构：在这里，RNAfold（Mol. Biol. AMB 2011; 6:26）用于预测RNA二级结构，并输出表示碱基配对模式的点括号序列。bpRNA（Nucleic Acids Res.2018; 46:5381–5394）对预测的RNA二级结构进行注释，将RNA折叠输出的碱基配对信息解析为一个详细的结构，提供相关的上下文注释信息，包括茎（S）、发夹环（H）、多环（M）、内部环（I）、凸起（B）和末端（E）。这允许RNA二级结构表示为序列，并输入TextCNN模型进行训练。上一节介绍了TextCNN的设置。

图的多层感知器（MLP）：该方法在RNA图的节点特征上应用具有ReLU激活的三层MLP，然后使用全局和池层来学习图级嵌入。它与GCN模块的不同之处在于，它在学习节点特征时不使用节点之间的关系，并且是一个与图形拓扑无关的模型。

XGBoost：该模型用于开发IRESpy工具（BMC Bioinformatics 2019; 20:409）。本发明使用340个序列k-mer（1-mer、2-mer、3-mer和4-mer）特征来训练XGBoost模型，并设置与IRESpy相同的XGBoost 超参数。scale_pos_weight参数在XGBoost模型中提供，通过控制正负权重的平衡，对于不平衡类很有用。本发明使用环状RNA IRES训练数据集训练了两个XGBoost模型，对于circ称为XGBoost_weight，对于cirk称为XGBoost，将scale_pos_weight参数分别设置为3，scale_pos _weight设置为1（与IRESpy相同）。本发明使用10倍CV和提前停止方法来获得最佳num_boost_round参数，并选择最佳模型与最终模型进行比较。

进一步地，本发明采用五个评估指标，包括准确性（Acc）、敏感性（Sn）、特异性（Sp）、精确度和马修斯相关系数（MCC）。指标定义如如下公式- 18到公式-22所描述：

公式- 18；

公式- 19；

公式- 20；

公式- 21；

公式- 22 ；

其中TP、TN、FP和FN表示真阳性、真阴性、假阳性和假阴性的数量。

进一步地，本发明使用ROC（接收器操作特征）曲线和P-R（精确召回）曲线直观地评估分类模型的整体性能。

附图说明

图1是本发明的流程框图。

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。图1是本发明的流程框图。

本实施例提供的一种基于DeepCIP计算框架将RNN和GNN的多模式深度学习方法用于提取RNA序列和结构的特征，并将这两种特征结合起来用于环状RNA-IRES预测过程，包括序列处理、RNA图构建、特征提取和模型集成四个阶段。具体地，在集成模块中，利用通过本发明方法开发的DeepIRESc工具采用软投票策略来集成由不同数据集训练的三个融合模型。每个融合模型包括四个模块，包括RNA处理模块、S-LSTM模块、GCN模块和特征融合模块。首先，RNA处理模块用于对输入的RNA序列进行预处理，以便进行序列编码、结构预测和RNA图构建。其次，利用S-LSTM模块和GCN模块分别提取RNA序列和RNA二级结构的特征。最后，在特征融合模块中，将从RNA序列和结构中提取的特征进行融合。融合特征通过完全连接层和乙状函数，最终输出一个概率，表示输入的RNA序列是一个环状RNA IRES的可能性。更具体地，序列处理包括高置信度数据筛选、独立测试数据集构建，三个训练子集构建和RNA序列的表示方法选择四个步骤。

针对序列处理：

第一步，筛选出用于训练模型的高置信度数据；

 首先，从Weingarten-Gabbay等人鉴定的具有IRES活性的序列<Do ZoteroRefresh: [@weingarten-gabbaySystematic2016]> ，选择环状IRES作为阳性样本，同时将排除了线性序列和IRES活性大于背景活性的序列（不包括启动子活性大于0.2和剪接活性小于-2.5的序列）后的eGFP（-）样本作为阴性样本。结果，获得24525个RNA序列，其中4582个阳性样本和19943个阴性样本；

为了去掉所筛选样本中收集的用于测试特定突变对IRES活性影响的合成序列，进一步采用Weingarten-Gabbay等人的方法对数据集进行优选，仅保留标记为“CDS_screen”、“Genome_Wide_Sceen_Elements”、“High_Priority_Genes_Blocks”、“Hight_Priurity_Virus_Block”、“Human_5UTR_screen”，“IRESite_blocks”、“Viral_5UTR_Screen”和“rRNA_Matching_5UTRs” 的天然序列；

结果，筛选出了筛选出用于训练模型的一个包含4531个阳性样本和9616个阴性样本的高置信度数据集。鉴于IRES活性测定的寡核苷酸文库是人工构建的，故RNA序列的长度选为174nt。

第二步，为了构建环状IRES的独立测试数据集，将第一步所获得的数据集分为训练数据集和测试数据集；

首先，从CircAtlas＜Do Zotero Refresh:[@wuCircAtlas2020]>收集人类环状RNA序列，并从Viruscircbase收集病毒环状RNA序列，然后使用Blastn将数据集中的RNA序列映射到数据库中的环状RNA。结果，获得与上述所获得的数据集具有100%同一性的582个阳性样本。随后，随机选择582个阴性样本，与所获得的阳性样本一起构建独立的测试集，用于所构建模型的后续评估；

进一步地，将第一步所获得的数据集中除开用于独立测试集的582个阳性样本和582个阴性样本后，其中剩余的3949个阳性样本和9034个阴性样本将用作训练集。

第三步，通过结合下采样和模型集成相结合的方法将第二步中的训练和测试集分割成三个训练子集；

为了解决由于阳性样本和阴性样本数量不平衡导致模型偏向于具有更多样本的类从而降低模型泛化能力的问题，采用一种结合下采样和模型集成的方法，对阴性样本进行随机抽样，来获得与阳性样本数量相同的阴性样本的三个子集。随机抽样的要求是不浪费数据且数据在子集之间没有重复；

然后，将所有阴性子集与阳性数据集组合，以获得具有3949个阳性样本和3949个阴性样本的三个训练子集。

第四步，RNA序列的表示方法选择：

首先，采用序列信息表示RNA序列。具体地，通过热编码将碱基a、T/U、C和G分别表示为[1,0,0,0], [0,1,0,0], [0,0,1,0]和 [0,0,0,1]，每个序列表示为维度为[L×4]的特征矩阵，其中L表示序列的长度；在本发明L=174;

其次，采用结构信息表示RNA序列。具体地，采用ViennaRNA包（版本2.5.1）中的RNApfold（Bioinformatics 2006, 22:614–615）来捕获动态RNA二级结构,同时通过McCaskill算法(Biopolymers 1990; 29:1105–1119)计算RNA的局部稳定二级结构，并输出RNA碱基配对的概率。在RNApfold运行过程中的参数设置为W=150，c=1e-3，并且不允许产生具有孤立对（--noLP）的结构。其他参数采用默认设置。

针对RNA图构建：

将序列处理阶段第四步获得的碱基配对概率用于构建RNA加权图G=（V，E，W）：所构建RNA加权图中每个基表示为图中的节点V，其中每个节点的特征类似地使用一个热向量进行编码；所构建RNA加权图中边缘E包含两种不同类型的化学键信息：沿着RNA骨架连接连续核苷酸的共价键和连接成对碱基的氢键，其中共价键的重量W定义为1，氢键的重量W等于其所连接的两个碱基配对的概率。

针对特征提取：

首先，用于RNA序列特征提取的模型结构采用句子状态LSTM（S-LSTM）。在每个时间步t的S-LSTM状态定义为：

公式-1： ，其中 是每个字的子状态，是句子的子状态。

随着时间步长t的增加，每个h_i捕获越来越大的n元语法，同时与g交换信息，使得和g学习的上下文信息越来越丰富。最终g可用于任务分类。在默认情况下，每个词节点仅与相邻词节点交换信息，其窗口大小为1。但增加窗口大小可以允许更多信息交流。

序列模型的性能受到许多超参数的影响，例如句子级节点的数量、窗口大小、时间步长和S-LSTM模块中的隐藏层大小。本发明在设计实验过程中重点关注S-LSTM模型中窗口大小和时间步长的影响，其中窗口大小从{1、2、3}不等，时间步长从{5、7、9}不等。通过比较分析各种S-LSTM设置的实验结果，发现将窗口大小设置为3，时间步长设置为9，可获得平均AUC的最佳值。对在不同层数下的S-LSTM、TextCNN和BiLSTM的性能进行比较，发现S-LSTM在平均AUC值方面优于其他模型。上述实验结果表明S-LSTM可以更有效地提取环状RNA IRES序列特征。在本发明中，窗口大小为3，时间步长为9。

与此同时，S-LSTM中词状态和句子状态的更新过程由以下公式-2到公式-7给出：

公式- 2；

公式- 3 ；

公式-4；

公式-5；

公式-6；

公式-7；

其中，是上下文窗口的串联向量，表示不同的门。和控制信息分别来自输入, 左侧上下文单元格,右侧上下文单元格, 和句子上下文单元格。为输出门。为实际输入。 W, , U, V和  b为可训练参数。σ是sigmoid函数。

与此同时，模型中S-LSTM句子状态的更新过程可以通过以下公式-8到公式-14来描述：

公式- 8；

公式- 9；

公式- 10 ；

公式- 11；

公式- 12；

公式- 13；

公式- 14 ；

其中， 和 为标准化的门用于控制来自和。 为输出门。 W,U, 和  b为可训练参数。

其次，每个RNA二级结构采用加权RNA图表示。通过堆叠多个GCN层进行核苷酸节点的特征向量学习。为了便于对RNA序列进行分类，将每个RNA图中的核苷酸节点特征聚集到图嵌入中。节点特征用矩阵表示，节点之间的连通性用具有权重的邻接矩阵 表示，其中N是节点数，D是每个节点特征向量的维数。矩阵X和A作为模型输入，而Kipf和Welling中的GCN分层传播规则由公式-15给出（Semi-SupervisedClassification with Graph Convolutional Networks. 2017）：

公式- 15其中，

 和是单位矩阵。是的对角度矩阵。是激活矩阵， 是层中 的可训练矩阵， 是非线性激活函数。

因此，在几个GCN层之后，矩阵X被转换为公式- 16：（F是滤波器的数量），每个节点 可以用公式-17表示。

公式- 16；

公式- 17；

这里， 表示滤波器参数矩阵，其中，, 是从节点j到节点i的边缘权重。

GCN层数和读出功能对从图形中提取RNA特征至关重要。在这里，本发明在实验设计过程中通过改变GCN层的数量和选择不同的读出函数来构建GCN模块，以研究不同配置对模型性能的影响。其中，GCN层数从{1、2、3}不等，读取策略从{全局和池、全局最大池、全局平均池}中选择。实验结果表明，在AUC得分方面，GCN的三层和全局和池层获得了最佳性能。

所以，本发明随后选择通过三个GCN层堆叠来学习节点级嵌入，因为GCN层太少会导致接收场较小，而太多的GCN层可能导致过度平滑。此外，在最后的GCN层之后添加了一个全局和池层提取图级以表示任务分类，因为这样的处理具有相对更好的性能。

为了研究在本发明的模型中使用加权RNA图是否可以提高环状RNA-IRES结构特征的表示能力，本实施例比较了加权RNA图+GCN、TextCNN+注释二级结构（TextCNN_nannotate_struct）和图的MLP的性能。比较结果表明GCN模块达到了最佳性能。上述这一切表明加权RNA图比顺序注释结构能更好地表示RNA的二级结构（其中，GCN比TextCNN_nannotate_struct性能更好）。与此同时，使用核苷酸之间不同的化学键关系是一个更好的选择，结果表明GCN比MLP在图形方面表现更好。

另，为了结合S-LSTM和GCN模块中的信息序列特征和结构特征，本发明在最后一个维度中将它们连接起来。然后，使用带有sigmoid函数的分类器输出预测概率。

针对模型集成：

首先，且，对内部CV使用90%的训练分割和10%的验证分割。即，本发明训练每个外部训练折叠，并将10%随机数据作为验证集，如果在执行n个周期后没有改善性能，执行就提前停止，然后在测试折叠中进行测试。最终超参数评估分数为所有测试折叠分数的平均值。在完成超参数选择后，使用10倍交叉验证在所有数据上重新训练和验证模型，最后在所构建的独立测试集上评估模型性能。

为了评估所提出模型及其各个子模块的有效性，本发明将其与使用嵌套CV的其他基线方法进行了比较。不同基线方法的描述如下：

序列的TextCNN：TextCNN是用于文本分类的卷积神经网络，它由嵌入层、卷积层、最大池层和完全连接层组成。其中序列被编码为嵌入层中的一个热向量，并且在卷积层中设置n个filters=64和filters_size=[2,3,4]。

BiLSTM：BiLSTM用于从序列中学习有效特征，并且可以多层叠加。单层、双层和三层BiLSTM用于模型比较，其中hidden_size设置为64。

TextCNN+注释二级结构：在这里，RNAfold（Mol. Biol. AMB 2011; 6:26）用于预测RNA二级结构，并输出表示碱基配对模式的点括号序列。bpRNA（Nucleic Acids Res.2018; 46:5381–5394）对预测的RNA二级结构进行注释，将RNA折叠输出的碱基配对信息解析为一个详细的结构，提供相关的上下文注释信息，包括茎（S）、发夹环（H）、多环（M）、内部环（I）、凸起（B）和末端（E）。这允许RNA二级结构表示为序列，并输入TextCNN模型进行训练。前面已介绍了TextCNN的设置。

图的多层感知器（MLP）：该方法在RNA图的节点特征上应用具有ReLU激活的三层MLP，然后使用全局和池层来学习图级嵌入。它与GCN模块的不同之处在于，它在学习节点特征时不使用节点之间的关系，并且是一个与图形拓扑无关的模型。

XGBoost：该模型用于开发IRESpy工具（BMC Bioinformatics 2019; 20:409）。本发明使用340个序列k-mer（1-mer、2-mer、3-mer和4-mer）特征来训练XGBoost模型，并设置与IRESpy相同的XGBoost超参数。scale_pos_weight参数在XGBoost模型中提供，通过控制正负权重的平衡，对于不平衡类很有用。本发明使用环状RNA IRES训练数据集训练了两个XGBoost模型，对于circ称为XGBoost_weight，对于cirk称为XGBoost，将scale_pos_weight参数分别设置为3，scale_pos _weight设置为1（与IRESpy相同）。本发明使用10倍CV和提前停止方法来获得最佳num_boost_round参数，并选择最佳模型与最终模型进行比较。

具体地，当模型训练完毕后，固定模型参数，以预测值和真实值之间的均方误差作为优选目标，重复五次，且每次用做验证集的子集不同。在训练过程中，如果当前模型在验证集上的准确度或均方误差优于之前的模型，则保存并替换之前的模型，该过程重复100次。

训练完模型后，对获得模型进行测试，具体地，利用测试集对模型进行测试。受试者特性曲线下面积(AUC)是模型参数选择的主要依据。使用测试集对目前文献报道的XGBOOST模型进行准确性和精确度比较，发现本发明所构建模型的结果要优于现有的方法。

采用五个评估指标，包括准确性（Acc）、敏感性（Sn）、特异性（Sp）、精确度和马修斯相关系数（MCC）。指标定义如如下公式- 18到公式-22所描述：

公式- 18；

公式- 19；

公式- 20；

公式- 21；

公式- 22 ；

其中TP、TN、FP和FN表示真阳性、真阴性、假阳性和假阴性的数量。

使用ROC（接收器操作特征）曲线和P-R（精确召回）曲线直观地评估分类模型的整体性能;

所有训练子集分别用于训练，得到三个模型。本发明采用软投票策略来整合这三个模型。投票机制是最典型的模型集成方法，其基本思想是输出所有分类器中投票最多的类。分类器直接给出称为硬投票的最终预测标签，并给出称为软投票的输出标签的预测概率。本发明使用软投票来计算三个模型概率的加权和，然后根据默认阈值0.5将预测标签确定为0或1。优选地，本发明将三个模型的权重设置为相等。

此外，本实施例表明特征融合和软投票机制能提高利用本发明方法开发的DeepIRESc工具对环状RNA IRES的预测性能。为了考察RNA结构在环状RNA IRES识别中的重要性，将融合模型与仅基于序列特征的模型在独立测试集上进行了比较。结果表明，AUROC、AUPRC、ACC、SN、Precision和MCC的性能均得到了提高。这也说明结构特征在环状RNA IRES预测中非常重要。

随后，本实施例对独立测试数据集进行了测试，以验证软投票机制是否可以提高预测性能。在本发明，软投票方法被用于由三个训练子集训练的三个集成模型。从三个单一模型与集成模型预测性能的比较结果表明，集成模型具有最佳性能，其AUROC分别为0.723，AUPRC为0.697。由于用于集成的模型很少，以及训练子集之间存在部分重叠，导致集成模型预测性能的提高并不是特别地显著。

进一步地，为了评估本发明模型的预测性能，在本实施例中，首先使用独立的测试数据集对利用本发明方法开发的工具和XGBoost模型进行了实验研究比较。比较结果表明，利用本发明方法开发的工具比XGBoost_weight模型具有更高的性能，其中AUC为0.723（vs.0.624），ACC为0.646（vs. 0.624），SP为0.512（vs. 0.503），SN为0.780（vs. 0.744），精度为0.615（vs. 0.600），MCC为0.303（vs. 0.255）。这归因于本发明的模型更好地代表了RNA结构。此外，与手工绘制的k-mer特征相比，通过深度学习自动提取的特征具有更好的表示能力。上述这些进一步表明深度学习可能是解决IRES共同特征尚未完全明确的现状的一个有益选择。

为了进一步确定利用本发明方法开发的DeepIRESc工具预测环状RNA IRES的有效性，本实施例探索了预测的环状RNA IRES概率与环状RNA IR实验活动之间的相关性。环状RNA IRES活性由Chen等人实验中的eGFP表达定义。在独立测试集中，环状RNA IRES的活性在0到6.9327之间，其中背景eGFP表达为3.466387。结果表明，具有较高活性的环状RNAIRES通常也具有较高的预测概率。

此外，本发明从独立测试集可视化了利用本发明方法开发的DeepIRESc工具对人类序列和病毒序列的预测概率。一般来说，人类和病毒环状RNA IRES的预测概率都高于非环状RNA的IRES。这进一步证明了利用本发明方法开发的DeepIRESc工具对人类和病毒环状RNA IRES预测具有较好的能力。

针对预测应用：

从文献中收集了不同长度的经实验验证的环状RNA IRES序列，因为阳性样本组成了一个基本事实数据集，其中包含10个智人环状RNA的IRES序列和4个黑腹果蝇环状RNA中的IRES顺序。真实有效数据的详细信息如表1所示。

表1

,

与此同时，将利用本发明方法开发的DeepIRESc工具预测结果与其他模型（包括XGBoost_weight、XGBoost和IRESpy）进行了比较。预测结果的比较如表2所示。对于本发明模型来说，如果预测阈值设置为0.5，那么智人的10个RNA序列中有6个被预测为环状RNAIRES，果蝇的4个RNA序列当中有1个被预测是环状RNA-IRES，这比其他方法的模型表现好。其中，利用本发明方法开发的DeepIRESc工具对黑腹果蝇RNA序列的不良表现可能是由于不同物种的环状RNA-IRES特征存在差异。因为本发明所采用的训练数据集仅包含人类和病毒的RNA序列。

表2

,

此外，通过核密度估计（KDE）比较了使用不同方法对人类环状 IRES进行预测概率分布的性能。结果表明，利用本发明方法开发的DeepIRESc工具比XGBoost的预测性能要更好。进一步地，IRESpy（用于线性RNA IRES识别）的预测性能低于实际环状RNA IRES数据中用于环状RNA IRES预测的其他方法。这进一步证实本发明开发高效预测环状RNA IRES方法的重要意义。

由此可见，利用本发明提供的方法，结合训练好的模型，能够快速较准确地预测具有介导翻译功能的环状RNA IRES序列。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种预测环状RNA IRES的方法，其特征在于，所述方法是依据所提取RNA的序列和结构特征结合序列处理、RNA图构建、特征提取和模型集成四个阶段来进行环状RNA IRES预测的，所述RNA的序列和结构特征是基于DeepCIP计算框架结合RNN和GNN多模式深度学习方法来进行提取的，所述RNA的序列和结构特征的提取工具是DeepIRESc工具，所述DeepIRESc工具采用软投票策略来集成由不同数据集训练的三个融合模型，所述融合模型包括RNA处理模块、S-LSTM模块、GCN模块和特征融合模块，所述RNA处理模块用于预处理输入的RNA序列，所述S-LSTM模块用于提取RNA序列特征，所述GCN模块用于提取RNA二级结构特征，所述特征融合模块用于融合所述提取的RNA序列和RNA二级结构特征，所述融合是通过完全连接层和乙状函数输出一个概率，所述概率是输入的RNA序列为环状RNA IRES的可能性。

2.根据权利要求1所述的一种预测环状RNA IRES的方法，其特征在于，所述预测环状RNA IRES方法中的序列处理阶段包括以下步骤：第一步，筛选出用于训练模型的高置信度数据；第二步，将第一步中筛选出的数据集分为训练数据集和测试数据集；第三步，通过结合下采样和模型集成相结合的方法将第二步中的训练和测试集拆分成三个训练子集。

3.根据权利要求1所述的一种预测环状RNA IRES的方法，其特征在于，所述预测环状RNA IRES方法中的序列处理包括采用序列信息和结构信息表示RNA序列的表示方法；所述序列信息的表示方法为：采用热编码将碱基a、T/U、C和G分别表示为[1,0,0,0], [0,1,0,0], [0,0,1,0]和 [0,0,0,1]，每个序列表示为维度为[L×4]的特征矩阵，其中L表示序列的长度；所述结构信息的表示方法为：采用RNA折叠来捕获动态RNA二级结构，通过计算RNA的局部稳定二级结构输出RNA碱基配对的概率。

4.根据权利要求1所述的一种预测环状RNA IRES的方法，其特征在于，所述预测环状RNA IRES方法中特征提取阶段模型结构的句子状态为S-LSTM，所述模型结构的节点聚合和消息传递的方式为文本GNN；所述模型结构的核苷酸节点的特征向量通过堆叠多个GCN层来学习。

5.根据权利要求1所述的一种预测环状RNA IRES的方法，其特征在于，所述预测环状RNA IRES方法中的模型集成阶段采用嵌套CV方法进行超参数调整和模型选择，所述的超参数调整包括对外部CV使用k倍技术，对内部CV使用90%的训练分割和10%的验证过程，所述的模型选择包括采用10倍交叉验证在所有数据上重新训练，验证模型和在独立测试集上评估模型性能的过程。

Images