CN111161793B - 基于stacking集成的RNA中N6-甲基腺苷修饰位点预测方法 - Google Patents

Abstract

本发明公开了基于stacking集成的RNA中N⁶‑甲基腺苷修饰位点预测方法，涉及系统生物学领域。通过6种特征提取方法提取酿酒酵母、智人和拟南芥3个物种的RNA序列特征，通过特征融合得到原始数据集的初始特征空间；利用弹性网络对其进行降维，剔除冗余、噪声特征，保留对模型分类相关的重要特征，得到最佳特征集合；把最优特征子集以及所对应的类别标签输入到stacking集成中进行模型训练，并结合评价指标评估模型的预测性能，得到预测模型；将测试集中待预测的RNA序列，输入至预测模型中，预测m⁶A位点并输出。此模型在测试集上的预测准确率分别达到92.30％和87.06％，在跨物种预测方面具有很好的发展潜力，可成为鉴定m⁶A位点的有用工具。

Description

基于stacking集成的RNA中N6-甲基腺苷修饰位点预测方法

技术领域

本发明属于系统生物学领域，涉及一种基于stacking集成的RNA中N⁶-甲基腺苷修饰位点预测方法。

背景技术

人类基因组计划的完成大大促进人们对遗传组织信息、传递和表达规律的认识，同时使我们意识到细胞内遗传信息表达机制的异常复杂性。RNA作为中心法则的关键环节，将遗传物质DNA和生命活动的执行者蛋白质紧密相连。研究表明，RNA中存在100多种化学修饰，这些化学修饰以甲基化修饰为主，包括：5-甲基胞嘧啶(m⁵C)、N¹-甲基腺嘌呤(m¹A)、5-羟甲基胞嘧啶(hm⁵C)、假尿嘧啶(ψ)，N⁶-甲基腺嘌呤(m⁶A)等。在20世纪70年代发现的N⁶-甲基腺苷(m⁶A)是一种常见的转录后RNA修饰，随后的研究陆续证实了该修饰广泛存在于动植物、细菌以及病毒在内的多个物种中。

作为一个动态和可逆的过程，m⁶A发生在腺嘌呤的第六个氮原子上，其动态变化可以通过调节多种与RNA相关的细胞信号通路影响基因表达和细胞命运，在mRNA剪接、输出、稳定、免疫耐受、RNA转录、加工、细胞分裂以及细胞分化等生物过程中发挥重要作用。此外，m⁶A修饰与人类疾病密切相关，包括：癌症、病毒感染和脑部发育异常。因此，准确识别m⁶A位点对于RNA甲基化修饰基础研究，对了解疾病机制和药物开发都至关重要，而且具有推动生物信息学发展的作用。用于鉴定RNA序列中m⁶A位点的方法有：二维薄层色谱，高效液相色谱和高通量方法(如m⁶A-seq和MeRIP-Seq)。然而，用于靶向m⁶A位点的纯生化实验方法非常昂贵、耗时，随着先进测序技术和基因组计划的发展，积累了大量的RNA序列，许多研究人员提出基于机器学习算法的有效计算方法，用于快速准确的预测m⁶A位点。

目前为止，一系列基于机器学习的m⁶A位点计算方法被提出。Huang等人通过集成深度学习算法和随机森林提出一种跨物种分类器BERMP来预测m⁶A位点。Zhao等人提出模型HMpre，用成本敏感的学习思想来解决人类mRNA中m⁶A位点预测问题中的不平衡数据问题，独立的测试数据集精度、F1和MCC的值分别达到0.3035，0.3961，0.3329。Chen等人提出iRNA-PseDNC模型，使用假二核苷酸组成鉴定N⁶-甲基腺苷位点，10折交叉验证表明iRNA-PseDNC的预测性能优于RAM-NPPS。Chen等人提出预测模型RAM-ESVM，使用集成支持向量机预测RNA转录组中的N⁶-甲基腺苷位点，刀切法测试结果表明RAM-ESVM优于单支持向量机分类器。Wang等人提出预测拟南芥中m⁶A位点的新工具RFAthM6A。Akbar等人基于SVM提出计模型iMethyl-STTNC，该模型通过将SAAC的概念融合到Chou的PseAAC编码RNA序列来鉴定N⁶-甲基腺苷位点，数据集1和数据集2基于STTNC特征的预测准确率分别达到69.84％，91.84％。Zhang等人提出一种计算方法鉴定大肠杆菌基因组中的RNAN⁶-甲基腺苷位点，10折交叉验证和独立测试集上的预测准确率均大于90％。Zhang等人引入一种新的启发式核苷酸物理化学性质选择(HPCS)算法构建M6A-HPCS预测器，通过启发式选择核苷酸物理化学性质提高N⁶-甲基腺苷位点预测性能。Xiang等人基于支持向量机构建模型RNAMethPre预测mRNA中的m⁶A位点。Qiang等人基于eXtreme Gradient Boosting(XGBoost)提出预测模型M6AMRFS，使用二进制编和局部位置特异性二核苷酸频率编码序列。Chen等人基于支持向量机，提出鉴定N⁶-甲基腺苷位点的预测工具MethyRNA，H.sapiens和M.musculus的预测准确率分别达到90.38％和88.39％。

尽管研究者在RNA甲基化修饰做出很大的贡献，提出了众多模型预测m⁶A位点，但预测模型多数是基于单个传统分类器或者受限于单一物种。因此，有必要设计一种新的跨物种m⁶A位点预测方法。

发明内容

本发明的目的在于提供一种基于stacking集成的RNA中N⁶-甲基腺苷修饰位点预测方法，实现跨物种预测，解决上述问题，该预测方法易于实施，预测准确率高。

为解决上述的技术问题，本发明采用以下技术方案一种基于stacking集成的RNA中N⁶-甲基腺苷修饰位点预测方法，其特征在于包括以下步骤：

1)收集RNA中N⁶-甲基腺苷修饰(m⁶A)位点信息：获得3个不同物种RNA的N⁶-甲基腺苷修饰位点数据集，包括正负数据集的RNA样本序列以及所对应的类别标签，确定训练集与测试集。

2)特征编码：使用序列衍生信息编码RNA序列，将字符信息转化成机器学习可识别的数值向量，融合6种特征提取方法对应的不同类型信息的特征向量，得到原始数据集的初始特征空间。

3)特征选择：对初始特征空间，利用弹性网络Elastic Net对其进行降维，剔除冗余、噪声特征，保留对模型分类相关的重要特征，得到最佳特征集合。

4)预测模型算法：把弹性网络Elastic Net对应的最优特征子集以及所对应的类别标签输入到stacking集成中进行模型训练，并结合评价指标评估模型的预测性能。

5)m⁶A位点预测：将独立测试集中待预测的RNA序列，输入至步骤4)中的预测模型中，预测该RNA序列是否包含m⁶A位点并输出。

更进一步的技术方案是所述步骤1)包括如下子步骤：

1-1)选择S.cerevisiae(酿酒酵母)、H.sapiens(智人)和A.thaliana(拟南芥)三个物种真实的数据；数据集S.cerevisiae作为训练集，数据集H.sapiens和A.thaliana作为测试集。

1-2)训练集和测试集中均包含N⁶-甲基腺苷修饰位点的RNA(正样本)和非N⁶-甲基腺苷修饰位点的RNA(负样本)。

1-3)S.cerevisiae数据集样本序列是以共同基序GAC为中心，窗口大小为51。如果样本序列在某些位置不存在核苷酸，缺少核苷酸使用其对称位置的核苷酸填充。数据集中包括1307条正样本，为了避免不平衡数据集对构建稳健模型的影响，从包含33,280条负样本的数据集中随机选取1307条负样本。

1-4)H.sapiens数据集样本窗口大小为41，若序列样本实际长度小于窗口大小，则缺少的核苷酸用其最临近的相同核苷酸填充，为了避免生成高度倾斜的数据集，该测试集中正负数据集数目保持一致，均为1130；A.thaliana数据集使用CD-HIT程序剔除序列相似性高于60％的样本，减少序列同源性偏差，数据集中包含394条正样本和394条负样本，样本窗口大小为25。

更进一步的技术方案是所述步骤2)6种特征编码方法包括二进制编码(binaryencoding)，化学性质(chemical property)，核苷酸频率(nucleotide frequency)，k-mer核苷酸频率(k-mer nucleotide frequency)，伪二核苷酸组成(pseudo dinucleotidecomposition,PseDNC)和位置特异性三核苷酸倾向(position-specifictrinucleotidepropensity,PSTNP)特征提取方法，其中：

所述binary encoding特征提取方法将能够准确描述样本序列中每个位置的核苷酸，该特征提取方法将RNA序列中包含的四种核苷酸：腺嘌呤(A)，鸟嘌呤(G)，胞嘧啶(C)和尿嘧啶(U)，依次编码成4维二进制向量。

所述chemicalproperty特征提取方法考虑到四种核苷酸具有不同的化学性质，依据环数、化学功能以及二级结构将A,C,G,U划分为三种不同类型的组，每种核苷酸用三维向量表示。

所述nucleotide frequency特征提取方法计算RNA序列中特定位置核苷酸的密度，反映m⁶A位点附近核苷酸的组成和频率；对于长度为L的样本，将生成L-1维的特征向量。

所述k-mer nucleotide frequency特征提取方法基于相邻核苷酸对会影响RNA序列的结构和功能的原理，计算相邻核苷酸在样本序列中出现的频率，反映m⁶A位点和非m⁶A位点序列背景之间的差异性，每条样本序列生成16维的特征向量。

所述pseudo dinucleotide composition(PseDNC)特征提取方法考虑到RNA序列的局部和全局序列信息，生成16+λ维的特征向量，前16个元素反映了RNA序列中的短程或者局部序列顺序信息，而剩余的元素表示远程或者全局的序列顺序信息，本文确定PseDNC中的最佳参数λ为23，对于每条样本序列生成39维的特征向量。

所述position-specifictrinucleotidepropensity(PSTNP)特征提取方法基于单链的位置特异性三核苷酸倾向描述RNA的统计意义，对于RNA序列，将生成4³＝64种三核苷酸，其三核苷酸位置特异性可以用64×(L-2)的矩阵表示，对于长度为L的RNA序列样本生成L-2维的特征向量。

更进一步的技术方案是所述步骤3)特征选择利用ElasticNet设置惩罚参数λ₁为0.1，惩罚参数λ₂为0.05，共剔除原始特征集合中346维冗余特征，保留165维对模型识别具有重要意义的特征，得到最优特征子集。

更进一步的技术方案是所述步骤4)预测模型是训练模型主要包括两个阶段的学习组成，分别为第一阶段的基分类器和第二阶段的元分类器；第一阶段根据初始数据集的特征训练出初级学习器，针对最佳特征集合，选择两个LightGBM和两个SVM作为基分类器，得到预测概率的输出值；第二阶段将初级分类器的概率输出作为输入特征，以纠正第一层的不准确训练，减少泛化误差，用于训练次级学习器，将第一阶段的概率输出值和最佳特征集合进行组合得到组合特征，选择SVM作为元分类器。

总体而言，本发明所提出的不同物种RNA中m⁶A位点预测方法与其他现有方法相比，具有以下技术优点：

本发明从序列信息、物理化学信息等方面对RNA序列进行特征编码，充分反映m⁶A位点与非m⁶A位点RNA序列之间的差异性。首次利用Elastic Net对原始特诊空间进行筛选和优化，得到最佳特征子集，提高模型计算速度和质量。将最佳特征子集输入到基分类器LightGBM和SVM中。元分类器集成第一阶段多个基分类器的概率输出值，学习不同的预测变量和真实类别之间的关系，以增强模型的预测性能。在第二阶段，我们将最佳特征集合和第一阶段的概率输出作为新的组合特征，分别输入到SVM和LR，根据预测结果选择SVM作为元分类器。

本发明通过在已知数据集上测试，并与其他方法比较，表明该发明在m⁶A位点识别方面具有较好的预测性能。本发明能够不仅能够帮助深入了解RNA甲基化修饰，而且能为生物学家进行跨物种m⁶A位点识别相关实验验证进一步提供有价值的参考信息。

附图说明

图1为本发明基于stacking集成的RNA中m⁶A位点预测方法的流程图。

图2为本发明不同预测方法关于训练集S.cerevisiae的ROC和PR曲线对比图，图中的A为ROC曲线，B为PR曲线。

具体实施方式

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

一种基于stacking集成的不同物种RNA中m⁶A位点预测方法，如图1所示，包括以下步骤：

1)收集RNA中N⁶-甲基腺苷修饰(m⁶A)位点信息：获得3个不同物种RNA的N⁶-甲基腺苷修饰位点数据集，包括正负数据集的RNA样本序列以及所对应的类别标签。

1-1)选择S.cerevisiae、H.sapiens和A.thaliana三个物种真实的数据；数据集S.cerevisiae作为训练集，数据集H.sapiens和A.thaliana作为测试集。其中S.cerevisiae数据集来源于W.Chen,P.Feng,H.Ding,H.Lin,K.C.Chou,iRNA-Methyl:Identifying N⁶-methyladenosine sites using pseudo nucleotide composition,Anal.Biochem.490(2015)26-33；H.sapiens数据集来源于P.Feng,H.Ding,H.Yang,W.Chen,H.Lin,K.C.Chou,iRNA-PseColl:identifying the occurrence sites ofdifferent RNA modificationsby incorporating collective effects of nucleotides into PseKNC,Mol.Ther.-Nucl.Acids,7(2017)155-163；A.thaliana数据集来源于W.Chen,P.Feng,H.Ding,H.Lin,IdentifyingN⁶-methyladenosine sites in the Arabidopsis thalianatranscriptome,Mol.Genet.Genomics 291(6)(2016)2225-2229。

1-2)训练集和测试集中均包括含有N⁶-甲基腺苷修饰位点的RNA(正样本)和非N⁶-甲基腺苷修饰位点的RNA(负样本)。

1-3)S.cerevisiae数据集样本序列是以共同基序GAC为中心，窗口大小为51。如果样本序列在某些位置不存在核苷酸，缺少的核苷酸将填充为它的镜像。数据集中包括1307条正样本，为了避免不平衡数据集对构建稳健模型的影响，从包含33,280条负样本的数据集中随机选取1307条负样本

1-4)H.sapiens数据集下样本窗口大小为41，若序列样本实际长度小于窗口大小，则缺少的核苷酸用其最临近的相同核苷酸填充，为了避免生成高度倾斜的数据集，该测试集中正负数据集数目保持一致，均为1130；A.thaliana数据集使用CD-HIT程序剔除序列相似性高于60％的样本，减少序列同源性偏差，数据集中包含394条正样本和394条负样本，样本窗口大小为25。

2)特征编码。使用序列衍生信息编码RNA序列，将字符信息转化成机器学习可识别的数值向量，确定PseDNC的最佳参数，融合6种特征提取方法对应的不同类型信息的特征向量，得到原始数据集的初始特征空间。

6种特征提取方法包括：binary encoding，chemical property，nucleotidefrequency，k-mer nucleotide frequency，pseudo dinucleotide composition(PseDNC)和position-specific trinucleotide propensity(PSTNP)

2-1)Binary encoding特征提取方法将能够准确描述样本序列中每个位置的核苷酸，该特征提取方法将RNA序列中包含的四种核苷酸：腺嘌呤(A)，鸟嘌呤(G)，胞嘧啶(C)和尿嘧啶(U)依次编码成4维二进制向量，分别为(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)。对于窗口大小为L的样本，将生成4×L维的特征向量。

2-2)Chemical property特征提取方法考虑到四种核苷酸具有不同的化学性质，依据环数、化学功能以及二级结构将A,C,G,U划分为三种不同类型的组。对环数而言，腺嘌呤和鸟嘌呤具有两个环，胞嘧啶和尿嘧啶只有一个环；对化学功能而言，腺嘌呤和胞嘧啶含有氨基团，而鸟嘌呤和尿嘧啶含有酮基团；对形成二级结构而言，鸟嘌呤和胞嘧啶具有强氢键，而腺嘌呤和尿嘧啶具有弱氢键；RNA序列中的每个核苷酸按照公式(1)编码成不同的向量s_i＝(x_i,y_i,z_i)。

根据以上三种划分方式，'A'可以用向量(1,1,1)表示，'C'可以用向量(0,1,0)表示，'G'可以用向量(1,0,0)表示，'U'可以用向量(0,0,1)表示。所以，对于窗口大小为L的样本，将生成3×L维的特征向量。

2-3)Nucleotide frequency特征提取方法计算RNA序列特中定位置核苷酸的密度，反映m⁶A位点附近核苷酸的组成和频率，定义如下：

其中N_i表示第i个子序列的长度，l表示子序列中每个核苷酸位置，q∈{A,C,G,U}。例如，RNA序列'GAUCACCG'，'A'在序列位置2和5的密度为1/2和2/5，'C'在序列位置4、6和7的密度分别为1/4，1/3，3/7，'G'在序列位置1和8的密度分别为1和1/4，'U'在序列位置3的密度为1/3。在计算样本序列每个位置的密度时，无论第一个位置以哪种核苷酸类型存在，其出现的密度均为1。因此对于长度为L的样本，将生成L-1维的特征向量。

2-4)K-mer nucleotide frequency特征提取方法基于相邻核苷酸对会影响RNA序列的结构和功能的原理，计算相邻核苷酸在样本序列中出现的频率，反映m⁶A位点和非m⁶A位点序列背景之间的差异性，运用K-mer算法计算相邻核苷酸在样本序列中出现的频率，将生成4^K维特征向量。随着K的增加，特征向量维数呈指数增加，会导致预测模型产生过拟合问题。因此选择2-mer特征提取方法编码样本序列，计算RNA序列中AA,AC,AG,AU,CA,CC,CG,CU,GA,GC,GG,GU,UA,UC,UG,UU出现的频率，生成16维的特征向量。

2-5)Pseudo dinucleotide composition(PseDNC)特征提取方法考虑到RNA序列的局部和全局序列信息，生成16+λ维的特征向量，其中λ值为23。定义如下：

R＝[d₁,d₂,…,d₁₅,d₁₆,,d₁₆₊₁,…,d_16+λ] (3)

其中：

公式(4)中，前16个元素反映了RNA序列中的短程或者局部序列顺序信息，而剩余的元素表示远程或者全局的序列顺序信息。λ是反映远程或者全局序列信息的伪分量总数，ω是权重因子，f_u表示RNA序列中第u个出现的非重叠二核苷酸的归一化出现频率，θ_j表示第j层序列相关因子。

其中θ₁称为第一层相关因子，反映所有沿着RNA序列最相邻的二核苷酸的顺序相关性，θ₂是第二最连续的二核苷酸之间的相关因子，依此类推。上式(5)的耦合因子C_i,i+j定义如下：

其中u表示RNA物理化学特性的数量，物理化学特性用于计算远程或者全局的序列顺序信息，在将其带入公式(6)时，需要对P_g(D_i)进行标准化，定义如下：

其中，符号<>表示取数量的平均值，SD表示相应的标准偏差。

2-6)Position-specific trinucleotide propensity(PSTNP)特征提取方法基于单链的位置特异性三核苷酸倾向描述RNA的统计意义，对于RNA序列，将生成4³＝64种三核苷酸，对于长度为L的RNA序列样本，其三核苷酸位置特异性可以用64×(L-2)的矩阵表示，定义如下：

其中z_i,j＝F⁺(3mer_i|j)-F^-(3mer_i|j),i＝1,2,…,64；j＝1,2,…,L-2。

F⁺(3mer_i|j)和F^-(3mer_i|j)分别表示正数据集(S⁺)和负数据集(S^-)中第i个核苷酸在第j个位置的出现频率。3mer₁等同于AAA，3mer₂等同于AAC，…，3mer₆₄等同于UUU。

对于给定的样本序列可表示为：

P＝[p₁,p₂,…,p_L-2]^T (9)

其中T是转置运算符，p_u定义如下：

3)特征选择。对初始特征空间，利用弹性网络Elastic Net对其进行降维，剔除冗余、噪声特征，保留对模型分类相关的重要特征，得到最佳特征集合，Elastic Net的最小化目标函数为:

其中，λ₁和λ₂是非负的惩罚参数，λ₂表示Elastic Net中l₁和l₂惩罚之间的缩放，取值范围在0和1之间。对于λ₂＝0，惩罚是l₂，即为Ridge。对于λ₂＝1，此时损失是关于l₁的惩罚，即为Lasso。当0＜λ₂＜1时，惩罚是l₁和l₂的组合。

Elastic Net惩罚参数λ₁设置为0.1，惩罚参数λ₂设置为0.05，共剔除原始特征集合中346维冗余特征，保留165维对模型识别具有重要意义的特征。

4)预测模型算法。把Elastic Net对应的最优特征子集以及所对应的类别标签输入到stacking集成中进行模型训练，并结合评价指标评估模型的预测性能。

模型训练主要包括两个阶段的学习组成，分别为第一阶段的基分类器和第二阶段的元分类器；第一阶段根据初始数据集的特征训练出初级学习器，针对最佳特征集合，选择两个LightGBM和两个SVM作为基分类器，得到预测概率的输出值；第二阶段将初级分类器的概率输出作为输入特征，以纠正第一层的不准确训练，减少泛化误差，用于训练次级学习器，本发明将第一阶段的概率输出值和最佳特征集合进行组合得到组合特征，选择SVM作为元分类器。

Stacking集成通过层级学习可以挖掘表征RNA中m⁶A位点的本质抽象特征，预测性能优于单分类器，伪代码见Algorithm 1。

5)m⁶A位点预测：将独立测试集中待预测的RNA序列，输入至步骤4)中的预测模型中，预测该RNA序列是否包含m⁶A位点并输出。

6)模型评估

本文选择最准确严谨的刀切法作为验证方法，数据集中的每个样本都被挑选出来作为独立的测试样本，而剩余的样本作为训练集训练模型。为了公正客观的评估预测模型的性能，选择四个常见指标：灵敏度(Sn)，特异度(Sp)，准确度(ACC)，Mathew's相关系数(MCC)，计算如下：

其中TP，TN，FP，FN分别表示真正例，真反例，假正例和假反例。TP表示正确预测真m⁶A位点的数量，TN表示正确预测非m⁶A位点的数量，FP表示非m⁶A位点被预测为真m⁶A位点的数量，FN表示预测非m⁶A位点为真正m⁶A位点的数量。ROC曲线和PR曲线也用来评估模型的鲁棒性及预测性能。ROC曲线是基于真阳性率(灵敏度)和假阳性率的曲线，PR曲线描绘精度(TP在所有预测的阳性结果中的比例)和召回率(灵敏度)的关系，该曲线比ROC曲线对误报更敏感。ROC曲线和PR曲线下的面积值分别记作AUC和AUPR，面积值越接近于1，模型预测性能越好。

6-1)6种特征提取方法融合后的特征集合预测准确性高：将单个特征Nucleotidefrequency(ANF)，Binary encoding(Binary)，Chemical property(NCP)，K-mernucleotide frequency(K-mer)，Pseudo dinucleotide composition(PseDNC)，Position-specific trinucleotide propensity(PSTNP)和融合后的特征集合All分别输入到基分类器LightGBM和SVM中，得到不同特征提取方法关于训练集S.cerevisiae中m⁶A位点的预测准确率如表1所示。

表1训练集S.cerevisiae关于不同特征提取方法的预测准确率

从表1可以看出，整体而言，分类器LightGBM和SVM关于S.cerevisiae的不同特征提取方法的预测准确度不同，即不同特征提取方法在对m⁶A位点识别的贡献程度是不一样的。对于单个特征提取方法PSTNP，基分类器LightGBM和SVM的预测准确率均达到最高，分别为79.42％和80.72％，分别比特征提取方法ANF高18.13％和19.17％，比特征提取方法K-mer高17.83％和17.60％。其次是特征提取方法Binary和Chemical的预测准确率较高，分类器LightGBM的预测准确率分别达到72.07％和72.72％，分类器SVM的预测准确率分别达到72.88％和73.22％。PseDNC关于数据集S.cerevisiae的基分类器预测准确率分别达到63.16％和64.38％。而基分类器LightGBM关于融合后的特征集合All的预测准确率均高于6种单独特征提取方法，表明多信息融合能够整合多种类型的信息，在一定程度上提高模型的预测准确率。而基分类器SVM关于All的预测准确率仅低于PSTNP的对应值，表明多信息融合会产生冗余特征，导致模型预测准确率的降低。

6-2)Elastic Net的降维方法能有效提高预测准确性：采用选择局部流行嵌入(LLE)，最大相关最大距离(MRMD)，谱聚类(SE)，奇异值分解(SVD)，互信息(MI)，ET等6种降维方法与Elastic Net进行对比。对S.cerevisiae数据集融合后的特征集合All进行特征优化，使用上述7种方法进行降维，将不同降维方法所对应的特征子集分别输入基分类器LightGBM和SVM，得到不同降维方法关于训练集S.cerevisiae预测准确率如表2所示。

表2训练集S.cerevisiae关于不同降维方法的预测准确率

从表2可以看出，对于训练集S.cerevisiae，基分类器LightGBM和SVM对于不同的降维方法的预测准确率有所不同，而且对于相同的降维方法基分类器对应的预测准确率之间的差异非常微小。其中，局部流行嵌入的降维效果最差，预测准确率分别达到69.01％和69.51％，甚至低于单个特征提取方法。相比于单个特征提取方法，降维方法MRMD，SE，SVD和MI能够保留对模型分类重要的特征，在一定程度上提高了模型的预测准确率，但是这四种降维方法对应的预测准确率仍低于原始特征集合All的对应值。相比于其他6种降维方法，Elastic Net具有最佳的降维效果，基分类器LightGBM和SVM对应的预测准确率分别达到80.07％和81.18％，比ET对应的预测准确率分别高0.69％和1.15％。Elastic Net能有效过滤对模型分类不相关、不重要的特征，保留对预测模型有贡献的特征，有效区分真实的m⁶A位点和非m⁶A位点，提高模型的预测性能。

6-3)为验证本发明中所述的stacking集成(下文简称为StackRAM)的预测准确性，对比AdaBoost，ERT，KNN，XGBoost，RF，LightGBM和SVM等7种分类器对训练集S.cerevisiae的m⁶A位点的预测性能。

AdaBoost通过不断调整训练集，设置学习率为0.1，将弱学习器联合得到强分类器。极端随机树分类器构建1000棵决策随机树，并且使用基尼指标来分裂节点。KNN通过学习距离其最近的50个训练样本的特征来完成位点识别任务。随机森林通过随机选择某些样本以及特征，避免模型过拟合并且获得很好的抗噪能力，本文构建具有1000棵的森林。XGBoost设置学习率为0.01，树的最大深度为10。LightGBM是基于树的梯度提升算法，设置最大深度为15识别RNA序列中的m⁶A位点。支持向量机通过径向基核函数将原始特征空间映射到高维空间，对N⁶-甲基腺苷位点进行识别。通过将Elastic Net对应的最佳特征集合分别输入到以上8种预测算法中，得到8种预测算法关于训练集S.cerevisiae的ROC和PR曲线对比图如图2所示。

根据图2可以直观的看出，StackRAM关于训练集S.cerevisiae的ROC和PR曲线均包含了其它7种分类器所对应的曲线，表明集成算法能够提高模型的鲁棒性。此外，对比ROC和PR曲线下所对应的面积AUC和AUPR值发现，StackRAM的AUC为0.9021，分别比AdaBoost，ERT，KNN，XGBoost，RF，LightGBM和SVM高8.26％，5.35％，4.25％，3.45％，3.41％，1.74％和1.63％。StackRAM的AUPR值为0.9022，分别比AdaBoost，ERT，KNN，XGBoost，RF，LightGBM和SVM高9.09％，5.68％，4.43％，3.85％，3.27％，2.26％和1.5％。相比于其它学习器，StackRAM通过将单个分类器拟合起来以获得具有高泛化性能的组合学习，学习不同的预测变量和真实类别之间的关系，有效挖掘表征RNA中m⁶A位点的序列特征。

6-4)为了对比预测模型，评估模型的鲁棒性以及预测性能，利用测试集进行验证。为了检验StackRAM的泛化性能，选择H.sapiens和A.thaliana作为独立测试集评估本文提出的新方法，得到StackRAM与其它预测方法关于独立测试集的结果对比如表3所示。

表3 StackRAM与其它方法关于测试集的结果对比

其中Feng's method来源于P.Feng,H.Ding,H.Yang,W.Chen,H.Lin,K.C.Chou,iRNA-PseColl:identifying the occurrence sites ofdifferent RNA modificationsby incorporating collective effects of nucleotides into PseKNC,Mol.Ther.-Nucl.Acids,7(2017)155-163；Chen's method来源于W.Chen,P.Feng,H.Ding,H.Lin,Identifying N⁶-methyladenosine sites in the Arabidopsis thalianatranscriptome,Mol.Genet.Genomics 291(6)(2016)2225-2229。

从表3可以看出，相比于其它方法关于独立测试集的预测结果，StackRAM具有识别m⁶A位点的优势。StackRAM关于H.sapiens预测准确率为92.30％，MCC值为0.8496，AUC值为0.9617，分别比Feng's method高1.92％，2.96％和11.27％。对于数据集A.thaliana，尽管Chen's method的Sp达到100％，但是ACC，Sn，MCC和AUC分别比StackRAM低2.67％，14.98％，2.27％和10.71％。总之，这些结果进一步验证了StackRAM的有效性和鲁棒性，表明StackRAM是一种功能强大的预测方法，不仅对训练集m⁶A位点识别具有竞争力，而且在跨物种位点识别方面具有更好的预测性能。

尽管这里参照本发明的多个解释性实施例对本发明进行了描述，但是，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的范围之内。对于本领域技术人员来说，其他的用途也将是明显的。

Claims (1)

1.一种基于stacking集成的RNA中N⁶-甲基腺苷修饰位点预测方法，其特征在于包括以下步骤：

1)收集RNA中N⁶-甲基腺苷修饰位点信息：获得3个不同物种RNA的N⁶-甲基腺苷修饰位点数据集，包括正负数据集的RNA样本序列以及所对应的类别标签，确定训练集与测试集；

2)特征编码：使用序列衍生信息编码RNA序列，将字符信息转化成机器学习可识别的数值向量，融合6种特征提取方法对应的不同类型信息的特征向量，得到原始数据集的初始特征空间；

3)特征选择：对初始特征空间，利用弹性网络对其进行降维，剔除冗余、噪声特征，保留对模型分类相关的重要特征，得到最佳特征集合；

4)训练预测模型：把弹性网络对应的最优特征子集以及所对应的类别标签输入到stacking集成中进行模型训练，并结合评价指标评估模型的预测性能，得到预测模型；

5)N⁶-甲基腺苷修饰位点预测：将测试集中待预测的RNA序列，输入至步骤4)中的预测模型中，预测该RNA序列是否包含N⁶-甲基腺苷修饰位点并输出；

所述步骤1)包括如下子步骤：

1-1)选择酿酒酵母、智人和拟南芥三个物种真实的数据；其中酿酒酵母的数据集S.cerevisiae作为训练集，智人的数据集H.sapiens和拟南芥的数据集A.thaliana作为测试集；

1-2)训练集和测试集中含有N⁶-甲基腺苷修饰位点的RNA作为正样本，非N⁶-甲基腺苷修饰位点的RNA作为负样本；

1-3)S.cerevisiae数据集样本序列是以共同基序GAC为中心，窗口大小为51，当样本序列在某些位置不存在核苷酸时，缺少核苷酸使用其对称位置的核苷酸填充；数据集中包括1307条正样本以及从总数为33,280条负样本的数据集中随机选取的1307条负样本；

1-4)H.sapiens数据集中样本窗口大小为41，若序列样本实际长度小于窗口大小，则缺少的核苷酸用其最临近的相同核苷酸填充；该测试集中正负数据集数目保持一致，均为1130；A.thaliana数据集使用CD-HIT程序剔除序列相似性高于60％的样本，减少序列同源性偏差，数据集中包含394条正样本和394条负样本，样本窗口大小为25；

所述步骤2)中，6种特征编码方法为二进制编码、化学性质、核苷酸频率、k-mer核苷酸频率、伪二核苷酸组成和位置特异性三核苷酸倾向特征提取方法；

所述二进制编码特征提取方法能够准确描述样本序列中每个位置的核苷酸，该特征提取方法将RNA序列中包含的四种核苷酸：腺嘌呤(A)，鸟嘌呤(G)，胞嘧啶(C)和尿嘧啶(U)，依次编码成4维二进制向量；

所述化学性质特征提取方法考虑到四种核苷酸具有不同的化学性质，依据环数、化学功能以及二级结构将A,C,G,U可以划分为三种不同类型的组，每种核苷酸用三维向量表示；

所述核苷酸频率特征提取方法计算RNA序列中特定位置核苷酸的密度，反映N⁶-甲基腺苷修饰位点附近核苷酸的组成和频率；对于长度为L的样本，将生成L-1维的特征向量；

所述k-mer核苷酸频率特征提取方法基于相邻核苷酸对会影响RNA序列的结构和功能的原理，计算相邻核苷酸在样本序列中出现的频率，反映N⁶-甲基腺苷修饰位点和非N⁶-甲基腺苷修饰位点序列背景之间的差异性，每条样本序列生成16维的特征向量；

所述伪二核苷酸组成特征提取方法考虑到RNA序列的局部和全局序列信息，生成16+λ维的特征向量，前16个元素反映了RNA序列中的短程或者局部序列顺序信息，而剩余的元素表示远程或者全局的序列顺序信息，λ为23，对于每条样本序列生成39维的特征向量；

所述位置特异性三核苷酸倾向特征提取方法基于单链的位置特异性三核苷酸倾向描述RNA的统计意义，对于RNA序列，将生成4³＝64种三核苷酸，其三核苷酸位置特异性可以用64×(L-2)的矩阵表示，对于长度为L的RNA序列样本生成L-2维的特征向量；

所述步骤3)中，弹性网络方法使用l₁和l₂范数正则化进行训练，其中惩罚参数λ₁设置为0.1，惩罚参数λ₂设置为0.05，共剔除原始特征集合中346维冗余特征，保留165维对模型识别具有重要意义的特征，得到最优特征子集；

所述步骤4)中，训练模型主要包括两个阶段的学习组成，分别为第一阶段的基分类器和第二阶段的元分类器；第一阶段根据初始数据集的特征训练出初级学习器，针对最佳特征集合，选择两个LightGBM和两个SVM作为基分类器，得到预测概率的输出值；第二阶段将初级分类器的概率输出作为输入特征，以纠正第一层的不准确训练，减少泛化误差，用于训练次级学习器，将第一阶段的概率输出值和最佳特征集合进行组合得到组合特征，选择SVM作为元分类器。

Images