基于改进的MRMD算法和DF模型的植物pri-miRNA编码肽预测方法

基于改进的mrmd算法和df模型的植物pri-mirna编码肽预测方法
技术领域
1.本发明属于生物信息学和计算机科学技术领域，涉及一种改进的mrmd特征选择算法在生物序列数据上应用以及结合深度森林(deepforest，df)模型构建植物pri-mirna编码肽预测模型的方法。


背景技术：

2.生物体内存在数目巨大、种类繁多的rna序列，研究人员根据其是否能够编码蛋白将其分为两大类，编码rna和非编码rna(non-coding rna，ncrna)。随着研究的不断深入，发现部分ncrna也具备编码蛋白质的能力。微小rna(microrna，mirna)作为一类ncrna，在生命体的病毒防御、器官生成等生命体关键活动中发挥了重要作用。成熟的mirna由较长的初级转录物经过一系列核酸酶的剪切加工而产生，初级转录物称为mirna初级体(primary mirna，pri-mirna)。由于pri-mirna在被转录后，存在的时间较短，易被剪切形成mirna前体(precursor mirna，pre-mirna)，之前对它的研究并不多。然而近年来的研究表明，pri-mirna包含小开放阅读框(small open reading frames，sorf)，其编码的肽(mirna-encoded peptide，mipep)能够促进对应mirna积累，进而间接调节靶基因的表达，对于生物体具有重要的调节作用。mipep的发现为现代农业等领域的发展提供了新思路，使得对其的研究成为了热点。如何正确识别mipep对于深入研究mipep的功能有着至关重要的作用。
3.目前，对植物mipep采用计算方法进行预测属于较为新颖的研究方向，取得的成果相对较少。对它的研究主要有生物实验、生物信息学等方法。生物实验结果可信度高，然而费时费力，不能够进行大规模验证；生物信息学方法可以通过构建机器学习模型快速识别mipep，为生物学实验鉴定提供参考，从而节省大量的时间与花费。目前，针对小开放阅读框编码肽(sorf-encoded peptide，sep)的研究较少，已有的识别工具难以鉴定出高可信度的mipep，所以迫切需要一个快速、准确的生物信息学工具来辅助识别mipep。然而使用生物信息学对植物mipep进行预测面临很大挑战。首先，由于没有pri-mirna数据库，需要采用一定的方法对pre-mirna序列进行扩充从而得到数据集。同时，如何采用合适的特征提取与特征选择方法，以及选取合适的模型对pri-mirna中的sorf进行预测，这都是采用生物信息学方法挖掘植物mipep必须要解决的问题。


技术实现要素：

4.基于以上所述需要解决的问题，本发明针对mrmd特征选择算法对于特定生物序列数据处理效果欠佳的情况，采用更适用的评价系数，提出一种特征选择算法immrmd(improved max-relevance-max-distance based dimensionality reduction)，提高特征数据的代表性，同时结合df模型，实现对植物pri-mirna上的sorf的预测。本发明能够判断输入的rna序列是否符合pri-mirna上首个sorf的相关特征，从而为判断该序列是否编码肽提供依据，为生物实验提供参考，节省大量的人力物力。
5.本发明的技术方案为：
6.一种基于改进的mrmd算法和df模型的植物pri-mirna编码肽预测方法，包括以下步骤：
7.步骤1：构建实验数据集。其中，正集采用从拟南芥pri-mirna序列上得到的首个sorf，负集采用拟南芥非编码序列(non-coding sequence，ncds)数据。
8.具体过程如下：
9.1.1：首先，从mirbase数据库(https://www.mirbase.org/)中下载拟南芥的所有mirna数据，得到其前体mirna(precursor mirna，pre-mirna)。然后，通过ensemble plants(http://plants.ensembl.org/index.html)工具进行数据扩展，将得到的所有pre-mirna向前扩充1000个碱基，然后利用softberry软件(http://www.softberry.com/berry.phtml？topic＝tssplant&group＝programs&subgroup＝promoter)查找其转录起始位点(transcription start site，tss)，从而得到相对应的pri-mirna。再利用orf finder(https://www.ncbi.nlm.nih.gov/orffinder/)工具获取每条pri-mirna的第一个sorf序列，并通过cd-hit工具去除相似度高于80％的序列，得到候选的正集sorf数据。
10.1.2：收集拟南芥ncds数据，将其作为负集数据源。使用orf finder工具获取其中的sorfs序列，得到的结果同样去除相似度高于80％的sorf序列，得到候选的负集数据。
11.1.3：选取总数据量的10％，即从正集随机取出24条，负集随机取出25条，作为测试集，测试最终模型预测性能。
12.步骤2：提取sorf序列数据的相关特征。
13.根据正、负集数据在序列组成和理化特性方面存在的差异，从正、负集数据中均提取k-mer特征、短序列模体(short sequence motifs，ssm)特征、信噪比、gc碱基含量(gc_con)、gc碱基比例(gc_ra)和序列长度等特征。
14.其中，k-mer特征描述sorf序列的组成信息，它充分考虑了相邻核苷酸之间的依赖关系，保存了相邻核苷酸的信息，从而避免了信息的损失。k代表序列中相邻的碱基数，k个相邻碱基会有4k种组合。本发明中k为1,2,3,4，共340维。采用步长为1的滑动窗口进行滑动匹配，公式如下：
[0015][0016]
lk＝l-k 1,k＝1,2,3,4
ꢀꢀꢀ
(2)
[0017]
其中，lk为滑动次数，si为每个k-mer出现的次数，fi为最终每个k-mer的频率，l为序列长度，最终得到了340维特征。
[0018]
k-mer特征考虑了连续碱基的性质，然而不连续的碱基之间也存在差异，因此采用ssm特征进行补充。不同的碱基对可以记为a*b，a**b，a***b等，其中a和b均是a、t、c、g中的任意一种碱基，*代表任意碱基，每一组ssm都有16维特征表达，提取a*b、a**b、a***b三组ssm，共48维。
[0019][0020]
uj＝l-j-1,j＝1,2,3
ꢀꢀꢀ
(4)
[0021]
其中，uj为滑动次数，vi′
为每个短序列模体出现的次数，ssmi′
为每个短序列模体的出现频率。
[0022]
根据序列的本身表达，提取序列长度l、gc碱基含量gc_con、gc碱基比例gc_ra以及信噪比特征，共4维。信噪比的大小能够表示序列中碱基使用的偏向性，通过碱基在三个相位的分布计算sorf在三分之一处的功率谱来获取信噪比。具体公式如下：
[0023][0024][0025][0026][0027]
其中，g和c分别代表序列中碱基g和碱基c的数目；l代表序列长度；a
x
、b
x
和c
x
均为一行三列的数组，分别代表碱基x在三个相位出现的频数；w表示功率谱；snr代表信噪比。
[0028]
最终，提取的上述特征共组成392维特征数据。
[0029]
步骤3：设计immrmd算法，对步骤2所得到的特征数据采用immrmd算法进行特征排序和选择。
[0030]
由于本发明自身数据的特异性，对mrmd算法进行了相应改进。mrmd算法采用皮尔森(pearson)系数作为特征重要性的判断标准，然而其要求数据呈正态分布，并且无异常点，对数据要求相对较高，本发明的实验数据并不符合此情况，因此选用斯皮尔曼系数(spearman)作为判断标准，spearman系数要求的数据量少，符合本实验的数据情况，同时其评价时对数据是否呈现正态分布没有提出要求，更符合现实数据情况。spearman相关系数计算公式如公式(9)所示，定义x和y为两组数据，
[0031][0032]
其中，dy为xy和yy之间的等级差，rs位于-1和1之间；n为样本总数。
[0033]
另一方面选择调整余弦相似度代替余弦相似度作为最大距离的判断标准。余弦相似度是内积空间中两个非零向量之间的相似度的度量，度量它们之间夹角的余弦值，注重的是从方向上区分的差异，而对绝对的数值不敏感，因此无法衡量各维数值的差异。调整余弦相似度是基于向量相似度的一种改进形式，在考虑方向上的差异的同时，也考虑了数值上的差异性，其目的是解决余弦相似度仅考虑向量维度方向上的相似而忽略各个维度的量纲的差异性的问题，所以在计算相似度的时候，做了每个维度减去均值的修正操作，弥补了余弦相似度的不足。
[0034]
综上所述，本发明在原有mrmd特征选择算法的基础上，对其原有评价标准进行了相应改进，设计了immrmd算法；immrmd算法的具体步骤如下：将步骤2提取出的原始特征数据作为immrmd算法的输入，首先利用spearman相关系数来计算子特征集与目标类别之间的
相关性，选择与目标类别相关性最高的子特征集；其次，考虑到减少特征之间的冗余性能够更好的提取相关特征，采用最大距离来衡量两个特征向量之间的相似度。计算euclidean距离、调整余弦相似度和tanimoto系数，得出任意一个特征向量与其余特征向量的冗余性程度。最后，通过对两部分的计算结果进行累加并排序，将排序后的特征数据逐维输入到逻辑回归模型中进行评价，选择其中分类效果最好的前m维，从而得到同时具备最大相关性和最小冗余性的m维子特征集。
[0035]
最终，原mrmd3.0算法将392维特征数据降维到314维，而本发明的immrmd算法在改变评价标准后，筛选出382维。
[0036]
步骤4：利用df模型得到预测结果。
[0037]
将步骤3中筛选出的数据输入到df模型中，得到最终的预测结果。
[0038]
本发明的有益效果：
[0039]
(1)综合考虑了基因序列中核苷酸的依赖联系，提取了k-mer、ssm等相关特征。
[0040]
(2)对mrmd特征选择算法进行相应改进，提出immrmd特征选择算法，使得选择后的特征数据更具有代表性，提高了分类效果。
[0041]
(3)在多种模型比较下选择df作为分类模型，获取关键特征信息，从而进一步提高了本发明对于植物pri-mirna编码肽的预测性能。
附图说明
[0042]
图1为本发明的整体框架流程图。
[0043]
图2为本发明原始特征数据在不同模型的分类效果对比图。
[0044]
图3为本发明在经过特征选择算法后特征数据在svm上的分类效果图。
[0045]
图4为本发明在经过特征选择算法后特征数据在lr上的分类效果图。
具体实施方式
[0046]
以下结合附图附表和技术方案，进一步说明本发明的具体实施方式。
[0047]
本发明提供了一种结合特征选择算法与机器学习的pri-mirna编码肽预测方法，通过对生物信息学中常用的mrmd算法进行改进提取更具有代表性的特征，然后选择预测效果更好的df作为预测模型，从而实现对于pri-mirna编码肽的预测，该方法整体框架如图1所示，步骤如下：
[0048]
步骤1：获取拟南芥相应的pri-mirna中的sorf序列作为正集以及ncds序列作为负集。
[0049]
1.1：通过mirbase数据库下载拟南芥的所有mirna数据，得到其pre-mirna数据；然后，通过ensemble plants工具进行数据扩展，将所有的pre-mirna向前扩充1000个碱基，并通过softberry软件对其进行数据处理，从而得到所需要的pri-mirna的首个sorf序列，通过cd-hit工具去除相似度高于80％的序列，得到候选的正集sorf数据，共245条。
[0050]
1.2：下载拟南芥的ncds数据，使用orf finder工具获取其中的sorfs序列，同样通过cd-hit工具去除相似度高于80％的序列，得到候选负集数据，共1670条。因为正负集数据差距过大，所以从1670条中随机抽取245条作为候选负集数据。
[0051]
1.3：在正负集数据中随机挑选10％，即正集挑选24条，负集挑选25条，作为模型测
试集。
[0052]
步骤2：提取sorf序列数据的相关特征。
[0053]
研究表明，提取的特征越丰富，模型的性能就越好。首先，为充分考虑相邻核苷酸之间的依赖联系，引入了k-mer特征，该特征描述sorf序列的组成信息，它保存相邻核苷酸之间的联系，避免了信息的损失。k代表序列中相邻的碱基数目，k个相邻碱基就会有4k种组合方式，本发明分别采用k为1,2,3,4，共得到特征340维。采用步长为1的滑动窗口进行滑动匹配。
[0054]
k-mer特征考虑了连续碱基的性质，然而不连续的碱基之间也存在差异，为了防止丢失这一部分信息，采用ssm特征进行补充。不同的碱基对组成可以获得不同信息，可以记为a*b，a**b，a***b等，其中a、b分别记为a、t、c、g中的任意一种碱基，*代表任意碱基，每一组不同的组成都由16维特征表达，提取了a*b，a**b，a***b三组的ssm特征，总共48维。
[0055]
除了考虑碱基对的不同信息之外，还根据序列的本身表达，提取了序列长度l、gc碱基含量、gc碱基比例以及信噪比等特征。以上组成4维特征。
[0056]
将以上特征进行组合，共组成392维特征集。
[0057]
步骤3：设计immrmd算法，进行特征排序与选择。
[0058]
mrmd算法常被用于生物数据降维方面。mrmd特征选择方法主要由两部分决定：其一是特征和实例类标之间的相关性，mrmd用pearson相关系数来计算特征和类标之间的相关性；其二是特征之间的冗余性，用三种距离函数(euclidean距离，cosine距离和tanimoto系数)计算特征之间的冗余性。pearson相关系数越大说明特征与类标关系越紧密，距离越大说明特征之间的冗余性越低。最后，mrmd算法选出来的是和类标具有强相关并且特征之间具有低冗余性的特征子集。
[0059]
然而pearson相关系数所处理的数据要求必须符合正态分布并且无异常点，与本实验数据不符，因此本发明选取更常用于核苷酸序列研究的spearman相关系数进行计算，spearman系数对于数据量的要求较低，并且运行数据不要求呈现出正态分布，更符合现实数据情况。同时，计算距离时的cosine距离仅考虑了数据在方向上的差异，尽管其对个体间存在的偏见可以进行一定的修正，但是因为只能分辨个体在维度方向之间的差异，没法衡量每个维数值的差异，因此会因为对于绝对的数值差异敏感度较差从而导致结果的偏差，因此本发明为了修正这种仅考虑向量维度方向上的相似而不考虑各个维度量纲的差异的不合理性，采用调整余弦相似性进行计算，即在所有维度上的数值都做一个减去均值的修正操作，使其计算结果更加符合现实。
[0060]
具体地，将步骤2中所提取出的原始特征数据作为immrmd特征选择算法的输入，一方面计算spearman相关系数获取子特征与分类类别之间的相关性，每一维特征数据与其余特征的相关性都由(-1,1)之间大小的数字表示；另一方面通过计算euclidean距离、调整余弦相似度和tanimoto系数得到每一个特征向量与其余特征向量的冗余性程度，同样采用数字表示；并将所得两方面的计算结果进行累加并排序，逐维输入到逻辑回归模型中对其进行分类效果的评价，采用f1值等评价标准，选择出效果最好的前m维特征，从而得到同时具备最大相关性和最小冗余性的m维子特征集。
[0061]
步骤4：采用svm、lr、df等模型进行对比实验，选择效果最好的df模型作为预测模型。
[0062]
df是随机森林的延展性算法，其综合了深度学习思想，是一种决策树集成方法。df具有低超参数、高效率，低数据量也可获得良好学习性能，适用于多种领域数据等优势。
[0063]
df由以下两个部分组成：级联森林和多粒度扫描。多粒度扫描是受卷积神经网络启发，产生级联森林的输入特征向量。滑动窗口用于扫描原始特征，采用不同采样窗口可以增强对输入数据的特征表达能力。例如，对于400维的输入特征数据，如果采用大小为100维的滑动窗口对输入特征进行处理，在步长为1的情况下，最终就会得到301个特征向量。经过窗口扫描生成的类向量将输入随机森林和完全随机森林生成变换特征向量。假设是一个二分类问题，则每个森林将会产生301个2维类向量，最后对两个森林产生的类向量进行拼接，输出对应于原始400维特征向量的1204维变换特征向量。为了提高模型特征的多样性，多粒度扫描通常使用多个不同大小的窗口对样本数据进行多尺度采样，从而生成尺度不一的特征向量，获得更多的样本信息。df模型规定了默认的多粒度扫描结构参数，对于具有d维特征的原始数据，将使用大小分别为d/16、d/8和d/4的特征窗口进行扫描，在保证生成的样本的多样性的情况下，使计算效率最大化。
[0064]
级联森林结构的每层都通过级联的方式堆叠若干个随机森林，使df具有了表征能，增强了模型的特征表示能力以及泛化能力。每个随机森林都包含若干个决策树，每个决策树都输出一个以类向量形式的结果，然后对每个决策树的输出类向量结果取均值。第一层级联森林以多粒度扫描的结果作为输入，此后每层级联的输入都是由经过级联后的增强特征以及多粒度扫描的结果拼接而成。考虑到深度结构通常出现过拟合问题，在每一层结束后都会在测试集上对结果进行验证，如果效果有所提升，则将增强向量继续传递给下一层级联，从而产生新的拼接向量，重复以上过程；如果没有提升，则终止训练。因此，级联森林的层数能够自动确定，df也正是通过该结构来自动确定深度结构的复杂度从而实现对不同大小的数据集的适应性，也正是这个特性，使得它与其他深度神经网络相比，不仅适用于大规模的训练数据，在小规模的训练数据集上也能够获得不错的结果。
[0065]
将实验数据集在svm、lr以及df模型上分别进行实验，本发明以准确率(accuracy，acc)、精确率(precision，p)、召回率(recall，r)和f1值(f1_score，f1)作为评价标准，四种评价标准的计算公式如下：
[0066][0067][0068][0069][0070]
其中，tp，fp，tn，fn的含义如表1所示。
[0071]
表1分类含义表
[0072][0073]
结果如附图2及附表2所示，在df上对于测试集的分类效果最好。在经过特征选择算法mrmd3.0以及immrmd处理后的特征数据集，通过df之后的分类效果在对于正负集的预测能力均衡，且acc、p值和f1值上均有所提升，具有更好的预测效果。
[0074]
表2基于拟南芥数据集采用不同的特征处理方式后在df下的分类性能表
[0075][0076]
同时，为了更好的说明immrmd在不同情况下对于本实验分类效果的提升作用，本发明还分别做了以下实验：图2表明利用原始特征数据，通过df、svm和lr分别进行分类预测时，df便表现出相比较更好的分类效果，在acc、p值、r值和f1值上都比其他模型表现优异。图3、图4则分别表示在经过mrmd3.0、immrmd和未经过特征选择算法的情况下，在svm和lr上的分类效果，由图3、图4可知，immrmd能够使得经特征选择后的数据在正负集标签预测上更加均衡，说明其更加适用于该序列数据，并且在acc和f1值上均有所提升。而表2则是采用df作为分类模型，通过不同方法对原始特征数据进行处理后的分类比较，表2中mrmd3.0-sc表示只改变mrmd3.0中的pearson系数为spearman系数，而不改变余弦相似度，autoencoder则表示采用自动编码器处理数据；由表2可知，immrmd在acc、p值和f1值上都有显著优势，且对于正负集的分类效果相比较其他更加均衡，说明其相比较其他方法分类效果更好。
[0077]
本发明的核心主要有两个方面，一是对mrmd特征选择算法进行改进，优化其评价标准，使用spearman系数代替pearson系数，使之更适用于本发明本身数据，同时以综合考虑数值与方向影响的调整余弦相似度代替之前的余弦相速度，对数据考虑更加全面；另一方面是结合深度学习思想，比较后采用更适合的df作为分类模型，提高分类效果。