多物种特征选择及鉴定未知基因的方法

摘要

本发明公开一种多物种特征选择及鉴定未知基因的方法，属于生命科学领域。所述多物种特征选择的方法，包括对覆盖全基因组的小片段区域进行特征赋值和贴注标签处理及物种内、物种间特征选择部分。本发明依靠整合不同物种间的基因共性来构建高效、准确的计算方法，用于准确鉴定和描述未知基因。

说明书

技术领域

本发明涉及生命科学领域，特别是指一种多物种特征选择及鉴定未知基因的方法。

背景技术

目前已发表了多个预测蛋白编码转录本概率的工具，包括CONC、CPC、PhyloCSF、RNAcode、PLEK、CNCI、CNCTDiscriminator、CPAT、HMMER和lncRNA-ID(1-10)等，但是这些工具的绝大部分只使用了转录本的序列信息。这些序列信息包括但不限于：开放阅读框(Open reading frame,ORF)特征，如ORF长度和覆盖率等(1,2,4,7,9)；碱基频率(nucleotide frequencies)特征，如k-mer序列模式、密码子使用频率(codon usage)等(1,2,5,7-9)；保守性分数(conservation score)特征如碱基序列比对或蛋白序列比对等(1-4)；进化相关特征如碱基替换频率(substitution rate)和系统发生树分数(phylogenic score)(7,10)以及模拟性特征(in silico features)等,如预测RNA二级结构和核糖体脱离分数(1,2,7)。

然而，这些序列得到的特征分数存在以下诸多限制，使预测出的数据准确度低。开放阅读框长度(ORF length)：蛋白编码基因中存在的ORF长度远远长于随机产生的长度，因此一个相对较长的ORF(比如大于300nt)的存在可以作为一个较高编码潜能的指标。因为ORF的长度会随着转录本的长度增长而增加，部分的编码潜能分数(4,9)工具会计算ORF长度占转录本的比例。但是需要注意的是，ORF长度本身并不具有绝对的预测效力，比如，多个著名的长链非编码RNA，Xist,Meg3,Hotair,Kcnq1ot1和H19等均含有长度超过300nt的ORF。此外，这个特征要求拼接得到的转录本必须是完整的，但完整的转录本的拼接需要额外的实验的证明。

碱基、密码子或短文字频率(Nucleotide,codon or short word frequen-cies)：在蛋白编码基因中的碱基频率并不是随机的，这一点可以用来区分蛋白编码和非编码基因。现有的很多个编码潜能分数工具都使用了这个特征。而且由于其计算量低的特点，很多的工具都很依赖这类特征。需要注意的是，这类特征依赖转录本的长度，如果转录本的长度较短的话，预测结果就存在偏差。较长的转录本可以提供更加准确地评估碱基频率。

替换模式(substitution patterns)：蛋白编码基因的序列为了保持开放阅读框并保留特定的氨基酸序列而承受进化压力。这种进化压力的存在会体现在多基因组序列比对文件中：在同一个阅读框中的不同位置的替换频率不一致、在开放阅读框中的插入和剔除事件频率很低、即使发生插入和剔除时氨基酸序列也会得以保留。但是这类特征在进化上刚刚产生功能的长链非编码RNA并没有用处。对于没有多重比对序列的未知的长链非编码RNA也没有用处。

已知的蛋白功能域的存在：蛋白编码基因基本上都会包含一些常用的功能域，一些概率模型就记录了这些功能域，从而可以搜索一条未知类型的转录本中是否含有记录的功能域。这类特征的固有缺点也比较明显，即有些蛋白编码基因本身也不具有特别保守的功能域，也会被预测为长链非编码RNA，而存在预测结果失误的问题。

上述这些编码潜能分数计算工具多数没有考虑到实验数据的整合，而仅仅利用了转录本序列信息。幸运的是，很多的实验数据提供了鉴定非编码RNA的额外信息。比如，多数的蛋白编码基因会富集在有polyA尾的RNA测序数据中，而经典的非编码RNA则不具有polyA尾。此外，长链非编码RNA在不同组织中具有更高的表达特异性。还有Ribosome profil-ing数据提供了全基因组翻译过程的快照。

得益于高通量测序方法的迅猛发展，生物学家发现了基因组中存在大量未知的转录本，这其中的很大一部分是新型的非编码RNA(noncoding RNA)。非编码RNA在生物体内行使着复杂且精细的调控功能，其相关功能的研究吸引了众多生物学家的广泛兴趣，但不同于蛋白质编码基因，我们对非编码RNA的表达调控缺乏综合性的了解，不同非编码RNA的生物物理特性也不尽相同，这使得在物种内以及不同物种间鉴定非编码RNA仍缺乏一个统一的鉴定标准，高通量的鉴定所有类型的非编码RNA依旧是一个极富挑战的问题，相应的计算方法亟需发展。

发明内容

本发明要解决的技术问题是提供一种即可以跨物种预测又高效、准确的多物种特征选择及鉴定未知基因的方法。

为解决上述技术问题，本发明提供技术方案如下：

一方面，提供一种多物种特征选择的方法，包括如下步骤：

步骤1：选取不同物种全基因组区域并将其切割成100nt长的基因组小片段，在基因组小片段的基础上计算结构和序列、表达水平、组蛋白修饰水平和转录调控因子结合水平、上下游影响值的特征值，并对上述特征值进行归一化处理；

步骤2：根据基因组小片段所在的基因元件对每一个基因组小片段贴注释标签，用随机森林算法作为分类器，对基因组小片段进行分类，并根据不同注释标签将基因组小片段分到对应的特征集内，其中选取经典的非编码RNA、确定的蛋白编码区域、5’和3’端非翻译区域和负对照区域(即表达水平极低的基因间区)4种基因元件做为黄金标准集；

步骤3：物种内的特征选择，以随机森林重要性指标作为排序指标，对特征集内基因组小片段进行排序，使用递归特征剔除算法对物种内的选择特征集进行预筛选来去掉非必需的特征，使用贪婪后向算来进一步筛选物种内特征集；

步骤4：物种间的特征选择，将不同物种的特征集取并集和交集，并集去掉交集剩余的特征构成了补集，将补集中误删除的必需特征添加到交集特征集中，得到最终的共有特征集。

其中，步骤1中，提取了622套高通量数据。这些数据来自于5个不同的物种，即人类(Homo sapiens,h19),小鼠(Mus musculus,mm10),线虫(Caenorhabditis elegans,ce10),果蝇(Drosophila melanogaster,dm3)和拟南芥(Arabidopsis thaliana,TAIR10)。对全基因组每100nt的区域，即每一个基因组小片段(bin)，本发明不光从大量高通量数据中计算了其表达水平、组蛋白修饰水平、转录调控因子结合水平的信号强度，还为其计算了GC含量、DNA序列保守性、蛋白质序列保守性、RNA二级结构稳定性、同源RNA二级结构、RNA二级结构的保守性、开放阅读框特性的特征值。由于基因组本身的特性，对每一个基因组小片段，本发明还考虑了其相应特征值的上下游信号，尤其是组蛋白修饰水平、转录调控因子结合水平信号值的上下游信号。

其中，所述步骤1中，GC含量(GC content)是指序列中的碱基G和C的含量比例；DNA序列保守性(DNA sequence conservation)是使用软件计算或者下载得到的；对于人类数据，下载了UCSC数据库里的phastCons46-way的PhastCons scores；对于小鼠数据，下载了UCSC数据库里的phastCons30way的PhastCons scores；对于果蝇数据，使用BLASTn的默认参数在Flybase(包含11种其他的果蝇物种)数据库里搜索计算得到；对于线虫数据，使用BLASTn的默认参数在Wormbase(包含19种其他的线虫物种)数据库里搜索计算得到；对于拟南芥数据，使用BLASTn的默认参数在EnsemblePlants(包含31种其他的植物物种)数据库里搜索计算得到；蛋白质序列保守性(protein conservation)的计算使用了和DNA序列保守性相同的物种，不同的是，使用了NCBI的nr数据库里的蛋白质序列来建BLAST库，对所有物种本发明均使用了BLASTx来计算保守性数值，除了线虫，本发明使用了tBLASTx软件；RNA二级结构稳定性(RNA secondary structure stability)是通过Randfold软件计算的结构的自由能的p-value来表示的；为了计算p-value，随机模拟了1000次碱基对出现的频率来作为背景噪音信号；同源RNA二级结构(RNA secondary structure homologs)是通过使用INFERNAL软件的默认参数来搜索Rfam得到的0-1值；RNA二级结构的保守性(RNA secondary structure conservation)是由RNAz软件在默认参数下计算各物种的多重比对序列得到的结构保守性指标structure conservation index(SCI)表示的；ORF特性(ORF property)是通过输入各物种的多重比对序列到RNAcode软件中，使用默认参数得到的。

在UCSC数据库上下载了人类和46个脊椎动物物种的多重比对序列、小鼠和30个脊椎动物物种的多重比对序列、线虫和7个线虫物种的多重比对序列和果蝇和15个果蝇物种的多重比对序列。VISTA：

(http://pipeline.lbl.gov/downloads.shtml)上下载了拟南芥和5个植物物种的多重比对序列。

表达水平数据(expression data)包括RNA深度测序数(RNA-sequenceing)和覆瓦式阵列数据(Tiling array)。本发明从98套有polyA尾的RNA深度测序(poly(A)+RNA-seq)、41套无polyA尾的RNA深度测序(poly(A)-RNA-seq)、48套总RNA深度测序(total RNA-seq)和70套小RNA深度测序(small RNA-seq)数据中获得了160亿测序读段(reads)数据。并从101套有polyA尾的RNA覆瓦式阵列(poly(A)+tiling array)或总RNA覆瓦式阵列(total RNA tiling array)数据中获得了4亿探针(prob)数据。

对于RNA深度测序数据，用DEGseq软件计算了每个基因组小片段的RPKM值(reads per kilobase per million)来代表表达水平。对于覆瓦式阵列数据，使用了R语言软件包AffyTiling来计算在每一个基因组小片段上重合的探针强度的最大值。最近的实验证据表明，一些新的非编码RNA更加倾向于在特定的组织中特异性表达。为了能够更加灵敏地检测到这种特异表达的新型非编码RNA，使用了不同实验条件下(如不同的细胞系、组织、发育阶段，生长环境等)的最大值来代表表达水平。对于相同条件下的重复试验，本发明取了平均值。

人类的表达水平数据来自于ENCODE项目，包含5个不同的细胞系：K562、GM12878、H1-hESC、HeLa-S3、和HepG2。所有数据均为RNA深度测序数据，包括了全细胞RNA深度测序(whole cell RNA-seq)、细胞核RNA深度测序(nucleus RNA-seq)和细胞质RNA深度测序(cytosol RNA-seq)。根据所富集的RNA的种类，可以分为有polyA尾的RNA深度测序(poly(A)+RNA-seq)，无polyA尾的RNA深度测序(poly(A)-RNA-seq)和小RNA深度测序(small RNA-seq)。每组实验均有两个重复试验。

小鼠的数据来自于ENCODE项目和GEO数据库，包含了不同的组织和细胞系：肝组织、心肌组织、ES-Bruce4细胞系和CH12细胞系。所有数据均为RNA深度测序数据，根据所富集的RNA的种类，可以分为有polyA尾的RNA深度测序(poly(A)+RNA-seq)、无polyA尾的RNA深度测序(poly(A)-RNA-seq)、总RNA深度测序(total RNA-seq)和小RNA深度测序(small RNA-seq)。

线虫数据来自于modENCODE项目，包含了不同的发育阶段：胚胎(embryo)、L1到L3阶段的幼虫(L1-L3 larvae)和幼体成虫(young adult)。所有数据均为RNA深度测序数据，根据所富集的RNA的种类，可以分为有polyA尾的RNA深度测序(poly(A)+RNA-seq)、总RNA深度测序(total RNA-seq)和小RNA深度测序(small RNA-seq)。

果蝇数据来自于modENCODE项目，包含了不同发育阶段：胚胎(embryo)、L1到L3阶段的幼虫(L1-L3 larvae)、和成虫(adult)。RNA深度测序数据包括有polyA尾的RNA深度测序(poly(A)+RNA-seq)、无polyA尾的RNA深度测序(poly(A)-RNA-seq)、总RNA深度测序(total RNA-seq)和小RNA深度测序(small RNA-seq)。还有一部分总RNA覆瓦式阵列(total RNA tiling array)数据。

拟南芥数据来自于GEO数据库，即拟南芥幼苗在不同生长环境下的测序数据：正常环境、干旱环境或高光环境。RNA深度测序数据包括有polyA尾的RNA深度测序(poly(A)+RNA-seq)、无polyA尾的RNA深度测序(poly(A)-RNA-seq)和小RNA深度测序(small RNA-seq)。还有一部分的覆瓦式阵列数据，包括有polyA尾的RNA覆瓦式阵列(poly(A)+tiling array)和总RNA覆瓦式阵列(total RNA tiling array)数据。

整理了18种不同的组蛋白的修饰数据和2种转录调控因子的结合数据。对于这些数据，用MACS14来处理原始数据得到结合位点的相对强度。紧接着，使用UCSC的工具bigWigAverageOverBed来得到每一个基因组小片段的信号值。每套数据经历两步清理：首先去除了0.01％的异常值，再使用Z-score方法进行标准化(即中心化和归一化)。之后，对来自于不同实验条件下的数据取平均值来代表其特征值。

基因往往受到其上下游的影响，比如，一些组蛋白修饰在特异地富集在基因的启动子区域。因此,对每一个基因组小片段，除了计算其所在位置上的特征值，还考虑了其上下游特征值的影响，即上下游影响值(Context Influence score,CIS)。考虑到上下游的影响会随着距离的增加而衰退，本发明对每一个基因组小片段的上下游均赋值一个以指数形式随距离衰减的权重函数(17):其中，h_k代表当前基因组小片段第k个特征值，d_k代表距离当前基因组小片段的距离，d₀代表衰退距离。

对每一个基因组小片段的每一个特征，其上游影响值，所在位置的特征值和下游影响值总是作为一个整体来考虑。本发明对每一个物种的每一个特征均优化选出了合适的衰退距离参数d₀。

以上不同特征的初始值是范围极大的。为了减少数据异质性的影响，首先平滑处理了数据的最高0.005％和最低的0.005％的异常值。数据范围大于1000的特征经历了对数化处理：y＝log₁₀(x-minimum)+1。所有的特征值均转化其范围为0到1：y＝(x-minimum)/(maximum-minimum)。

其中，所述步骤2中，基因组小片段表达水平低于所有的基因间区的平均值，不含有任何不确定的碱基，并且所在的基因间区是与任何已知基因组元件相隔最少2000nt(对人类和小鼠基因组)或500nt(对果蝇，线虫和拟南芥基因组)远的区域时，则判定该基因组小片段为负对照基因组小片段。如果一个基因组小片段的50％落在已知的非编码RNA上，该基因组小片段就被注释为相应的非编码RNA基因组小片段；2)如果一个基因组小片段的90％落在已知的CDS，UTR，ancestral repeat或intergenic region，那么就会被相应地注释；3)如果一个基因组小片段的50％落在pseudogene,intronic,TE or ambiguous等区域，那么就会被相应地注释。

黄金标准集的基因组小片段均被贴上对应的注释标签用来训练和检验模型。特征集直接决定模型的表现。模型的表现是指模型在检验集中的分类准确度。因此，一种特征选择的原则是：分类准确度越高，模型越好，则特征集越好。基于此原则，使用了基于监督式机器学习(supervised machine learning)算法，并交叉验证模型表现。在这种框架下，优化了构建训练集和检验集的抽样方法，使得训练集和检验集更加独立。

首先将黄金标准集分为2/3的训练集和1/3的测试集。接着，在训练集上使用5倍交叉验证的策略挑选最优算法。考虑了基因组小片段之间的距离的影响。由于邻近的基因组小片段的特征值会高度相关，邻近的基因组小片段不能被分别分配在训练集和测试集中。在机器学习时，如果邻近的两个基因组小片段分别被分配在训练集和测试集里，那么在测试集的基因组小片段会被预测为与其相近的训练集的基因组小片段的标签，这种预测结果是过度拟合(over-fitting)的结果，不具有实用价值。

为了消除这种影响，本发明开发了区块抽样(block sampling)的策略：相近的基因组小片段会合并成基因组区块(genomic block)，而来自相同区块的基因组小片段就会被分配到同样的集合中(训练集或测试集)。为了避免来自于同一个基因中的两个外显子被分配到不同的基因组区块中，本发明选择大于90％的内含子的长度的距离作为基因组区块之间最小的距离。最终，选定在人类和小鼠基因组中，基因组区块之间的最小距离是15000nt；在线虫，果蝇和拟南芥基因组中，基因组区块之间的最小距离是5000nt。

其中，所述步骤3中，递归特征剔除算法从最差的特征开始一个一个剔除，直到剔除到某个特征时，随机森林算法分类器的准确度开始下降即停止剔除，贪婪后向算法可以以一种递归的形式剔除当前特征集中最差的特征，直到当前特征集中剔除任意一个特征均会导致随机森林算法分类准确度变差为止。递归特征剔除算法(RFE)是作为物种内特征选择的预筛选步骤的。通过对特征排序(排序越高，特征越优秀)，RFE算法不可逆地逐步剔除排序较低的特征。在误差允许范围内，RFE算法给出最小的特征集。贪婪后向算法(GBA)是作为物种内特征选择的严格筛选步骤的。贪婪后向算法递归地评估候选特征集，剔除当前候选特征集中最差的特征，更新候选特征集，直到从候选特征集中剔除任意特征均会导致模型变差为止。预先定义的容忍度阈值可以帮助贪婪后向算法跳出一些局部最优解。

由于分类器中使用了随机算法，这导致分类器在处理同一套数据时也会出现准确度的波动性。为了消除掉随机波动性带来的影响，首先通过重复建模100次来估计出分类器分类准确度的分布。该准确度的分布类似正态分布。基于此，将这个随机波动的大小定义为该分布的4倍标准差。在预筛选步骤中，对于两个准确度相差不大(即相差不超过随机波动)的特征集，本发明优先选择较小的特征集，而不是准确度更高的特征集。

其中，所述步骤4中，物种间的特征选择目的是找到在多个物种中均表现很好的特征集。在物种内特征选择之后，四个物种中均有的特征有7个，这7个特征构成了交集特征集intersecting feature set:I＝(I1,…,I7)。交集特征集中包括如下7个特征：GC含量(GC content)、DNA序列保守性、蛋白质序列保守性、RNA二级结构稳定性、同源RNA二级结构、RNA二级结构的保守性和开放阅读框特性。

由于之前的特征选择步骤中均不是全局最优算法，补集特征集中可能含有在之前的特征选择步骤中误删掉的必需特征，所以尝试将第二步中将补集特征集中的部分特征添加到交集特征集中。那么从补集特征集中挑选n个特征，可以有个子集G,其中的一个子集可以记为Gn,j(j＝1,…,m)。将子集Gn,j添加到交集特征集可以构成候选特征集Cn,j(Cn,j＝Gn,j+I)。对每一个物种s(s＝human,mouse,worm or fly)，相同大小的所有的候选特征集Cn,j均会被排序。以添加两个特征(n＝2)为例，本发明可以得到m＝105个子集(Gn,j(j＝1,…,105))，相应的，可以得到105个候选特征集Cn,j。对每个候选特征集Cn,j，选择其在四个物种中倒数第二的准确度作为其准确度，这样可以保证至少有三个物种的模型是优于该准确度的。对所有的候选特征集Cn,j的准确度进行排序，并使用该排序分数衡量候选特征集Cn,j的最终表现。

不论挑选2、3或4个特征到子集G中，DNA序列保守性、H3K36me3和H3K4me3组蛋白修饰特征均会以较高频率出现在排序较高的候选特征集中，故而将它们添加到交集特征集中，得到最终的四个物种中的共有特征集。最终选择得到的共有特征包括：GC含量(GC content)、DNA序列保守性(DNA sequence conservation)、蛋白质序列保守性(protein conservation)、同源RNA二级结构(RNA secondary structure homologs)、ORF特性(ORF property)、有polyA尾的RNA深度测序(poly(A)+RNA-seq)、无polyA尾的RNA深度测序(poly(A)-RNA-seq)、小RNA深度测序(small RNA-seq)、H3K4me3组蛋白修饰(H3K4me3 modification)和H3K36me3组蛋白修饰(H3K36me3 modification)。共有特征集选自人类、小鼠、线虫和果蝇4个物种。

由于基因组序列数量巨大，并且采集样本困难，建议使用高通量高性能的计算集群处理，本发明中使用64GB内存，CPU model:Intel(R)Xeon(R)CPU E5-2609 0@2.40GHz；GNU C(GLIBC):ldd(GNU libc)2.12；GCC version:gcc(GCC)4.4.6 20120305(Red Hat 4.4.6-4)；Red Hat Enterprise Linux Server release 5.5(Tikanga)系统进行数据收集并处理。

另一方面，提供一种基于上述特征选择方法的鉴定未知基因的方法，包括如下步骤：

步骤1：首先构建了人类基因组(Homo sapiens,hg19)的索引文件。全基因组按照染色体(1-22号染色体和XY染色体)和正负链分成了48大块，每一大块均被100nt的小窗按照50nt的步长切割成100nt长的基因组小片段。每一个基因组小片段均包含了其独一无二的位置信息：染色体，正负链，起始和终止位点。当对基因组小片段按照其所在位置进行排序之后(先按照染色体排序，再按照正负链排序，最后按照起始位点排序)，基因组小片段的排序与其所在的位置构成了一一对应的索引关系，此时可以反过来利用基因组小片段的排序找回其所在的基因组位置。这就是构建索引的过程。通过相同的流程，构建了小鼠(Mus musculus,mm10)、线虫(Caenorhabditis elegans,ce10)、果蝇(Drosophila melanogaster,dm3)和拟南芥(Arabidopsis thaliana,TAIR10)的基因组索引文件。

对每一个基因组小片段，本方法首先计算了九个共有特征的特征值(RNAfeature)算法得到的十个共有特征中除去“ORF特性”特征，因为这个特征对于转录本完整性的要求很高，而实际拼接得到的转录本很不精确)。这九个特征值包括GC含量(GC content)，DNA序列保守性(DNA sequence conservation)，蛋白质序列保守性(protein conservation)，同源RNA二级结构(RNA secondary structure homologs)，有polyA尾的RNA深度测序(poly(A)+RNA-seq)，无polyA尾的RNA深度测序(poly(A)-RNA-seq)，小RNA深度测序(small RNA-seq)，H3K4me3组蛋白修饰(H3K4me3 modification)和H3K36me3组蛋白修饰(H3K36me3 modifica-Tion)。由于不同物种特征选择方法中基因组小片段已经完成了索引构建，所以对应的特征值就可以按照基因组小片段的索引，顺序存储在允许随机访问的HDF5格式的文件中。

步骤2：将转录本拆分，获取该转录本的多重特征值的向量；

步骤3：对每个转录本的所有的特征值的向量进行最大值、均值和方差计算，并使用最大值、均值和方差来代表该转录本的特征值，构建整个转录组的特征值矩阵；

步骤4：使用9个共有特征集及上述转录本的特征集，对确定的蛋白编码区域的基因组小片段和长链非编码RNA的基因组小片段建模；由于样本集的不平衡性，本发明采用了平衡的随机森林算法(Balanced Random Forest)。平衡的随机森林算法首先将不平衡的样本集经过多次向下采样，生成多个平衡的小样本集，并保证所有的样本被采样至少一次。对每个平衡的小样本集，应用随机森林算法(Random Forest)训练之后，会生成多个子模型，再将多个子模型使用bagging的方法合并，成为最终的模型。以人类数据为例，每次，随机抽取各类基因组小片段各10000个，并用着20000个基因组小片段训练一个子模型。重复这一过程100次，可以得到100个子模型。将这100个子模型使用bagging的方法合并，可以得到最终模型，其预测结果是100个子模型的平均结果。

步骤5：将未知基因置于上述模型中，采用随机森林分类器算法计算转录本的编码潜能分数，鉴定未知基因。

其中，所述步骤2中，转录本拆分的步骤为先将其外显子拆分，并得到与其外显子区域重合(50％)的基因组小片段的索引，再将其上游区域拆分，得到与其上游区域重合(50％)的基因组小片段的索引，其中对于启动子标记物H3K4me3组蛋白修饰特征，与转录本上游区域重合的索引并通过索引得H3K4me3组蛋白修饰特征的向量；对于其他特征而言，与转录本外显子重合的索引并通过索引得到每个特征的向量。

其中，所述步骤2中，对于人类和小鼠数据中上游区域是5000nt；对于线虫，果蝇和拟南芥数据中上游区域是2000nt。

随机森林分类器算法得出的转录本潜能分数数值范围应置于0～1之间，当大于0.5时，未知基因则为编码蛋白序列，小于0.5时，未知基因则为非编码RNA序列。

由于鉴定未知基因的方法(定义为COME)需要使用多物种特征选择得到的特征矩阵，所以本发明提供了五种模式生物的特征矩阵和COME的可执行程序:http://github.com/lulab/COME。可以在任何类型的计算机上实现，本发明还提供了COME的web server：

http://RNAfinder.ncrnalab.org/COME。

本发明具有以下有益效果：

本发明中，开发了一套多物种特征选择的算法，定义为RNAfeature，来筛选得到一套在多个物种中的多种非编码RNA所共有的特征集。这些特征是由RNAfeature算法在大于600套不同组织、不同细胞系、不同发育阶段的基因组数据、转录组数据和表观遗传数据经过多重筛选得来。RNAfeature算法以机器学习模型来区分全基因组的经典的非编码RNA(canonical ncRNAs)、确定的蛋白编码区域(confirmed coding sequence,CDS)、5’和3’端非翻译区域(untranslated regions,UTRs)和负对照区域(netative control，即表达水平极低的基因间区)，并且考虑了每个基因组区域的局部特征和上下游影响；同时开发了一套基于多重特征鉴定未知基因的方法，定义为COME。通过机器学习模型整合RNAfeature挑选得到的共有特征，COME可以计算出蛋白编码转录本的概率值。由于拼接得到的转录本经常是不完整的，对这些新型长链非编码RNA，COME使用拆分合并策略来消减不完整转录本带来的影响，极大的提高了算法的鲁棒性。

附图说明

图1为本发明的多物种特征选择方法流程图；

图2为本发明的鉴定未知基因的计算方法；

图3跨物种预测非编码RNA的准确度；

图4人类数据中跨种类预测非编码RNA；

图5四个物种中跨种类预测非编码RNA的准确度；

图6使用共有特征描述长链非编码RNA MALAT1中潜在的功能域；

图7 COME算法输入文件界面；

图8 COME算法输出鉴定结果界面。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明针对现有技术中预测蛋白编码的工具由于转录本序列信息不完整预测未知基因片段是否为编码蛋白序列或非编码序列准确度低，存在偏差的问题，提供一种多物种特征选择及鉴定未知基因的方法。

如图1所示，一种多物种特征的选择方法，包括如下步骤：

本发明收集了四个物种的高通量数据、序列和结构特征，并把不同的特征按数字标记为f1，f2，f3等。

图1中使用不同颜色来区分来自不同物种的特征：红色代表人类数据提取的特征，绿色代表小鼠数据提取的特征，紫色代表果蝇数据提取的特征，蓝色代表线虫数据提取的特征。这些特征被用来区分四类基因组元件：确定的蛋白编码区域(confirmed coding sequence,CDS)、经典的非编码RNA(canonical ncRNAs如rRNA、tRNA、miRNA等)、5’和3’端非翻译区域(untranslated regions,UTRs)和负对照区域(即表达水平极低的基因间区(intergenic regions with weak expression signals),negative control)。

特征选择过程包括了物种内特征选择部分：即在同一个物种数据中剔除不能将四种基因组元件区分开来的特征，和物种间特征选择部分：即从落选的特征中找到四个物种共有的特征，并添加到最终的共有特征集中(common features)。在物种内特征选择部分结束之后，所得到7个特征是被所有四个物种选出来的，这个特征集被标记为特征交集(intersect set)。将特征交集对比特征全集可以得到特征补集(remaining set)。通过物种间特征选择步骤，RNAfeature从特征补集中再次选出3个特征，可以帮助提高模型分类能力。

图2所示，基于上述多物种特征选择方法的鉴定未知基因的方法，包括如下步骤：

(a)首先构建全基因组索引文件，即将全基因组切割成固定长度的基因组小片段，然后使用基因组小片段的排序信息来获取其对应的基因组位置信息。同时，对每一个基因组小片段计算多个特征的特征值，从而实现对多重特征的索引构建。

(b)将转录本进行拆分，获取转录本的多重特征值的向量。

(c)通过计算每个特征值向量的最大值、均值和方差实现转录本特征值合并。

(d)COME使用监督式分类器算法将未知基因序列置于整合多重特征的模型中来区分蛋白编码转录本和非编码RNA。预测未知序列为蛋白编码转录本的概率值。

本发明中，整合不同物种间的基因共性来构建高效、准确的计算方法，可以准确鉴定并由于拼接得到的转录本经常是不完整的，对一些新型长链非编码RNA，COME使用拆分合并策略来消减不完整转录本带来的影响，可以极大地提高本发明的鲁棒性。

实施例1

对于每一个物种而言(对角线所在元素位置)，非编码RNA与确定的蛋白编码区域、非翻译区域和负对照区域均可以很好地区分开来。为了验证RNAfeature所选出的共有特征的鲁棒性(robustness)，测试了其跨物种预测的准确性。

如图3所示，跨物种预测非编码RNA的准确度，测量了四个物种(H代表人类、M代表小鼠、F代表果蝇、W代表线虫)模型区分四类基因组元件的准确度(accuracy ACC)。物种内预测的准确度展示在对角线上，跨物种预测的准确度以非对角线元素展示。图中，红色代表确定的蛋白编码区域(CDS)，紫色代表5’和3’端非翻译区域(UTRs)，绿色代表经典的非编码RNA(canonical ncRNAs)，蓝色代表负对照区域(negative control。

跨物种预测是指对特征数据标准化处理之后，使用共有特征构建各个物种模型，并使用物种A的模型去预测其他物种的数据。

如图3非对角线元素所示，所得到的跨物种预测准确度非常高(70％到90％)，这说明RNAfeature所选共有特征是在四个动物物种中普适的。

除此之外，将四个动物物种数据混合，生成了混合模型。接着将混合模型应用到了植物物种拟南芥的数据上。该预测的准确度，即将非编码RNA与确定的蛋白编码区域、非翻译区域和负对照区域区分开来的准确度，可以达到82％。除了准确度之外，以非编码RNA为正样本，其他样本为负样本，还计算了精确度分数为0.84,敏感性分数为0.85,特异性分数为0.95。如此高的准确度预示着RNAfeature所选出的共有特征不只局限在动物物种之内。

实施例2

RNAfeature只使用了经典的非编码RNA(包括rRNA,tRNA,snRNA,miRNA,Y RNA等)作正样本集，所以很有必要检验共有特征预测新型非编码RNA的能力。为此，本发明设计了跨种类的交叉验证，即首先将某一类型的非编码RNA(包括rRNA,tRNA,snRNA,miRNA,Y RNA等)从RNAfeature的训练集中剔除，再使用共有特征训练模型，使用该模型去预测剔除掉的非编码RNA类型。

如图4所示，以人类数据为例，箱线图展示了特定类型的基因组元件(每个窗口的题目)被预测为不同的类别(x轴)时的概率分布(y轴)。不同颜色代表的不同类别：红色代表确定的蛋白编码区域(CDS)，紫色代表5’和3’端非翻译区域(UTRs)，绿色代表经典的非编码RNA(canonical ncRNAs)，以及蓝色代表负对照区域(negative control)。

在人类数据中，跨种类预测模型可以很好地将任意种类的非编码RNA与确定的蛋白编码区域、非翻译区域和负对照区域区分开来(平均的敏感性分数是0.89)。即使在被剔除掉的非编码RNA和剩余的非编码RNA性质非常不同的情况下，比如剔除掉的非编码RNA是rRNA时，预测的灵敏度分数也高达0.88。

实施例3

对其余物种中也执行了跨种类的交叉验证。如图5所示，展示了4个物种(人类、小鼠、果蝇和线虫)数据中，RNAfeature模型区分确定的蛋白编码区域(CDS)、5’和3’端非翻译区域(UTRs)、经典的非编码RNA(canonical ncRNAs)和负对照区域(negative control)的准确度。NA表示该物种中不含有该类型的非编码RNA。更深的颜色代表更高的准确度。

结果显示，四个物种的跨种类交叉验证的准确度均很高，而且远远高于随机猜测的准确度。相比较而言，人类和小鼠模型的准确度要比果蝇和线虫的模型准确度高一些，这可能是由于人类和小鼠的已经注释的非编码RNA数目更多一些，模型的训练更加充分。

实施例4

以一个功能比较清楚的长链非编码RNA，MALAT1(Metastasis associated lung adenocarcinoma transcript 1)，为例来表明如何使用共有特征来寻找长链非编码RNA中的潜在功能域。

如图6所示，IGB画出非编码潜能分数及共有特征值在MALAT1的其转录本全长上的信号值。各个特征及其缩写是：GC含量(GC％)，DNA序列保守性(DNA Cons)，蛋白质序列保守性(Protein Cons)，同源RNA二级结构(RNA structure)，ORF特性(ORF property)，有polyA尾的RNA深度测序(poly(A)+)，无polyA尾的RNA深度测序(poly(A)-)，小RNA深度测序(small RNA)，H3K4me3组蛋白修饰(H3K4me3)和H3K36me3组蛋白修饰(H3K36me3)。

MALAT1的5’端有很强的H3K4me3组蛋白修饰信号，这表明MALAT1被特定的组蛋白修饰调控，3’端的非编码潜能较高，暗示潜在的功能域可能在其3’端位置。

实施例5

随机选取已知的13段序列使用COME进行鉴定验证，但在使用COME算法鉴定非编码RNA需要实验人员提供转录本信息和特征矩阵信息。由于特征数据的采集的困难，本发明提供了五种模式生物的特征矩阵和COME算法的可执行程序:http://github.com/lulab/COME。

出于用户体验的考虑，本发明还提供了COME的web server：http://RNAfinder.ncrnalab.org/COME。

COME web server的输入文件如图7所示，是转录本的坐标文件(gtf格式)。输出的鉴定结果如图8所示，显示了每一条转录本的编码潜能分数和预测结果。鉴定结果(如图8所示)，准确率达100％。

本发明中，构建了一套多物种特征选择的方法(RNAfeature)，并得到了在多个物种中具有较高的保守性非编码RNA共有的特征。使用这些共有特征，可以准确地跨物种鉴定非编码RNA，包括经典的非编码RNA和长链非编码RNA。为了提升长链非编码RNA鉴定的准确度和鲁棒性，设计了COME，来整合共有特征并计算蛋白编码转录本的概率值。使用这些共有特征，还可以进一步描述非编码RNA，包括对长链非编码RNA进行分类，以提供验证的候选集和寻找长链非编码RNA中潜在的功能区域。

上述的对实施例的描述是为了便于该技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对这些实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于这里的实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。