CN108614955A - 一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法 - Google Patents

Abstract

本发明是一种新型的lncRNA鉴定方法。本发明提出了基于对数度量，多尺度二级结构，电子‑离子互作赝势的三种新型特征设计提取方式，并进而利用机器学习算法构建分类器。对数度量特征可大幅降低基于序列组成特征的维数，在保证高准确度的同时提升模型的效率；多尺度二级结构特征则可在结构层面挖掘更加保守的特征；电子‑离子互作赝势通过序列理化特征进一步提升本发明的跨物种稳定性。经实验验证，本发明针对人类数据集准确度高达97.28％，针对小鼠数据集准确度达93.47％，且本发明仅需56.01秒即可完成5000条序列的预测。相比其他算法，本发明有着更良好的准确度与效率，也具有更优秀的容错性与跨物种稳定性。

Description

一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法

技术领域：

本发明属于生物信息领域，尤其涉及lncRNA的鉴定。

背景技术：

在人类基因组中，有80％的序列具有生物学功能，而编码蛋白质的序列只占基因组的不到 2％。我们将长度大于200碱基，不能编码蛋白质的RNA称为长非编码RNA，简称lncRNA。 lncRNA与众多生物通路以及多种疾病，特别是恶性肿瘤，有着密切的联系。更有研究认为，lncRNA是解决肿瘤问题的最终落脚点。但lncRNA保守性较差，表达量也较低，与编码RNA十分相似，甚至与编码RNA存在序列上的重合区域。这些因素为lncRNA研究的第一步，即lncRNA的鉴定，带来了诸多困难。

随机器学习技术的快速发展，目前已有几种基于机器学习算法的lncRNA鉴定方法被提出。其主要原理是首先针对序列进行特征提取，而后利用机器学习算法构建lncRNA分类器。一个算法表现的优劣与特征提取有着非常直接的关系，关键的特征可以更加准确地对lncRNA进行鉴定。在lncRNA鉴定问题中，通常将lncRNA标记为正类(positive class)，编码RNA标记为负类(negativeclass)；并使用敏感度(Sensitivity)，特异度(Specificity)，准确度(Accuracy)与F-度量(F-measure)这几种指标来对算法的表现进行评估：

公式中P为正样本，N为负样本，TP是预测为正实际也为正的样本，TN为预测为负实际也为负的样本，而FP为预测为正而实际为负的样本，FN为预测为负而实际为正的样本。准确度可用来衡量正确预测所占的比例；敏感度与特异度可分别用来检测算法判断 lncRNA与编码RNA的能力；而F-度量则可综合评价算法对两类序列的判断能力。

目前主要的基于机器学习方法的lncRNA鉴定方法主要包括CPC，CNCI，PLEK与CPC2。这四种方法均为当前学界使用最为广泛的或代表最新科研进展的几种lncRNA鉴定手段。CPC由北京大学生物技术中心高歌课题组于2007年开发完成(Kong L,Zhang Y,Ye Z Q,et al.CPC:assess the protein-coding potential of transcripts using sequencefeatures and support vector machine[J].Nucleic Acids Research,2007,35(WebServer issue):W345.)。CPC 是基于序列比对的lncRNA鉴定方法的代表，其特征主要从开放阅读框信息与蛋白质序列比对信息两方面提取而得。CPC首先将待检测RNA序列翻译为蛋白质序列，而后将其与数据库中的蛋白质序列进行比对来提取比对信息特征。CPC认为由编码RNA翻译而来的蛋白质序列与数据库中的蛋白质序列之间，往往有数量更多，质量更优的匹配片段。但作为基于序列比对的lncRNA鉴定方法，CPC有诸多难以避免的缺陷：首先，大量lncRNA 与编码RNA非常相似，因此lncRNA翻译而得的蛋白质序列与数据库中的蛋白质序列之间，同样易于出现匹配片段，故CPC非常容易将lncRNA判断为编码RNA，造成敏感度较低。其次，CPC严重依赖于序列比对，对待检测序列与比对数据库的质量都有着较高的要求。然而目前测序技术得到的序列时常因信号较弱而出现测序误差，同时数据库中大量物种的注释信息又十分有限，很难为CPC提供充足的数据进行比对。因此CPC在对这类序列进行预测时，其结果不可避免地会产生较大误差，甚至因程序错误而无法进行预测。最后，序列的比对过程非常耗时，CPC可能需要数十小时才能完成几千条序列的预测，因此难以将CPC应用于大规模数据计算任务。值得关注的是，目前lncRNA的研究已经越来越聚焦于冷门物种的研究，而物种的序列往往是通过高通量测序技术得来，序列数量巨大，且碱基误差难以避免。因此CPC已越来越难以满足当前学术界对lncRNA研究领域提出的更高要求。

其他三种鉴定lncRNA的方法CNCI，PLEK与CPC2均不需要进行序列比对，一定程度上改善了CPC方法中因序列比对带来的不足，特别是在计算效率上有了较大的提升。其中CNCI由中国科学院计算技术研究所赵屹团队于2013年开发完成(Liang S,Luo H,Bu D, etal.Utilizing sequence intrinsic composition to classify protein-coding andlong non-coding transcripts[J].Nucleic Acids Research,2013,41(17):e166.)。CNCI主要利用的特征包括：序列最近似编码区片段上的相邻密码子碱基信息，以及3联碱基频率。与CPC相比，CNCI有着更加良好的敏感度，即在敏感度与特异度上取得了更好的平衡。但CNCI主要基于序列碱基频率进行预测，因不同物种序列的碱基频率千差万别，CNCI在不同物种上的表现也会有所波动；此外，由于需要寻找序列的最近似编码区片段，CNCI依然对序列的质量有较高要求。与CPC一样，CNCI无法判断某些存在碱基错误的序列。PLEK由西安电子科技大学计算机学院张军英团队于2014年开发完成(Li A,Zhang J,Zhou Z.PLEK:atool for predicting long non-coding RNAs and messenger RNAs based on animproved k-mer scheme[J]. Bmc Bioinformatics,2014,15(1):311.)。PLEK利用1至5联碱基频率进行lncRNA预测。与 CPC，CNCI相比，PLEK的特征设计相对简单，计算效率也得到了进一步提升。同时，PLEK 不需要像CPC一样进行数据比对，亦不需要像CNCI寻找序列的最近似编码片段，因此 PLEK对序列有良好的容错性，序列中存在的碱基错误不会明显地干扰PLEK的表现。但 PLEK较CNCI更加严重地依赖于序列碱基频率，因此PLEK在不同物种下的表现更加不稳定。为弥补方法CPC的不足，高歌团队于2017年提出了新的lncRNA鉴定算法，即CPC2(Kang Y J,Yang D C,Kong L,et al.CPC2:a fast and accurate coding potentialcalculator based on sequence intrinsic features.[J].Nucleic Acids Research,2017,45(W1).)。这一方法亦不再依赖于序列比对，其主要利用开放阅读框、等电点信息与密码子中碱基含量的偏好性来对lncRNA进行鉴定。CPC2的效率相比CPC，CNCI与PLEK有了大幅的提升，同时在敏感度与特异度上有了更加良好的平衡。但CPC2仍然未能取得优秀的跨物种稳定性。

综上，机器学习技术已逐渐成为学术界针对lncRNA鉴定进行研究的主流手段。但基于序列比对方法的CPC，无论表准确度还是效率，其表现均不优秀；而诸多非基于序列比对的方法，则几乎都单一地从序列组成这一角度对lncRNA进行考虑，仅CPC2从等电点这一理化性质特征来判断lncRNA与编码RNA之间的差异。不同物种间的序列组成往往存在差异，仅利用序列组成对lncRNA进行鉴定，其准确度往往会产生较大的波动；同时，仅从单一角度设计而得的特征，对lncRNA的识别能力相对有限，非常易于陷入瓶颈，难以进一步提升lncRNA鉴定的准确度。此外，因部分方法较差的容错性，使得其对高通量测序而得的序列无法进行计算。当前，面临着高通量测序技术带来的海量序列，研究人员亟需一种高准确度，高效率，容错度强，表现稳定的新型lncRNA鉴定方法。

发明内容：

本发明旨在解决lncRNA鉴定问题，提出一种能够广泛应用于高通量测序序列与冷门物种序列的，高准确度，高效率，高容错度且具有优秀跨物种稳定性的lncRNA鉴定方法。

有益效果：通过本发明，可以得到高准确度，高效率，同时稳定性与容错性优异的lncRNA分类器。本发明综合表现优秀，操作便捷，其关键技术在于三类异源特征的设计与提取。依据以上步骤说明，可以使用R或Python等语言高效便利地实现本发明，并进一步用于lncRNA的鉴定。具体优势如下：

1.创新性提出了基于对数度量方案的序列组成特征设计方法。这一设计方法包括特征的提取手段，以及特征提取过程中，k的取值与滑动步长均经过仔细设计与实验验证。现有的序列组成特征通常针对序列全长使用k-merscheme策略进行特征提取，k值取值为1至5，每次滑动1个碱基，且每一种碱基频率组合都作为特征使用。这一特征设计方式，不仅跨物种表现稳定性不佳，同时特征数目巨大，例如当k为6时，特征数目高达4096 条。使用本发明中的对数度量方案，特征将在最长开放阅读框中进行计算，k取值为6，且每次滑动3个碱基；当开放阅读框不存在时，k取值为6，每次滑动1个碱基。对数度量特征首先融合了序列组成与生物学意义，因此有着更高的准确度；其次利用对数度量方案提取特征，特征数目均为3条而不随于k值变化而改变。整体而言，本发明的对数度量方案有着更优秀的准确度与跨物种表现，大幅度提高了lncRNA的鉴定效率与稳定性。

2.创新性提出了基于多尺度二级结构的特征提取方式。这一提取方式包括多尺度二级结构序列的构造方式以及特征提取方法。目前鲜有使用二级结构为特征进行lncRNA鉴定的方法。而综合多尺度二级结构序列与对数度量特征提取方案进行lncRNA鉴定亦是首次。

3.创新性提出了基于理化性质的特征提取方式。包括以电子-离子互作赝势为基础的特征以及特征提取过程中的区间设计方式。目前仅有少数方法将RNA序列转换为0-1编码的数值序列，而后通过快速离散傅里叶变换得到的能量谱进行外显子定位或lncRNA分析，该方法需要对4条数值序列进行分析。而本发明设计的利用电子-离子互作赝势特征进行lncRNA鉴定尚属首次提出。电子-离子互作赝势特征不仅融入序列的理化性质，更仅需针对1条数值序列进行分析即可。且以能量谱不同区间的分位数统计信息作为特征，进行lncRNA分析，亦是首次提出。本发明中的这一类特征设计方式以编码序列的理化性质与周期3特性为基础，与现有的以数值信号为基础的傅里叶变换分析进行外显子定位的方法有本质区别。

4.本发明融合了以上三类特征对lncRNA进行鉴定。从方法设计方面而言，本发明的特征提取步骤简单快速；从鉴定表现而言，本发明有着优秀的准确度与鲁棒性，并同时具有出色的序列容错性与跨物种稳定性。

附图说明

图1为本发明流程图；

图2为多尺度二级结构序列的构造方法；

图3为编码RNA的电子-离子互作赝势能量谱示意图；

图4为lncRNA的电子-离子互作赝势能量谱示意图；

图5为不同区间长度上基于理化性质的分位数统计量特征的表现。

具体实施方式

以下结合附图对本发明的流程进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

本发明从人类基因入手，构建用于预测人类lncRNA的分类器，并基于此分类器，进一步构建用于预测其他物种lncRNA的分类器。以下将结合附图详细说明此其流程。

本发明用以对lncRNA进行鉴定。主要包括如下几个步骤：

第一步：构建训练集。

如图1示，从注释信息全面的数据库中收集序列，筛选高质量序列用以构建训练集。训练集的要求为正负类样本数量尽可能接近，同时样本来源可靠。因此可以通过以下方法构建高质量的训练集：

步骤1.1：选择GENCODE数据库，下载人类与小鼠的lncRNA序列与编码RNA序列。选择Ensembl数据库，下载诸如原鸡，斑马鱼等物种的lncRNA序列与编码RNA序列。同时，下载Ensembl数据库中相应物种的编码区序列。GENCODE为当前使用最为广泛，收录信息最为可靠的人类、小鼠数据库之一。而Ensembl数据库则包括大量物种的不同类别的基因序列，这一数据库也包括大量物种的编码区序列。编码区序列指编码RNA序列去除非编码区部分而得到的序列区域。

步骤1.2：对序列进行质量筛选。可筛去含有非正常碱基的序列(正常碱基包括a，c， g，t)，重复序列等；也可每条基因仅保留一条序列，或筛去相似度大于80％的序列。

步骤1.3：选择相同数目的数量足够的lncRNA，编码RNA或编码区序列，用以构建训练集。在本实例中，12190条人类lncRNA与相同数目的人类编码RNA序列，被用于构建训练集。且任意两条序列均来自于不同基因，且训练集序列均不含非正常碱基。

第二步：设计提取基因序列组成特征，这一特征大类主要将基因序列组成与生物学性质相结合，进而计算对数度量特征。具体包括以下步骤：

步骤2.1：对训练集中的序列查找最长开放阅读框。开放阅读框将在三种序列读码方式上分别进行计算，并选择当中最长的开放阅读框作为序列最终的最长开放阅读框。开放阅读框的判断方式可以是严格的以起始密码子开始，终止密码子结束的部分；也可以更加广泛地，当序列无终止密码子出现时，以起始密码子开始至此读码方式的最后一个密码子结束的部分作为开放阅读框。在得到最长开放阅读框后，计算最长开放阅读框的长度与覆盖率。覆盖率为最长开放阅读框长度与序列长度的比值。当开放阅读框不存在时，长度与覆盖率均为0。

步骤2.2：利用本发明所定义的对数度量方案，针对最长开放阅读框提取LogDist.LNC， LogDist.PCT与LogDist.Ratio三条对数度量特征。特征计算公式如下：

公式中的freq.seq为待检测序列的k联碱基组合频率，freq.lnc为lncRNA的平均k联碱基频率，freq.pct为编码RNA的平均k联碱基频率，i表示不同的k联碱基组合，n表示一条待检测序列中的k联碱基频率组合总数。公式中的freq.lnc与freq.pct可通过计算数据库中全部的lncRNA序列与编码区序列的平均k联碱基频率得到。如果编码区序列数量有限，可将编码RNA序列的最长开放阅读框代替编码区序列，计算平均k联碱基频率。在本实例中，使用编码RNA的最长开放阅读框进行平均k联碱基频率的计算。同时在本实例中，针对最长开放阅读框进行对数度量特征计算时，k取值为6；滑动窗口每次向后滑动 3个碱基。若序列中未能找到最长开放阅读框，则针对序列全长提取对数度量特征。此时 k取值为6；滑动窗口每次向后仅滑动1个碱基。

步骤2.3：将由训练集提取而得的五条特征与相应标签另存为特征文件，在本实例中，这五条序列组成特征可取得准确度96.30％与F度量0.9628的表现。

第三步：依据RNA序列的二级结构，设计提取多尺度二级结构特征设计提取多尺度二级结构特征。这一特征大类基于RNA序列的最小自由能与多尺度二级结构序列。提取得到的序列，以10折交叉验证的形式进行特征选择。

图2为6条多尺度二级结构序列的构造方法。这6条多尺度二级结构序列将用以提取多尺度二级结构特征。具体的序列构造方式与特征提取方法如下：

步骤3.1：计算RNA序列的二级结构序列，并用点括号形式表示二级结构序列。此处使用seq[n](seq[n]∈{a，c，g，u})表示一条长度为N的RNA序列，使用SS[n](SS[n]∈{.，(，)})表示序列seq[n]的点括号形式的二级结构序列。RNA序列的二级结构序列可通过程序RNAfold计算而得，也可通过数据库检索，实验验证等其他途径获得。本实例使用程序RNAfold进行二级结构的计算。

步骤3.2：计算RNA序列的最小自由能，并以此作为多尺度二级结构特征中的初级特征。RNA序列的最下自由能可通过程序RNAfold计算获得。

步骤3.3：根据SS[n]，构建如下三条中级尺度的二级结构序列。首先根据SS[n]判断如下四种二级结构子单元：茎stem，环loop，发卡环hairpin与突环bulge。将seq[n]碱基用对应的结构子单元进行替换后，即可得到一条由结构子单元组成的序列，命名为SSE.Full Seq。进一步将序列中相同且连续的结构子单元，使用一个结构子单元来代替，便可得到一条新的结构子单元序列，命名为SSE.Abbr Seq。根据二级结构序列SS[n]中的配对情况，可通过如下公式得到第三条中级尺度的二级结构序列，即配对-不配对序列(Paired-Unpaired Seq)。公式表示形式如下:

通过本步骤，最终得到SSE.Full Seq，SSE.Abbr Seq与Paired-Unpaired Seq三条中级尺度二级结构序列。

步骤3.4：根据seq[n]与SS[n]，构建如下三条高级尺度的二级结构序列。首先将RNA 序列中对应于SS[n]的不配对碱基使用字母D(dot)来代替，得到二级结构序列acguDSeq。利用公式可表示为：

而后，将RNA序列中对应于SS[n]的配对碱基使用字母S(stem)来代替，得到二级结构序列acguS Seq。公式即：

最后将RNA中对应于SS[n]的配对碱基使用相应的大写字母表示，得到二级结构序列 acgu-ACGU Seq。可用公式表示为：

本步骤最终得到acguD Seq，acguS Seq与acgu-ACGU Seq高级尺度二级结构序列。

步骤3.5：针对步骤3.4与3.5中提取的六条中高级尺度二级结构序列，分别使用k-merscheme以及步骤2.2中的对数度量方案来提取相关特征。k-merscheme特征计算公式如下：

k-mer frequencies表示k-mer特征，k为k联碱基组合，c表示某种碱基组合的出现次数，i表示某种碱基组合，l为序列长度。在本实例中，k值由1取值至准确度发生下降时的k值；滑动窗口扫描碱基组合时，每次均向后滑动1个元素。

步骤3.6：针对多尺度二级结构特征与步骤2.3中得到的序列特征组合，并一起进行特征选择。在特征选择阶段，以特征大类为基本的特征选择单位，针对同一多尺度二结构特征序列使用k-merscheme与对数度量方案提取的特征子类，选择能够提升基于序列组成特征的表现的特征子类作为最优的二级结构特征。例如针对高级尺度二级结构序列acguDSeq，利用k-merscheme提取的特征在k取值为1至3时，取得最优准确度86.17％，因此选择此时的k-mer特征与序列特征组合，进行评估。而使用对数度量提取的特征在k 取值为4时，取得最优准确度84.37％，则选择此时的对数度量特征与序列特征组合，进行评估。而评估结果显示，k-mer特征与序列组成特征相组合的准确度为95.87％，未能提高原有序列组成特征的表现；特征与序列组成特征相组合的准确度为96.50％，提升了原有特征的表现。因此选择acguDSeq的对数度量特征作为多尺度二级结构的特征子类。对于其他几种多尺度二级结构序列也以类似的方法进行特征评估与选择。本实例中几种多尺度二级结构特征的具体特征选择结果如下所示：

序列组成特征与多尺度二级结构k-mer特征的特征组合表现

序列组成特征与多尺度二级结构对数度量特征的特征组合表现

在对每一种多尺度二级结构序列进行如上操作后，可以确定k-merscheme与对数度量方案，哪种特征提取策略能够提升原有序列组成特征的表现。但当两种策略提取的二级结构特征均不能提高原有特征表现时，则对k-mer特征子类中权重最高的单条特征进行验证，并将可以提升原有序列组成特征表现的特征选取为多尺度二级结构特征。

例如在本实例中，无论中级尺度二级结构序列Paired-Unpaired Seq的k-mer特征还是对数度量特征均不能提升原有特征的表现。因此进而判断k-mer特征子类中权重最高的单条特征的表现。在k取值为1至4时，Paired-Unpaired Seq的k-mer特征取得最佳表现。本实例中，基于Paired-Unpaired Seq的前十名k-mer特征的重要程度如下所示：

Paired-Unpaired的前十名k-merscheme特征重要程度

UP频率特征的权重最高，意味着此特征发挥着最大的作用，因此将这一特征也纳入多尺度二级结构的候选特征集。在本实例中，多尺度二级结构的候选特征集包括最小自由能，UP频率特征，以及利用对数度量方案提取而得的基于多尺度二级结构序列acguDSeq，acguSSeq和acgu-ACGUSeq的对数度量特征。

候选特征集是初步得到的最有可能提升分类器表现的特征子类或特征的组合，但仍需进行最终的特征评估，以确定最佳的特征组合。本实例中，关于候选特征集的特征选择结果如下所示：

多尺度二级结构候选特征集的特征选择结果

根据特征选择结果，序列特征与多尺度二级结构特征共同可取得最优准确度96.79％。

步骤3.7：将由训练集提取而得的多尺度二级结构特征与相应标签另存为特征文件。在本实例中，多尺度二级结构特征包括最小自由能，UP频率，以及多尺度二级结构序列acguDSeq与acgu-ACGUSeq的对数度量特征。这八条多尺度二级结构特征可取得88.53％的准确度。

第四步：转换基因序列为数值序列，设计并提取基于电子-离子互作赝势的理化性质特征；这一特征大类基于电子-离子互作赝势的能量谱。提取得到的序列，以10折交叉验证的形式进行特征选择。而后判断第二步到第四步中提取而得的特征中是否有冗余特征。如果有，重新进行特征选择。

在筛选得到多尺度二级结构特征后，本发明进入第四步，以提取理化性质特征。这一类基于电子-离子互作赝势与快速傅里叶变换，其分辨lncRNA的原理为编码RNA经快速傅里叶变换后的能量谱，趋向于在三分之一处出现波峰。最早这一性质被用于外显子定位，而本发明则将其应用于lncRNA的鉴定当中，并以此设计提取特征。具体流程包括：

步骤4.1：使用电子-离子互作赝势将训练集中的序列转换为数值序列。其中碱基的电子-离子互作赝势分别为：{a→0.1260，c→0.1340，g→0.0806，t→0.1335}。

步骤4.2：将得到的数值序列使用如下公式，进行快速傅里叶变换：

S_e[kl＝|x_e[k]|²

公式中X_e[n]步骤4.1中得到的数值序列，n为序列长度，{S_e[k]}为经由傅里叶变换后得到的能量谱。

步骤4.3：计算能量谱1/3处波峰的能量值(S_e[N/3])、能量谱平均能量值以及信噪比(SNR)。并依次为基于理化性质的前三条特征。以及SNR的计算公式如下：

图3与图4分别为编码RNA与lncRNA的电子-离子互作赝势能量谱示意图。波峰通常出现在能量谱的1/3处，为进一步保证本发明的容错性，也可判断以能量谱的1/3处为中心的前后各2个位置，共5处的能量值，并选择最大的能量值作为能量谱1/3处的波峰值，即特征S_e[N/3]。

步骤4.4：将能量谱能量降序排列。并在得到的降序排列能量谱上，在不同的区间上提取能量的最大值(Max)，最小值(Min)，第一四分位数(Q1)，中位数(Q2)及第三四分位数(Q3)作为特征。区间可通过如下两种设计方法确定，第一种区间设计方法基于能量谱长度，以降序排列能量谱的前10，前10，前30，…，前100位置，设计区间长度；第二种区间设计方法基于能量谱百分比，以降序排列能量谱的前10％，前20％，前30％，…，完整能量谱，设计区间长度。此步骤在共20种不同区间上分别计算5条特征，最终得到 100条特征。而后对20种区间上的表现进行验证，判断最佳区间。

图5展示了本实例中，不同区间长度上基于理化性质的分位数统计量特征的表现。区间前10％提取而得的分位数特征取得最佳表现，最佳准确度为84.14％。因此以区间前10％为范围，提取分位数特征。

步骤4.5：将步骤4.3与步骤4.4中得到的8条特征进行特征选择，并选取最优表现的特征组合作为最终的理化性质特征。本实例，理化性质特征的特征选择结果如下：

理化性质特征的特征选择结果

根据以上实验，最终筛选得到S_e[N/3]，SNR以及分位数Q1，Q2，Min与Max在内的共6条理化性质特征，且并最终取得准确度88.53％。

步骤4.6：将由训练集提取而得的理化性质特征与相应标签另存为特征文件。

步骤4.7：将第二步到第四步中得到的特征进行组合，并验证特征的表现，以判断特征集中是否具有冗余特征。如果有冗余特征，重新进行特征选择。本实例中的三类特征大类中无冗余特征，并所有特征最终可取得准确度96.84％。本实例中关于三类异源特征的特征组合表现如下：

三类异源特征的综合表现

截至第4步，本发明最关键的特征提取部分已经结束。本实例得到的三类关键的异源特征包括：5条序列组成特征(最长开放阅读框长度与覆盖率，基于最长开放阅读框的对数度量特征)，8条多尺度二级结构特征(最小自由能，UP频率，基于多尺度二级结构序列acguD-Seq与acgu-ACGUSeq的对数度量特征)，与6条理化性质特征(S_e[N/3]，SNR，降序能量谱10％区间的Q1，Q2，Min与Max能量值)。

第五步：构建分类器，用以lncRNA鉴定；此步骤中，可以选择多种机器学习算法进行建模，且参数可以由10折交叉验证来判断。最终选择表现最优的机器学习算法构建最终的分类模型。具体步骤包括：

步骤5.1：在本实例中，三类异源特征最终并将最终特征与相应标签另存为特征文件。

步骤5.2：将步骤5.1中得到的特征用于分类器构建。并通过10折交叉验证的形式对分类器进行调优。在此步骤中，可以选用多种机器学习算法构建分类器，评估结果，并选取表现最优的机器学习算法训练最终的分类模型。特征的设计对lncRNA鉴定的表现起到至关重要的作用，同时也是本发明的关键。优秀的特征在不同的机器学习算法下的表现通常较为稳定，准确度不会有太大波动。在本实例中，logistic回归，支持向量机，随机森林，极限学习机与深度学习，这几种流行的机器学习算法被用于分类器构建；且支持向量机取得最佳表现，10折交叉验证的平均准确度为96.87％。但五种机器学习算法得到的分类器，在准确度上的差异并不明显，证实了本发明提取得到的特征具有良好的稳定性。几种机器学习算法构建而得的分类器具体表现如下：

使用不同机器学习算法构建分类器取得的表现

步骤5.3：得到最优分类器。并用于lncRNA的鉴定。

在详细解释各个步骤后，以下用实验结果证明本发明的有效性。

使用目前最为权威的数据库GENCODE与Ensembl中收集而得的人类(Homosapiens)，小鼠(Mus musculus)，斑马鱼(Danio rerio)与原鸡(Gallus gallus)序列构建测试集，来评估当前流行方法CPC，CNCI，PLEK，CPC2与本发明的表现。这三种物种的测试集中，任意一条序列与步骤1.3中训练集中的序列均不重复。

人类测试集(GENCODE)包括2500条lncRNA与相同数量的编码RNA。以下为四种方法与本发明在人类测试集上的表现：

CPC，CNCI，PLEK，CPC2与本发明在人类测试集上的表现

小鼠测试集(GENCODE)包括1800条lncRNA与相同数量的编码RNA。以下为四种方法与本发明在小鼠测试集上的表现：

CPC，CNCI，PLEK，CPC2与本发明在小鼠测试集上的表现

斑马鱼测试集(Ensembl)包括4000条lncRNA与相同数量的编码RNA。值得注意的是，这一测试集中的一些序列，因测序信号较弱而存在难以确定的碱基(即用X表示)。而对于包含此类序列的测试集，CPC无法进行计算，而CNCI则会自动忽略这类序列。在本实例中，PLEK，CPC2与本发明可以计算斑马鱼测试集序列，而CPC无法计算测试集中的编码RNA序列。此外，CNCI自动忽略了斑马鱼数据集中的13条编码RNA序列。以下为四种方法与本发明在斑马鱼数据集上的表现：

CPC，CNCI，PLEK，CPC2与本发明在斑马鱼测试集上的表现

原鸡测试集包括8000条lncRNA与相同数量的编码RNA。在这一测试集中，同样有些序列，因测序信号较弱而存在难以确定的碱基。CPC，PLEK，CPC2以及本发明可以计算本测试集中的序列，但CNCI忽略了本测试集中的7条lncRNA序列与6条编码RNA序列。

CPC，CNCI，PLEK，CPC2与本发明在原鸡测试集上的表现

从以上评估结果可以发现，本发明有着最佳的准确度，F-度量与跨物种稳定性。在人类数据集中，CPC2作为当前流行的四种方法中最新被提出的方法，取得了最佳的准确度，96.14％。而本方法准确度高达97.28％，准确度更是领先方法CPC超过14％。虽然CPC2 在人类数据集上超过CNCI，PLEK与CPC，但在小鼠数据集上，CPC2仅取得准确度86.11％的表现。本发明在小鼠数据集上准确度达93.47％，领先第二名CNCI超2％，更领先PLEK 近12％。在斑马鱼数据集上，有一部分序列CPC与CNCI无法进行计算，而PLEK，CPC2 与本发明有着较好的容错度，可以顺利完成鉴定。实验结果显示，CNCI的表现优于PLEK 与CPC2，而本发明以88.26％的准确度超过其他方法取得最佳表现。在最后的原鸡测试集上，PLEK以准确度92.35％的表现超过CNCI，CPC2与CPC。CPC虽然能够正常计算这一测试集的序列，但其准确度仅78.36％。而本发明则以准确度94.06％的准确度超过了其他几种方法。

可以发现，当前流行的四种方法在不同数据集上，表现的波动十分明显。PLEK随在原鸡测试集上优于CNCI，CPC2，但PLEK在小鼠与人类测试集上，均劣于CNCI，CPC2。 CPC2在人类数据集上，优于除本发明外的其他几种工具，但在其他几种测试集上的表现均不优秀。综合而言，CNCI的跨物种稳定性优于PLEK，CPC与CPC2，但CNCI与CPC 对序列的容错性较低。研究人员在对某一物种的序列进行计算时，无法提前知晓哪一种方法能够取得最佳表现，因此高准确度与高稳定性为lncRNA鉴定方法是否可行的关键因素。本发明在以上几种物种的测试集上，不仅准确度明显优于当前流行的几种方法，同时更有着优秀的跨物种稳定性与序列容错性。

下面将对当前流行方法CPC，CNCI，PLEK，CPC2与本发明进行效率评估。四种方法与本发明在统一硬件环境下，对5000条序列(2500条lncRNA与2500条编码RNA)进行鉴定，各自的计算时间如下：

CPC，CNCI，PLEK，CPC2与本发明的效率评估

根据效率评估结果，本发明领先于方法CPC，CNCI与PLEK，但劣于CPC2。但本发明仍然能够在1分钟内计算数千条序列，并取得优于CPC2等方法的准确度。

本发明提出的lncRNA鉴定算法，易于使用，方便快捷。本发明能够取得优于当前流行算法与最新算法的准确度与跨物种稳定性，可广泛应用于不同物种的lncRNA鉴定领域，并极大便利注释信息有限的物种研究。此外，本发明不需要进行序列比对，亦有良好的序列容错性。因此相比CPC等基于序列比对的方法更加高效快捷，并相比CPC与CNCI等对序列质量要求较高的方法，可以更加稳定地应用与冷门物种的lncRNA序列鉴定任务。

以上所述为本发明的较好实施例，并不用以限制本发明，凡在本发明精神和原则之内，所进行的任何修改，等同替换，改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法，其特征在于：包括以下步骤：

第一步：构建训练集；

第二步：设计提取基因的序列组成特征；

第三步：依据RNA序列的二级结构，设计提取多尺度二级结构特征；

第四步：转换基因序列为数值序列，设计并提取基于电子-离子互作赝势的理化性质特征；

第五步：构建分类器，用以lncRNA鉴定。

2.根据权利要求1所述的一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法，其特征在于：第一步所述的构建训练集主要进行高质量训练集的构建，具体步骤包括：

步骤1.1：获取某一物种大量的lncRNA序列与编码RNA序列；

步骤1.2：对序列进行质量筛选，即清除质量不佳的序列，以得到高质量lncRNA与编码RNA序列；质量不佳序列通常表现为含有非标准碱基的序列，序列间相似度大于40％的序列。

步骤1.3：选取相同数目的lncRNA与编码RNA构建训练集。

3.根据权利要求1所述的一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法，其特征在于：第二步所述的序列组成特征共五条，设计提取序列组成特征，具体步骤包括：

步骤2.1：提取最长开放阅读框的长度与覆盖率；

对训练集中的序列查找最长开放阅读框，并计算最长开放阅读框的长度；覆盖率为最长开放阅读框的长度与序列长度的比值；

步骤2.2：提取三条k联碱基组合的对数度量特征；

提取对数度量特征，首先需按如下公式定义对数度量方案，用以提取LogDist.LNC，LogDist.PCT与LogDist.Ratio三条对数度量特征，这一对数度量方案不仅可用于本步骤提取序列组成特征，亦可用于阶段3中的基于结构信息的特征的提取，具体公式如下：

公式中的freq.seq为待检测序列的k联碱基组合频率，freq.lnc为lncRNA的平均k联碱基频率，freq.pct为编码RNA的平均k联碱基频率，i表示不同的k联碱基组合，n表示一条待检测序列中的k联碱基频率组合总数；此步骤针对步骤2.1中提取而得的最长开放阅读框提取三条对数度量特征，公式中的k值可根据训练集数据，通过10折交叉验证的形式确定；

步骤2.3：将由训练集提取而得的五条序列组成特征与相应标签另存为特征文件。

4.根据权利要求1所述的一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法，其特征在于：第三步所述的设计提取多尺度二级结构特征，具体步骤包括：

步骤3.1：获取训练集RNA序列的二级结构序列，并用点括号的形式表示；

此处使用seq[n](seq[n]∈{a，c，g，u})表示一条长度为N的RNA序列，使用SS[n](SS[n]∈{.，(，)})表示序列seq[n]的点括号形式的二级结构序列；

步骤3.2：计算多尺度二级结构特征中的初级尺度特征；

这一特征为RNA序列的最小自由能；

步骤3.3：构建三条中级尺度的二级结构序列，以供后续提取中级尺度二级结构特征；

根据步骤3.1中定义的SS[n]，三条中级尺度的二级结构序列构造方式如下：

首先根据步骤3.1中定义的SS[n]判断如下四种二级结构子单元：茎stem，环loop，发卡环hairpin与突环bulge；将seq[n]中的碱基用对应的结构子单元进行替换后，即可得到一条由结构子单元组成的序列，命名为SSE.Full Seq；

进一步将序列中相同且连续的结构子单元，使用一个结构子单元来代替，便可得到一条新的结构子单元序列，命名为SSE.Abbr Seq；

最后，根据二级结构序列步骤3.1中定义的SS[n]中的配对情况，可通过如下公式得到第三条中级尺度的二级结构序列，即Paired-Unpaired Seq:

通过本步骤，可最终得到SSE.Full Seq，SSE.Abbr Seq与Paired-Unpaired Seq三条中级尺度二级结构序列；

步骤3.4：构建三条高级尺度的二级结构序列，以供后续提取高级尺度二级结构特征；

根据步骤3.1中定义的seq[n]与SS[n]，三条高级尺度的二级结构序列构造方式如下；

首先将RNA序列中对应于SS[n]的不配对碱基使用字母D，即dot代替，得到二级结构序列acguD Seq，可用公式表示为：

而后，将RNA序列中对应于SS[n]的配对碱基使用字母S，即stem来代替，得到二级结构序列acguS Seq，公式即为：

最后，将RNA中对应于SS[n]的配对碱基使用相应的大写字母表示，得到二级结构序列acgu-ACGU Seq，利用公式可表示为：

步骤最终可以得到acguD Seq，acguS Seq与acgu-ACGU Seq三条高级尺度二级结构序列，以供后续进行高级二级尺度特征的提取；

步骤3.5：针对步骤3.3与3.4中设计的共六条中级与高级尺度二级结构序列，提取中级与高级尺度二级结构特征，特征可使用k-mer scheme以及步骤2.2中定义的对数度量方案来提取相关特征，k-mer scheme特征的计算公式如下：

公式中，k-mer frequencies表示k-mer特征，k为k联碱基组合，c表示某种碱基组合的出现次数，i表示某种碱基组合，l为序列长度，在使用k-mer scheme以及对数度量方案进行特征提取的过程中，公式中的k值可根据表现评估的结果而确定；滑动窗口扫描碱基组合时，每次均向后滑动1个元素；

步骤3.6：针对多尺度二级结构特征进行特征选择，将准确度最高的特征组合作为最终的多尺度二级结构特征；

步骤3.7：将由训练集提取而得的多尺度二级结构特征与相应标签另存为特征文件。

5.根据权利要求1所述的一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法，其特征在于：步骤4提取基于理化性质的特征，具体步骤包括：

步骤4.1：利用电子-离子互作赝势将训练集中的序列转换为数值序列，其中碱基的电子-离子互作赝势分别为：{a→0.1260，c→0.1340，g→0.0806，t→0.1335}，转换得到的数值序列可用X_e[n]表示；

步骤4.2：将步骤4.1中的数值序列通过快速傅里叶变换转换为能量谱，以供后续进行理化性质特征的提取，快速傅里叶变换可通过如下公式进行：

公式中X_e[n]步骤4.1中得到的数值序列，n为序列长度，{S_e[k]}为经由傅里叶变换后得到的能量谱；

步骤4.3：此步骤可提取三条理化性质特征；

特征分别为能量谱1/3处波峰的能量值，即S_e[N/3]；能量谱平均能量值，即以及信噪比，即SNR。相关计算公式如下：

步骤4.4：此步骤进一步设计并提取其他理化性质特征；

首先将能量谱能量按降序进行排列，并在得到的降序排列能量谱上，按照两种区间设计方法，每种区间设计方法中选取10种不同长度的区间，并针对每种区间分别提取五条理化性质特征；

第一种区间设计方法基于能量谱长度，以降序排列能量谱的前10，前10，前30，…，前100位置，共10种长度以设计区间长度；

第二种区间设计方法基于能量谱百分比，以降序排列能量谱的前10％，前20％，前30％，…，完整能量谱，共10种长度以设计区间长度；

五条理化性质特征分别为不同区间长度上的分位数统计量，包括最大值Max，最小值Min，第一四分位数Q1，中位数Q2及第三四分位数Q3；此步骤在共20种不同区间上分别计算5条特征，而后针对每种区间上的五条特征进行评估，以确定在何种区间范围内得到的五条理化性质特征可取得最高准确度，并将该区间上的5条分位数统计量确定为此步骤得到的最终理化性质特征；

步骤4.5：将步骤4.3与步骤4.4中得到的共8条特征进行特征选择，并选取准确度最高的特征组合作为最终的理化性质特征；

步骤4.6：将由训练集提取而得的理化性质特征与相应标签另存为特征文件；

步骤4.7：将2至4阶段得到的特征进行组合，并进行10折交叉验证，以评估特征表现并判断是否存在冗余特征，当准确度无法继续提升或开始降低时，可判断为特征中无冗余特征，若特征集合中存在冗余特征，则重新进行特征选择。

6.根据权利要求1所述的一种基于序列组成，结构信息及理化特征的lncRNA鉴定方法，其特征在于：步骤五所述的构建分类器具体包括以下步骤：

步骤5.1：将步骤二至四阶段得到的基于序列组成、结构信息与理化性质的三大类特征与相应标签另存为特征文件；

步骤5.2：将步骤5.1中的特征用于分类器构建，并通过10折交叉验证的形式对分类器进行调优；

步骤5.3：得到最优分类器，并用于lncRNA的鉴定。