基于k茎的核糖核酸假结结构的预测方法及装置

摘要

本发明提供一种基于k茎的核糖核酸(RNA)假结结构的预测方法及装置，预测方法包括以下步骤：输入一段核糖核酸碱基序列；定义假结、k(k≥1)茎；从左向右查找RNA碱基和k茎，对查找出的所有k茎进行确定标记；根据k茎的交叉形成假结的特性，查找假结；计算出包含k茎的核糖核酸假结结构的最小自由能量；输出核糖核酸的假结结构。本发明所涉及的方法的搜索速度快、正确率高，敏感性和特异性等方面都优于其他相关算法，如PKNOTS算法等。本方法在平面假结的预测上比PKNOTS算法更有效。

说明书

技术领域

本发明属于生物信息工程领域，涉及一种对核糖核酸(以下简称 为RNA)的假结结构进行预测的方法，尤其涉及基于k茎的RNA假 结结构预测的方法及装置。

背景技术

RNA是生物系统内最为重要的大分子之一，它在生物体内行使 多种功能，是合成蛋白质的模板。RNA二级结构预测用于蛋白质功 能分析，是RNA三级结构预测的基础。假结(pseudoknot)是RNA中 最广泛的结构单元，是非常复杂和稳定的RNA结构，假结在RNA 分子中具有构造、催化和调节功能，是目前RNA结构预测研究的关 键点。

RNA二级结构预测采用的方法主要有两种：早期采用的是序列 对比分析方法，即对于在不同有机体中起相同生物功能的一级结构进 行比较，此方法的困难之处在于：许多RNA分子的同源序列不易得 到；需要大量人力，效率较低，所以目前主要采用的是最小自由能量 方法。

最小自由能量算法的理论依据是稳定的RNA二级结构的自由能 量最小。基于最小自由能量算法的PKNOTS算法使用O(n⁶)时间和 O(n⁴)空间计算任意的平面假结和部分非平面假结。PKNOTS算法仅 能计算长度短于140个碱基的RNA序列，不能满足较长RNA序列 结构预测的需要。PknotsRG算法计算由两个茎区构成的简单的嵌套 假结，其中任意两个假结为并列或嵌套关系。事实上，由内环和凸起 构成的假结在RNA中普遍存在，交叉假结也具有重要作用。因此， 两者都不能被忽略。平面假结是最广泛的假结子类，包含上面提到的 由内环和凸起构成的假结以及交叉假结的情况。PseudoBase数据库的 所有序列中仅一个序列折叠为一个非平面假结，其余序列都折叠为平 面假结。因此我们主要考虑任意平面假结的计算。

Zuker首次将动态规划算法用于最邻近邻居模型，提出了MFOLD 算法，经过二十多年的不断改进和发展，现己成为国际上广泛使用的 RNA二级结构预测方法，对于包含n个核苷酸的RNA序列，MFOLD 算法使用O(n³)时间和O(n²)空间预测其最优二级结构，目前对于长度 小于700个核普酸的RNA序列，MFOLD算法可正确预测大约73％ 的RNA基对，但对于长序列和部分子类的预测正确率会降低，该算 法仅仅给出了三级结构预测的粗略框架，另外由于算法本身的限制， MFOLD算法不能预测假结和更复杂的三级结构。

发明内容

本发明解决的技术问题是使得对RNA结构预测、尤其是对基于 k茎包含假结的RNA结构进行预测方法，降低预测的时间复杂性和 空间复杂性，提高预测准确性。

本发明涉及的一种基于k茎的核糖核酸假结结构的预测方法包括 以下步骤：

输入一段核糖核酸碱基序列；

定义假结、k茎，k≥1；

从左向右查找碱基和k茎，对查找出的所有k茎进行标记；

根据两个以上k茎碱基对的交叉构成假结结构特性，查找假结；

计算出包含k茎的核糖核酸假结结构的最小能量；

输出核糖核酸的假结结构。

1茎(记为S₁[i,j])由碱基对(i,j)和(r,s)∈S所封闭，设(k-1)茎由 碱基对(r’,s’)和(k,l)∈S所封闭，i<r<r’<k<l<s’<s<j，v＝r’–r+s-s’>2， 则由(i,j)和(k,l)∈S所封闭的结构称为k茎(记为S_k[i,j])。其中，两 个k茎中碱基对的交叉构成假结。从左向右查找碱基时，首先查找1 茎，若找到1茎，则对1茎中的所有碱基标记，同理，查找2茎、3 茎……k茎，若找到，则对k茎中的所有碱基标记。

一种基于k茎的核糖核酸假结结构的预测装置包括：

输入单元：其输入一段核糖核酸碱基序列；

定义单元：定义1茎、2茎……k茎；

查找单元：从左向右查找碱基，对查找出的所有1茎、2茎…k茎 中的碱基进行标记；

假结结构查找单元：根据两个以上k茎碱基对的交叉构成假结结 构特性，查找假结；

假结计算单元：计算出包含k茎的核糖核酸假结结构的最小能量；

输出单元：其根据最小能量原理，输出核糖核酸碱基序列的假结 结构。

本发明的方法的搜索速度、正确率、敏感性和特异性都优于 PKNOTS算法。因此本方法在平面假结的预测上比PKNOTS算法更有 效。

附图说明

图1是本发明的基于k茎的RNA假结结构的预测方法流程图；

图2是本发明的k茎处理的流程图；

图3是对应图1中用于预测RNA假结结构的预测装置；

图4是本发明的一个RNA假结结构的例子；

图5是本发明的计算RNA假结结构最小能量中W和V的表示图 示。

具体实施方式

首先说明关于RNA序列、碱基对、假结等的概念。

RNA一级结构：RNA分子侧链上四种碱基的排列顺序表示。一 般来说RNA碱基序列从5′开始到3′结束，这样整个序列s表示为 s＝s₁s₂…s_n，s_i表示RNA序列的第i个碱基，s_i∈{A,U,G,C}，RNA子 序列s_i,j是s的一个序列片段，表示为：s_i,j＝s_i…s_j。

碱基对：如果s_i·s_j∈{AU,CG,GU}，则s_i·s_j构成碱基对。碱基 对中堆叠的能量为负值。

RNA二级结构：RNA序列中的一组基对构成的集合，以S表 示。对于任意基对，如果s_i·s_j∈S、s_i′·s_j′∈S且若i＝i′，则j＝j′， 亦即，一个碱基不可同时与两个及两个以上的碱基构成基对。

假结：如果基对s_i·s_j与s_i′·s_j′∈S，如果i＜i′＜j＜j′，则 序列s_i...s_i′...s_j...s_j′构成假结结构。

图1是根据本发明的用于预测基于茎区RNA假结结构的预测方 法的流程图；本发明的方法包括以下步骤：输入一段核糖核酸序列； 定义假结、k茎(k≥1)；从左向右查找碱基，对查找出的所有k茎进 行标记；根据两个以上k茎碱基对的交叉构成假结结构特性，查找假 结；计算出包含k茎的核糖核酸假结结构的最小能量；输出核糖核酸 的假结结构。图3是对应图1中用于预测基于茎区的RNA假结结构的预 测装置。RNA假结结构的预测装置包括：输入单元：其输入一段核糖 核酸碱基序列；定义单元：其定义假结和定义k茎，k≥1；查找单元： 从左向右查找碱基，对查找出的所有1茎、2茎…k茎中的碱基进行标 记；假结结构查找单元：根据两个以上k茎碱基对的交叉构成假结结 构特性，查找假结；假结能量计算单元：计算出包含k茎的核糖核酸 假结结构的最小能量；输出单元：其根据最小能量原理，输出核糖核 酸碱基序列的假结结构。

图2是根据本发明的k茎处理的流程图：输入一段s＝s₁s₂…s_n序列， 从左向右查找碱基，如果存在i、j，使得s_i和s_j配对，j-i≥6，并且s中 存在三个以上连续的相邻基对s_i·s_j、s_(i+1)·s_(j-1)。。。、s_k·s_l，则基对s_i·s_j和s_k·s_l封闭的区间确定为1茎；对1茎中所有配对的碱基进行标记； 在1茎封闭的游离碱基中继续查找配对的碱基，如果存在三个以上基 对，确定为2茎；对2茎中所有配对的碱基进行标记；在1茎和2茎封闭 的游离碱基中继续查找配对的碱基，如果存在三个以上基对，确定为 3茎；对3茎中所有配对的碱基进行标记......直到查找到k茎。如果存在 两个以上k茎碱基对的交叉，则构成假结。

定义:RNA子序列S_i,j中，如果(i,j),(i+1,j-1),…，(k,l)都是基 对，i<k<l<j，则由(i,j)和(k,l)∈S所封闭的结构称为1茎，表示为S₁[i, j]。若1茎S₁[i,j]由(i,j)和(r,s)∈S所封闭，1茎S₁[r’,s’]由(r’,s’)和 (k,l)∈S所封闭，i<r<r’<k<l<s’<s<j，v＝r’–r+s-s’>2，则由(i,j)和(k, l)∈S所封闭的结构称为2茎,表示为S₂[i,j]。

同理，如果S₁[i,j]由(i,j)和(r,s)∈S所封闭，(k-1)茎由(r’,s’)和(k, l)∈S所封闭，i<r<r’<k<l<s’<s<j，v＝r’–r+s-s’>2，则由(i,j)和(k,l)∈S 所封闭的结构称为k茎,表示为S_k[i,j]，S_k[i,j]的最小能量表示为ES_k(i, j)，k茎S_k[i,j]的长度表示为LS_k(i,j)＝k-i+1或RS_k(i,j)＝j-l+1。

设2茎S₂[i,j]由两个嵌套的1茎和其内部未配对碱基构成。设 E₂(r,r’:s’,s)表示基对(r,s)和(r’,s’)构成的2环结构的能量，ES₁(i,j) 和ES₁(r’,s’)分别表示由基对(i,j)和(r’,s’)封闭的1茎的能量，则ES₂(i, j)＝ES₁(i,j)+E₂(r,r’:s’,s)+ES₁(r’,s’)。同理ES_k(i,j)＝ES₁(i,j)+E₂(r,r’: s’,s)+ES_k-1(r’,s’)。

设LS(i,j)∈{LS₁(i,j),LS₂(i,j)},ES(i,j)∈{ES₁(i,j),ES₂(i,j)}。在 本发明的方法中，1茎和2茎的自由能量和长度使用O(n³)的时间预 处理并分别存于三角矩阵ES(i,j)、LS(i,j)中，其计算过程见程序1。

同理，由ES_k(i,j)的计算公式可知，计算k茎的时间复杂度为 O(n³)，空间复杂度为O(n²)。k茎(k≥3)的计算由后面的动态规划算法 实现。

k茎由茎和2环构成，其自由能量为其所含的堆叠和环的能量之 和，任意假结可分解为k茎和多分枝环。

实施例1：

在RNA假结结构预测中,若k茎中k＝1或k＝2时，相关1茎和2 茎的程序计算如下所述。

程序1：1茎和2茎的能量和长度的计算

图4给出一个简单的假结。使用两个1茎(S₁[1,19]、S₁[7,30])和 三个子序列(s_6,6、s_13,14、s_20,24)构成一个假结。由于每个1茎由两个参 数确定，1茎的存储需要O(n²)空间，因此计算假结的时间复杂度为 O(n⁴)，空间复杂度为O(n²)。

由图4知：W(1,30)＝ES₁(1,19)+ES₁(7,30)+W(6,6)+W(13,14)+W (20,24)

实施例2：

给定一个序列s＝s₁s₂…s_n，序列片段s_i,j＝s_i…s_j，1＜i＜j＜n。设 W(i,j)是子序列S_i,j对应的包含假结的二级结构S的最小能量。设V(i,j) 是s_i和s_j构成基对(i,j)的情况下，子序列S_i,j对应的包含假结的二级 结构S的最小能量。

图5给出W(i,j)和V(i,j)的计算图式。包含假结结构的W(i,j)由 下列4种情况计算：

1)s_j是未配对碱基，碱基s_i和s_j-1配对关系未确定,如图5.1，计 算的W(i,j)＝W(i,j-1)；

2)s_i是未配对碱基，碱基s_i+1和s_j配对关系未确定，如图5.2， 计算的W(i,j)＝W(i+1,j)；

3)s_i和s_k，s_k+1和s_j不构成基对且在不同子序列S_i,k和S_k+1,j对 应的二级结构中，i<k<j，如图5.3，计算的$W(i,j)=$$\underset{i<k<j-1}{\min }(W(i,k)+W(k+1,j));$

4)s_i和s_j构成基对(i,j),如图5.4，计算的W(i,j)＝min(V(i,j))。

包含假结结构的V(i,j)由下列三种情况计算：

(1)S是由基对(i,j)封闭的1环，如图5.5，计算的 V(i,j)＝minE₁(i,j)；

(2)S是由基对(i,j)和(k,l)封闭的2环，如图5.6。

V(i,j)＝min(E₂(i,k:l,j)+V(k,l)),i<k<l<j,u＝(k-i+j-l)-2<U

(3)S是k环(k≥3)或者假结结构，i<k<j，如图5.6，计算的 $V(i,j)=\min \underset{i<h<j-1}{\min }(W_M(i+1,h)+W_M(h+1,j-1)+M+P),$其中M表示构成 一个多分枝环的惩罚值，P表示多分枝环中每一基对的惩罚值。W_M 和W的计算公式相同，但参数不同，W_M专门用于多分枝环内序列片 断的结构预测，而W仅用于无外部封闭基对时序列片断的结构预测。

本发明的方法与PKNOTS算法的实验比较

PknotsRG方法仅能预测由两个1茎构成的简单的嵌套假结，不 能预测交叉假结。使用C++实现本发明的方法，并与PKNOTS方法 进行比较。

表1本发明的方法与PKNOTS方法的计算时间比较

计算时间的比较见表1。本发明的方法在双核CPU:3.0GHz,内存 为4GB的PC机进行测试，而PKNOTS算法高性能计算机Silicon  Graphics Origin200进行测试。从表1可知，计算长度为75个碱基的 RNA序列，本发明的方法使用51秒，而PKNOTS算法使用20分钟。 计算长度为105个碱基的RNA序列，本发明的方法使用225秒，而 PKNOTS算法使用235分钟。计算长度为200个碱基的RNA序列， 本发明的方法使用72分钟，而PKNOTS算法不能计算。事实上，本 发明的方法可以成功预测长度超过1500个碱基的RNA序列的二级结 构。

虽然本发明所揭露的实施方式如上，但所述的内容只是为了便于 理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所 属技术领域内的技术人员，在不脱离本发明所揭露的精神和范围的前 提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明 的专利保护范围，仍须以所附的权利要求书所界定的范围为准。