CN110136776B - 一种从低质量核糖体印迹数据预测基因编码框的方法和系统 - Google Patents

Abstract

本发明公开一种从低质量核糖体印迹数据预测基因编码框的方法，本发明综合利用核糖体印迹和密码子使用频率进行蛋白编码框的预测，利用multitaper算法和复杂度对核糖体印迹数据质量进行描述，根据核糖体印迹数据的复杂度自动分配相应的权重，从而平衡数据质量的影响。具体的，本发明提取密码子使用频率，结合核糖体印迹数据的3碱基周期性，科学度量核糖体印迹的数据质量并合理分配相应权重，计算每个密码子位于核糖体P位点的概率，提取序列特征，通过统计分析综合评定编码框的预测概率值，进而预测新的编码框。本发明将大幅降低对核糖体印迹数据质量的要求，将极大促进核核体印迹技术应用的拓展，特别是在农作物研究中的应用。

Description

一种从低质量核糖体印迹数据预测基因编码框的方法和系统

技术领域

本发明属于生物技术领域，具体涉及利用低质量的核糖体印迹数据进行蛋白编码框的预测方法，即一种从低质量核糖体印迹数据预测基因编码框的方法，还涉及一种预测基因编码框的系统。

背景技术

随着第二代和第三代基因测序的不断发展，基因组数据近年来呈井喷式增长，极大的促进了生命科学的研究和应用。基因功能是一切生命活动的基础，对基因功能的研究有助于增进我们对疾病发生，以及农作物性状形成机理的了解，并进而帮助人们更加有效的预防和治疗疾病或改良农作物性状。在已有的很多基因组学和生物学研究中，人们主要关注基因组中较大的编码基因(长度>＝300bp)，而直接忽略基因组中的小编码框，认为其表达量低，编码能力弱，没有或者仅有非常次要的功能。随着人们对基因组研究和认识的不断深入，越来越多证据表明，基因组中的小编码框在基因表达和翻译中均发挥了重要的调控作用，对植物性状形成、酵母发育以及动物胚胎发育都具有非常关键的作用。由此可见，基因小编码框的研究在医学、工业和农业应用中都具有非常广泛的前景。与此同时，基因小编码框的研究对于全面了解生物过程和发生机理也至关重要。

基因编码框(Open reading frame,ORF)的准确预测是一切基因组研究和相关研究和应用的基础工作。目前，基因编码框的预测主要通过对DNA序列特征进行判断，从而确定蛋白编码基因的起始和结束位置，进而推测基编码的蛋白质序列。现有数据表明，这种传统的预测方法对长编码框的预测具有较高的准确度，但是对小编码框(small ORF)的预测却几乎无能为力。传统的方法通过实验方法对小编码框逐个进行确认和验证，这种手段效率耗时耗力，在绝大多数生物中都不具有可操作性。目前，仅有酵母基因组中完成了约300个小编码框的实验验证工作。近年来，核糖体印迹测序(Ribo-seq)技术的兴起，快速、准确地对全基因组中的小编码框进行预测成为了可能。其基本原理是，被翻译的RNA序列会受到核糖体的保护，将这些受保护的序列提出来以后进行测序，就可以获得被翻译的序列，从而预测小编码框的位置。随着核糖体测序技术应用范围的不断扩展，许多基于核糖体测序数据预测小编码框的方法和软件也随之被开发出来。然而，由于目前这些主要方法都是在模式物种的研究中开发出来的，因此它们都基于一个理想的假设，即核糖体测序数据均具有较高的质量(完全呈3个碱基的周期性分布)。这一先决条件在模式物种中相对比较容易达到，但是在其它非模式物种中并不总是如此。甚至，即使在模式物种中，对不同细胞器中的核糖体保护序列进行测序也并不总是能获得满足条件的高质量数据。因此，对高质量核糖体印迹数据的要求极大地阻碍了这一技术在非模式物种中的应用，同时也限制了其应用范围。开发新的能够用于低质量核糖体测序数据分析的方法和软件对于推进这一技术的应用以及小编码框的研究具有重要的意义。

发明内容

本发明的目的在于克服上述现有技术之不足而提供一种从低质量核糖体印迹数据预测基因编码框的方法，本发明引入密码子使用频率，结合核糖体印迹数据的3碱基周期性，科学度量核糖体印迹数据的质量并合理分配相应权重，计算每个密码子位于核糖体P位点的概率，提取序列特征，通过统计分析，综合评定编码框的预测概率值，进而预测新的编码框。这一方法的发明有助于降低核糖体数据分析对数据质量的要求，快速扩展其应用范围。

为实现上述目的，本发明提供一种从低质量核糖体印迹数据预测基因编码框的方法，包括如下步骤：

S1，将原始测序的核糖体印记数据去掉接头后与基因组参考序列进行比对；

S2，采用multitaper算法分析不同长度的核糖体印迹序列(RPF)的3碱基周期性，保留评估合格的RPF，用于后续分析；

S3，通过基因组注释文件信息，提取转录本和已知编码框的序列和位置信息，同时获得全基因组所有转录本和已知编码框序列；

S4，对步骤S2中保留的RPF进行特征训练，并依此进行权重分配；

S5，计算各转录本上每一个碱基或每一个三碱基组合正好位于核糖体P位点(P-site)的概率；

S6，根据已知的各编码框的序列信息以及步骤S5中计算得出的P-site概率，同时提取基因编码框特征；

S7，根据S5中计算得出的每个碱基或三碱基组合正好位于核糖体P位点的概率，以及S6得到的基因编码框特征，预测出未知的基因编码框。

需要指出的是，S6中的基因编码框特征是指已知编码框的密码子使用频率。

优选的，S2中，各长度RPF的3碱基周期性通过multitaper算法进行评估，频率显示为0.33Hz～0.34Hz，P值≤0.01的RPF得以保留，用于后续分析。

更为优选的，S2中，各长度RPF的3碱基周期性通过multitaper算法进行评估，频率显示为0.33Hz或0.34Hz，P值≤0.01的RPF得以保留，用于后续分析。

优选的，S4包括：

S41，统计各个长度的RPF的5’端与P-site之间不同距离的出现频率；

S42，权重分配：根据S41中得到的各个RPF在相位0，1和2位置出现的频率，计算分布集中度。

更为优选的，S41具体是：通过分析包含已知编码框启始密码子或终止密码子的RPF与对应启始或终止密码子的位置信息，计算每条RPF 5’端与核糖体P位点(P-site)和\或核糖体A位点(A-site)的距离，统计各个长度的RPF的5’端与P-site之间不同距离的出现频率。

优选的，S42具体是：根据S41中得到的各个RPF在相位0，1和2位置出现的频率，计算分布集中度；分布集中度由复杂度Entropy描述，公式一如下：

其中，i表示不同的相位，i的取值范围为0，1和2，Pi为各个RPF在i相位上分布的比例；根据公式一计算出复杂度Entropy的值，分配RPF的权重为1–Entropy，相应的，序列特征的权重分配为Entropy。

优选的，S5具体是：根据Ribo-seq得到各RPF的位置信息以及各RPF的5’端与P-site之间的距离信息，计算各转录本上每一个碱基或者每一个三碱基组合正好位于P-site的概率。

优选的，S6，根据各编码框的序列信息以及S5中计算得出的P-site概率，提取编码框特征，具体包括如下步骤：

S61，Z-score：将S5计算得到的P-site的概率转化为Z-score；

S62，密码子使用频率：根据基因组中所有编码框的密码子使用情况，计算每个密码子的出现频率，然后计算每个已知编码框中密码子出现频率的平均值。

优选的，S7具体包括：

S71，根据S3中所有转录本的序列信息，对基因编码框候选序列进行提取和搜索；

S72，按照S6中的方法提取经S71得到的候选编码框的特征，进行多组统计检验，得到多个P值；

S73，P值合并：将S72中的多个P值经加权卡平方算法合并成最终P值；

S74，预测结果输出：控制S73中的P以及P编码框输出错误发现率FDR的值，将满足输出标准的候选编码框进行输出。

更为优选的，S7具体包括：

S71，依据S3中所有转录本的序列信息，提取所有候选编码框序列，依据标准为，拥有启始密码子(NUG)、终止密码子(UAG，UAA，UGA)并且其长度为3的整数倍数；优先搜索AUG起始的候选编码框，由长到短，逐一进行计算，AUG起始的候选编码框全部搜索完全且不满足输出条件后，再进行NUG编码框的搜索和计算；

S72，按照S6中的方法提取这些候选编码框的特征，进行四组统计检验，分别是：单尾检验(a):位于相位0上的Z-score值极显著大于位于相位1上的Z-score；单尾检验(b)：位于相位0上的Z-score值极显著大于位于相位2上的Z-score；单尾检验(c)：位于相位0上的密码子的使用频率值极显著大于位于相位1上的密码子频率；单尾检验(d)：位于相位0上的密码子的使用频率值极显著大于位于相位2上的密码子频率；

S73，P值合并：将S72中的多个P值经加权卡平方算法合并成最终P值：

S74，将预测的基因编码框RPF结果输出：输出P值≤0.001的修选编码框并根据Benjamini和Hochberg法控制编码框输出错误发现率FDR≤0.0001，满足这一标准的候选编码框进行最后的结果输出。

优选的，S7中，预测出未知的基因编码框RPF包括小编码框和\或正常的基因编码框。

为实现本发明的另一目的，本发明还提供一种预测基因编码框的系统，包括计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有预测基因编码框的计算机程序，所述预测基因编码框的计算机程序被至少一个处理组件执行时，能够实现上述的从低质量核糖体印迹数据预测基因编码框的方法的步骤。

本发明的有益效果是：

1.本发明引入密码子使用频率，结合核糖体印迹数据的3碱基周期性，科学度量核糖体印迹数据的质量并合理分配相应权重，计算每个密码子位于核糖体P位点的概率，提取序列特征，通过统计分析，综合评定编码框的预测概率值，进而预测新的编码框。这一方法的发明有助于降低核糖体数据分析对数据质量的要求，快速扩展其应用范围。提高对噪音数据的耐受程度，有效降低了对数据质量的要求。本发明的预测方法适用于：在模式生物中，某些细胞器难以获得高质量的核糖体印记数据，可以采用本发明的预测方法；在非模式生物中，如果较难以获得高质量的核糖体印记数据，可以采用本发明的预测方法预测基因编码框。本发明将预测基因编码框的范围大大增加了，对于推进小编码框的研究具有重要的意义。

2.为了方便的应用本发明的预测方法，将本发明的方法步骤以计算机程序的形式呈现给用户。用户将核糖体印记数据等必要信息作为输入，计算机程序可以输出预测得到的基因编码框。有利于提升用户的处理效率，在将本发明的预测基因编码框的方法推广向各个物种时，采用计算机程序的实现方式有助于提升预测编码框的效率，使本发明的预测方法得以更快速的普及。

附图说明

图1是本发明的技术路线示意图，即本发明的工作流程图；

图2是本发明候选编码框搜索策略的示意图；

图3是本发明应用实例，其中：图3(A)是实例数据RPF长度分布情况；图3(B)是三碱基周期性评估结果；图3(C)是RPF分布集中度计算及权重分配情况；图3(D)是预测效果评估结果；图3(E)是小编码框的预测结果；图3(F)是蛋白质质谱数据的支持证据；图3(G)是预测所得的ncsORF的进化分析，其中，图3G是一个热图，方块里的颜色深浅表示值的大小；

对图3进一步的放大形成如下的附图，用以更清晰的显示图3中各视图的细节：

图4是图3中A视图的放大图；

图5是图3中B视图的放大图；

图6是图3中C视图的放大图；

图7是图3中D视图的放大图；

图8是图3中E视图的放大图；

图9是图3中F视图的放大图；

图10是图3中G视图的放大图；

图11是本发明的从低质量核糖体印迹数据中预测基因编码框的方法的示意图。

本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

具体实施方式

实施例1

本发明公开了一种从低质量核糖体印迹数据中预测基因编码框的方法，该方法可以准确度量核糖体印迹的数据质量，并依此对数据进行初步过滤并合理分配相应权重，随后整合密码子使用频率辅助蛋白编码框的预测，该发明方法对核糖体印迹数据质量不敏感，具有较强的容错性。不仅如此，该发明方法在高质量的核糖体印迹数据中也具有极好的表现，能够全面准确预测翻译的编码框。因此，该方法适用于所有核糖体印迹数据。本发明的要点如下：

1.综合利用核糖体印迹和密码子使用频率进行蛋白编码框的预测。

2.利用multitaper算法和复杂度(熵)对核糖体印迹数据质量进行描述。

3.根据核糖体印迹数据的复杂度(熵)自动分配相应的权重，从而平衡数据质量的影响。

如上所述，本发明主要针对现在核糖体印迹测序数据分析方法中对数量质量要求过高的问题，提出一种新的预测基因编码框的方法，提高对噪音数据的耐受程度，有效降低了对数据质量的要求。需要注意的是：本发明仅适用于有参考基因组序列和注释信息的物种。

请参阅图1和图4，本发明方法主要包括以下几个步骤：

(1)基因组比对

将原始的核糖体印记测序数据去掉接头后与基因组参考序列进行比对。基因组参考序列可以从公开的渠道获取。

步骤(1)进行基因组比对的目的是：获取核糖体印迹序列在基因组上对应的位置信息。基因组参考序列就是已知的基因组序列，将核糖体印迹数据与之进行比对是为了获取他们在基因组上的位置信息。若比对结果不对，则后续所有预测都不对。这也是本发明预测方法的实现需要参考基因组序列的原因之一。

(2)核糖体印迹数据质量评估

通过分析核糖体印迹数据不同长度RPF的3碱基周期性，对完全不具有周期性的数据进行过滤。具体做法是：将各长度的3碱基周期性通过multitaper算法进行评估，频率显示为0.33Hz～0.34Hz，P值≤0.01的RPF得以保留，用于后续分析。

上述的步骤(2)，包含了数据过滤的操作，具体是：将完全不可用的数据过滤掉，保留评估合格的数据。采用multitaper算法进行数据质量评估，质量评估的目的是为数据过滤提供一个明确的过滤标准。

(3)转录本和已知编码框组装

通过基因组注释文件信息，提取转录本和已知编码框的序列和位置信息，获得全基因组所有转录本和已知编码框序列。

上述步骤(3)的目的或意义是：编码框是从转录本的序列基础上进行预测的。已知编码框的序列信息用于训练密码子使用频率，其位置信息用于训练RPF 5’端与对应P-site的距离信息。

(4)核糖体印迹数据(RPF)特征训练和权重分配

①特征训练：通过提取比对到已知编码框启始或终止密码子的RPF比对信息，计算每条RPF 5’端与核糖体P位点(P-site)和\或核糖体A位点(A-site)的距离，统计各个长度的RPF的5’端与P-site之间不同距离的出现频率。

对上述步骤(4)①进行优化：选择计算每条RPF 5’端与核糖体P位点(P-site)的距离这一方案即可。这是因为：A和P相隔3个碱基，是确定的信息。

步骤(4)中的特征训练的目的是：获得每一条RPF的5’端距离其所对应的P位点的距离信息。

步骤(4)中的特征训练的意义或作用是：训练RPF 5’端与对应P-site的距离信息。该信息将用于确定每一条RPF所对应的P-site位置。需要注意：并不是每一条RPF都明确知道其对应的P-site，只有包含已知起始或终止密码子的RPF才可以获取这一信息；通过这一部分RPF训练获得该距离信息后再用于其他RPF。

②权重分配：根据各RPF在相位0，1和2位置出现的频率计算其分布集中度。此处的分布集中度是指相位分布的集中度。分布集中度由复杂度(熵)进行描述，公式如下：

其中，i表示不同的相位(0，1和2)，P_i为RPF在相位i上分布的比例。根据计算所得复杂度Entropy的值对RPF分配相应的权重为(1–Entropy)，相应的，序列特征的权重分配为Entropy。

步骤(4)中，“对RPF分配相应的权重”，该权重为一个系数，用于确定在后续预测过程中该证据的贡献度。具体而言：RPF质量越高，获得的权重越高，对后续预测的贡献越大；相反，RPF质量越低(噪音高)，其对预测的贡献越小，这时预测结果更依赖于其他证据的支持，从而降低RPF噪音对预测结果的不利影响。“序列特征”，指序列自身的特征，相对于RPF而言，RPF是非序列特征。这里具体指密码子的使用频率。

(5)计算P-site概率

根据核糖体印记测序(Ribo-seq)得到各RPF的位置信息以及其5’端与P-site之间的距离信息，还需指出：5’端与P-site之间的距离信息不是一个确定值，而是一系列值，我们这里采用3个值，每个值会对应一个概率。计算方法见第(4)步的特征训练部分：通过提取比对到已知编码框启始或终止密码子的RPF比对信息，计算每条RPF 5’端与核糖体P位点(P-site)或核糖体A位点(A-site)的距离；计算各转录本上每一个碱基或三碱基组合正好位于P-site的概率，并将其转化为Z-score，即进行数据的标准化。若采用计算各转录本上每一个碱基正好位于P位点的概率的方案，则：每一个碱基都会获得一个概率值，其代表的是以这个碱基为起始的“三碱基组合”位于P位点的概率值。

需要指出的是：

a)步骤(5)中的位置信息是指RPF 5’端的位置，是通过与基因组比对所获得的。

b)步骤(5)中的三碱基组合进一步的限定是：连续排布的三个碱基所构成的组合。

c)如果需要采用计算各转录本上每一个三碱基组合正好位于P位点的概率的方案，则应该将该方案应该理解为：若连续的三个碱基组合在当前所检测的物种所适用的遗传密码规则下对应某一密码子，则计算该密码子正好位于P位点的概率，按照上述方式，将当前转录本上所有可能的密码子组合都计算出P位点的概率，更进一步的，按照上述的方式，将各转录本都完成计算。

(6)已知编码框特征提取

根据各编码框的序列信息以及上一步中计算得出的P-site概率，提取编码框特征，如下：

①Z-score：计算各密码子正好位于P-site的概率，并转化为Z-score。

②密码子使用频率：根据基因组中所有编码框的密码子使用情况，计算每个密码子的出现频率，然后计算每个已知编码框中密码子频率的平均值。

需要指出的是：第(4)步训练的是RPF的特征，RPF里面包含有实际测量到的编码框信息。步骤(6)训练的是已知编码框的序列特征。步骤(4)的特征训练结果和步骤(6)的特征提取结果会共同用于未知编码框的预测。

(7)编码框的预测

①编码框候选序列提取和搜索(请参阅图2)：依据(3)中所有转录本的序列信息，提取所有候选编码框序列，依据标准为，拥有启始密码子(NUG)、终止密码子(UAG，UAA，UGA)并且其长度为3的倍数。优先搜索AUG起始的候选编码框，由长到短，逐一进行计算，AUG起始的候选编码框全部搜索完全且不满足输出条件后，再进行NUG编码框的搜索和计算。

②统计检验：按照(6)中的方法提取这些候选编码框的特征，进行四组统计检验，分别是(a)位于相位0上的Z-score值极显著大于(单尾检验)位于相位1上的Z-score；(b)位于相位0上的Z-score值极显著大于(单尾检验)位于相位2上的Z-score；(c)位于相位0上的密码子的使用频率值极显著大于(单尾检验)位于相位1上的密码子频率；(d)位于相位0上的密码子的使用频率值极显著大于(单尾检验)位于相位2上的密码子频率。

③P值合并：以上统计所得的4个P值(P value，P值是用来判定假设检验结果的一个参数)经加权卡平方算法(Weighted chi-square method)合并成最终P值，计算方法如下，

首先按照步骤(4)中分配的权重，将P值转化为卡平方值，公式如下：

其中M表示合并后的卡方值，i为第i个检验，Pi为第i个检验的P值，wi为第i个P值的权重，由于wi之和须为1，且RPF和密码子使用频率各进行了两次检验，因此，相对应P值的权重为上一步中计算所得RPF/密码频率权重的一半。

计算自由度(k)

k＝2{E(M)}²/var(M)

其中，

s_i为P_i单独转化后的卡方值，s_i＝-2×w_i×ln(P_i)

其中，w_i,w_j为相的权重，与以上公式等价。ρ_ij为第i个检验与第j个检验之间的相关性。ρ又可以从计算所得的P值间接估算得出。如下，

其中，

为s_i的平均值，由于q_t的期望值E(q_t)＝4–(0.75ρ²+3.25ρ)，所以计算可得

0.75ρ²+3.25ρ+E(q_t)–4＝0

最后可以求解ρ的近似值为-2.167+(10.028-4q_t/3)^0.5。

根据计算所得的自由度k和合并后的卡方值，根据卡方分布2χ² _k/k获取对应的P值。

④编码框输出错误发现率(FDR)控制

输出P值≤0.001的修选编码框并根据Benjamini和Hochberg法控制FDR≤0.0001，满足这一标准的候选编码框进行最后的结果输出。

实施例1主要涉及一种利用低质量的核糖体印迹数据进行预测蛋白编码框的方法。蛋白编码框(包括小编码框)的准确预测是所有基因相关研究和应用的基础。核糖体印迹测序技术的兴起使得能够更加准确的对蛋白编码框进行预测，特别是使小编码框的预测成为可能。虽然已有许多软件和流程可以用来从核糖体印迹数据中预测蛋白编码框，但是这些工具的使用都必须基于一个理想的条件，即核糖体印迹数据均具有较高的质量(完全呈3碱基的周期性分布)。这一条件的满足需要极高的实验技术和昂贵的试剂和设备，极大的制约了该技术的应用拓展。此外，高质量的核糖体印迹数据通常长度较短(28nt)，在基因组上会有多个比对位点，会引入大量的错误，不利于后续研究的开展。总的来说，目前已有的流程和工具对低质量的核糖体印迹数据完全无能为力。为了解决低质量核糖体印迹数据无法使用，而高质量核糖体印迹数据又容易引入错误的问题，本发明提取密码子使用频率，结合核糖体印迹数据的3碱基周期性，科学度量核糖体印迹的数据质量并合理分配相应权重，计算每个密码子位于核糖体P位点的概率，提取序列特征，通过统计分析综合评定编码框的预测概率值，进而预测新的编码框。本发明将大幅降低相关工作对核糖体印迹数据质量的要求，将极大的促进核核体印迹技术应用的拓展，特别是在农作物研究中的应用。

对于上一段的论述，还需要进一步的指出：

a)权重分配的多寡取决于数据质量，核糖体印迹数据质量越高，其分配的权重就越高。

b)本发明的预测方法并不只局限于在“农作物研究中应用”，在动物、植物、微生物领域都可以使用本发明的预测方法，且都表现很好。相对而言，动物，微生物和人里面通常数据质量比较高，现有的方法可以较好处理。核糖体印记数据质量低的情况通常会在植物物种出现，特别是在非模式物种中遇到。也就是说，本发明的基因编码框预测方法还可以处理现有的预测方法不能处理的低质量核糖体印记数据。

实施例2：拟南芥膜结合核糖体数据的分析

(1)实验数据从NCBI下载(GEO编号：GSE82041)，该数据由LiShengben等于2016年发表于elife，文章名称为“Biogenesis of phased siRNA on membrane-bound polysomesin Arabidopsis”。实验中通过分离结合于膜上的核糖体，并对其保护的mRNA片段进行测序，得到MBP(membrane-bound polysomes)Ribo-seq数据。在制备MBP保护的片段的过程中，由于裸露RNA的降解通常不够完全，导致核糖体印记数据(Ribo-seq)质量较低，不能呈现出很好的3碱基周期性。

(2)请参阅图3至图10，利用本发明方法，首先对该数据进行质量评估。结果显示该数据中RPF长度分布不集中(图3(A)和图4)，理论上来说，真核生物中核糖体的印迹长度为28个核苷酸(nt，nucleotide)，因此RPF长度理应集中出现在28nt。图3A显示该组数据中，RPF长度的分布范围从18nt到35nt不等，分布范围较广，虽然在32nt处出现一个峰值，但是总体占比并不高，仅10％左右，并且这一值也远远偏离了理论值(28nt)，这表明该数据的产生过程中，祼露mRNA的降解并不完全，导致剩余的由核糖体保护的片段(RPF)长度不一，这将导致RPF的分辨率和精确度不足。这也表现在该数据的三碱基周期性也不强(图3(B)和图5)，理论上讲，由于密码子长度为3个碱基，各个核糖体印迹之间的距离应该为3的倍数，最小距离为3个碱基，其在转录本序列上的分布呈现出3碱基的周期性，在multitaper的检验结果中，将表现为频率峰值为1/3，且P值极显著，周期性越好，其P值越小，通常理想情况下，-log10(P-value)>＝10。图3B显示大部分RPF的频率峰值并不出现在1/3处，且P值较大，图中深色线表示长度为32nt(丰度最高)的RPF，数据显示其值约为3，刚刚通过multitaper检验(cutoff＝2)。RPF对应P-site的分布集中度不强(图3(C)和图6)，该图展示了长度为32nt的RPF P-site分布的集中度，计算所得其熵值为0.862。理想情况下，如果RPF仅对应唯一的P-site，熵值计算将为0，而如果RPF对应3个P-site，且分布平均，熵值计算将为1。图3C显示该组数据熵值为0.862，接近1而远离0，因此表示该组数据的分布不够集中。我们据此，为RPF和密码子使用频率分配相应的权重(RPF：0.138，密码子频率：0.862)，由于该数据中RPF在0，1，2位上的分布集中度不够高，我们更多的使用密码子频率(权重为0.862)进行编码框的预测。利用本发明方法，有76％的已知编码框得以成功预测，并且准确率高达98％，综合评分也高达86％[综合评分＝2×召回率×准确率/(召回率+准确率)](图3(D)和图7)，并且成功预测了1471个小编码框，其中含有114个uORF，93个ouORF，245个dORF，232个odORF，653个teORF，121个pORF，13个ncsORF(图3(E)和图8)。已发表的蛋白质质谱数据分析显示这些预测得到编码框均受到了很好的支持(图3(F)和图9)，图3(F)和图9中，横向虚线表示基因组中所有已知编码框受到蛋白质质谱数据的支持率，我们以此为参考进行比较，由此图可见，利用本方法从该数据中预测所得到的annotated ORF的质谱支持率明显高于整体水平(虚线表示)，其他几类(uORF,ouORF,dORF,odORF,teORF,pORF和ncsORF)为小编码框，由于长度较短，能产生的肽段较少，因此不容易被检测到，所以支持率相对较低，特别是ncsORF，由于其数量少，因此在质谱数据中没有检出，这些都是正常现象。为了进一步验证ncsORF的准确性，我们对所预测得到的ncsORF序列进行进化分析，通过其序列保守性来证实该预测的准确性。图3(G)和图10显示，大部分预测的ncsORF都表现出了较强的保守性，具体来说，有5个ncsORF从苔藓中就开始出现，其序列在所有植物分支中都非常保守，另一部分(4个)ncsORF从十字花科植物开始出现，并且在这一分支中非常保守，据此我们可以推断这些ncsORF是具有重要生物学功能的，这些预测结果是正确的。。

实施例2是实施例1的一个具体实例。

实施例3

本发明还公开一种预测基因编码框的系统，包括计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有预测基因编码框的计算机程序，所述预测基因编码框的计算机程序被至少一个处理组件执行时，能够实现上述的从低质量核糖体印迹数据预测基因编码框的方法的步骤。

实施例3主要是解决问题是：现有的预测基因编码框的系统只能处理高质量核糖体印记数据，对于低质量核糖体印记数据无能为力。

所述存储介质存储器可以是ROM或可存储静态信息和指令的其他类型的静态存储设备，RAM或者可存储信息和指令的其他类型的动态存储设备，也可以是EEPROM、CD-ROM或其他光盘存储、光碟存储(包括压缩光碟、激光碟、光碟、数字通用光碟、蓝光光碟等)、磁盘存储介质(包括机械硬盘、固态硬盘、混合硬盘等)或者其他磁存储设备(包括磁带)、或者能够用于携带或存储具有指令或数据结构形式的期望的程序代码并能够由计算机存取的任何其他介质(包括SD卡等)，但不限于此。所述存储介质可以是存放于本地，也可以是设置在云端。

所述处理组件是处理器，处理器可以是CPU，通用处理器，DSP，ASIC，FPGA或者其他可编程逻辑器件、晶体管逻辑器件、硬件部件或者其任意组合。处理器也可以是实现计算功能的组合，例如包含一个或多个微处理器组合，DSP和微处理器的组合等。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (10)

1.一种从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，包括如下步骤：

S1，将原始测序的核糖体印记数据去掉接头后与基因组参考序列进行比对；

S2，采用multitaper算法分析不同长度的核糖体印迹序列(RPF)的3碱基周期性，保留评估合格的RPF，用于后续分析；

S3，通过基因组注释文件信息，提取转录本和已知编码框的序列和位置信息，同时获得全基因组所有转录本和已知编码框序列；

S4，对步骤S2中保留的RPF进行特征训练，并依此进行权重分配；

S5，计算各转录本上每一个碱基或每一个三碱基组合正好位于核糖体P位点(P-site)的概率；

S6，根据已知的各编码框的序列信息以及步骤S5中计算得出的P-site概率，同时提取基因编码框特征；

S7，根据S5中计算得出的每个碱基或三碱基组合正好位于核糖体P位点的概率，以及S6得到的基因编码框特征，预测出未知的基因编码框。

2.根据权利要求1所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S2中，各长度RPF的3碱基周期性通过multitaper算法进行评估，频率显示为0.33Hz～0.34Hz，P值≤0.01的RPF得以保留，用于后续分析。

3.根据权利要求1所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S4包括：

S41，统计各个长度的RPF的5’端与P-site之间不同距离的出现频率；

S42，权重分配：根据S41中得到的各个RPF在相位0，1和2位置出现的频率，计算分布集中度。

4.根据权利要求3所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S41具体是：通过分析包含已知编码框启始密码子或终止密码子的RPF与对应启始或终止密码子的位置信息，计算每条RPF 5’端与核糖体P位点(P-site)和\或核糖体A位点(A-site)的距离，统计各个长度的RPF的5’端与P-site之间不同距离的出现频率。

5.根据权利要求3所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S42具体是：根据S41中得到的各个RPF在相位0，1和2位置出现的频率，计算分布集中度；分布集中度由复杂度Entropy描述，公式一如下：

其中，i表示不同的相位，i的取值范围为0，1和2，Pi为各个RPF在i相位上分布的比例；根据公式一计算出复杂度Entropy的值，分配RPF的权重为1–Entropy，相应的，序列特征的权重分配为Entropy。

6.根据权利要求1所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S5具体是：根据核糖体印记测序Ribo-seq得到各RPF的位置信息以及各RPF的5’端与P-site之间的距离信息，计算各转录本上每一个碱基或者每一个三碱基组合正好位于P-site的概率。

7.根据权利要求1所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S6，根据各编码框的序列信息以及S5中计算得出的P-site概率，提取编码框特征，具体包括如下步骤：

S61，Z-score：将S5计算得到的P-site的概率转化为Z-score；

S62，密码子使用频率：根据基因组中所有编码框的密码子使用情况，计算每个密码子的出现频率，然后计算每个已知编码框中密码子出现频率的平均值。

8.根据权利要求1所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S7具体包括：

S71，根据S3中所有转录本的序列信息，对基因编码框候选序列进行提取和搜索；

S72，按照S6中的方法提取经S71得到的候选编码框的特征，进行多组统计检验，得到多个P值；

S73，P值合并：将S72中的多个P值经加权卡平方算法合并成最终P值；

S74，预测结果输出：控制S73中的P以及P编码框输出错误发现率FDR的值，将满足输出标准的候选编码框进行输出。

9.根据权利要求8所述的从低质量核糖体印迹数据预测基因编码框的方法，其特征在于，S7具体包括如下步骤：

S71，依据S3中所有转录本的序列信息，提取所有候选编码框序列，依据标准为，拥有启始密码子NUG、终止密码子包括UAG，UAA和UGA并且候选编码框序列长度为3的整数倍数；优先搜索AUG起始的候选编码框，由长到短，逐一进行计算，AUG起始的候选编码框全部搜索完全且不满足输出条件后，再进行NUG编码框的搜索和计算；

S72，按照S6中的方法提取这些候选编码框的特征，进行四组统计检验，分别是：

单尾检验(a):位于相位0上的Z-score值极显著大于位于相位1上的Z-score；

单尾检验(b)：位于相位0上的Z-score值极显著大于位于相位2上的Z-score；

单尾检验(c)：位于相位0上的密码子的使用频率值极显著大于位于相位1上的密码子频率；

单尾检验(d)：位于相位0上的密码子的使用频率值极显著大于位于相位2上的密码子频率；

S73，P值合并：将S72中的多个P值经加权卡平方算法合并成最终P值：

S74，将预测的基因编码框RPF结果输出：输出P值≤0.001的修选编码框并根据Benjamini和Hochberg法控制编码框输出错误发现率FDR≤0.0001，满足这一标准的候选编码框进行最后的结果输出。

10.一种预测基因编码框的系统，包括计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有预测基因编码框的计算机程序，所述预测基因编码框的计算机程序被至少一个处理组件执行时，能够实现如权利要求1至9任一项所述的从低质量核糖体印迹数据预测基因编码框的方法的步骤。

Images