CN109273047B - 一种基于模拟退火的核酸结构预测方法 - Google Patents

Abstract

本发明涉及一种基于模拟退火的核酸结构预测方法，通过最小茎区数和环中最小碱基数来确定随机点匹配列表和每组随机点对应的连续匹配数，生成茎区候选区域，再利用了模拟退火的思想，使得更快的产生有效匹配，最后引入能量评价函数，来提高RNA分子假结预测的准确率，从而达到降低了时间复杂度和空间复杂度，以及提高了RNA分子假结预测的准确率。

Description

一种基于模拟退火的核酸结构预测方法

技术领域

本发明属于生物信息工程领域，涉及一种RNA(核糖核酸)的二级结构预测的方法，尤其涉及基于模拟退火算法的RNA假结预测的方法。

背景技术

RNA序列的功能主要取决于它们的三维结构，假结又是RNA序列中最广泛的结构单元，是非常复杂和稳定的RNA结构，假结在RNA序列中具有构造、催化和调节功能。通过X射线衍射和核磁共振(NMR)的数据，可以推断出RNA长链功能的3D形式。这种方法虽然结果精确但是只有在拥有相关设备的环境下才能进行，所用设备也非常昂贵且非常耗时。因此，采用计算机和热力学模型预测RNA序列二级结构的方法被采用。研究人员也普遍认为寻找RNA的二级结构是理解生物功能的第一步，因此RNA假结预测是目前RNA结构预测研究的热点。

RNA序列设计问题是由药物设计应用领域激发的新兴研究课题。RNA设计的目的是根据自身的一级结构确定其二级结构，如图1，2所示。虽然使用暴力法可以从RNA一级结构中找到二级结构，但我们目前计算技术无法在合理的时间内得到一个理想结果。此外，RNA包含假结的二级结构预测问题已经被证明为是NP完全问题。

目前试验方法测定RNA二级结构成本高且非常耗时。比如：pknotsRG-mfe方法、NUPACK方法，Mfold方法等。其中Mflod方法不能预测假结结构，且它的时间复杂度为O(n³)。Rivas的pknotsRG-mfe方法方法处理了大量假结结构，提出一种基于最小自由能的动态规划算法模型，其时间复杂度为O(n⁴)，空间复杂度O(n²)，但是预测分子的最大长度不超过600。由Dirks和Pierce提出的NUPACK方法也是基于动态规划思想实现，该方法时间复杂度O(n⁵)，空间复杂度O(n⁴)。由此可见，这些预测方法存在时间复杂度和空间复杂度偏高、仅能预测严格符合严格限制的RNA二级结构且无法对长链的RNA序列有效。

由于对RNA二级结构预测方法的空间复杂度和时间复杂度直接影响预测成本，因此，如何使RNA二级结构、尤其是包含假结的RNA的二级结构的预测方法的时间复杂度和空间复杂度尽可能小，且确保预测结构的准确性，已经成为生物信息工程领域的一项重要的研究课题。

发明内容

本发明所要解决的技术问题在于提供一种使得对RNA假结结构进行预测方法的时间复杂度和空间复杂度尽可能小并提高预测结果的准确性的核酸二级结构预测算法。

本发明的技术方案可以通过以下技术措施来实现：

一种基于模拟退火的核酸结构预测方法，包括以下步骤：

S100，用户直接输入设置最小茎区数、环中最小碱基数、初始温度值T₀、临界温度T_临、最大假结交叉分组数，然后初始化分子新状态、分子旧状态；设置最小茎区数与环中最小碱基数目的在于，能够确定每组随机点的k连续匹配表；

S101，一般地，RNA序列S表示为x₁x₂x₃…x_n，其中x_i∈{A,C,G,U},序列的长度为n,1≤i≤j≤n；因此，将RNA序列的每个碱基，用该碱基所在的位置序号来替代，即表示为1，2，3,…,i,…,j,…,n，这种编码方式被称为长度编码；

S102，根据待预测的RNA序列，判断RNA序列中的碱基配对情况，当RNA序列中i，j位置发生碱基配对时，将长度编码中其对应的位置进行交换，即此时的RNA序列长度编码为：1，2，3,…,j,…,i,…,n；根据碱基配对情况计算符合碱基配对规则的随机点匹配列表：

(i，j，k)，

其中i，j分别表示RNA分子序列的第i个位置和第j个位置，k为随机点(i，j)的连续匹配数，即RNA序列长度编码为：1，2，3,…,j，j-1,…j-k+1,…,i+k-1,…,i+1,i,…,n；当k≥2时，该随机点组合k连续匹配可称为茎区。

S103，计算每组随机点(i，j)的连续匹配数列表：

该连续匹配数列表的计算方法为：根据随机点匹配列表中的任意两个数组(i，j，k)和(m，n，l)，若i＝m，j＝n，k≠l，那么随机点组合(i，j)对应的连续匹配数列表中有k、l；

S104，判断当前系统是否达到临界温度T_临，若没有达到，进入步骤S105，若达到进入步骤S113；

S105，判断当前环中最小碱基数是否达到上限，上限为分子长度n-2*最小茎区数，茎区由碱基对构成，因此，剩余游离的碱基(未配对的)个数为n-2*最小茎区数。若没有达到，则进入步骤S106，若达到，则进入步骤S112；

S106，根据随机点匹配列表，计算新的分子状态，具体过程为首先从随机点匹配列表中随机挑选一组随机点(i，j)，然后从该组随机点的k连续匹配数列表中随机挑选一个连续匹配数k，根据挑选的随机点组合和挑选的连续匹配数生成新的RNA分子长度编码，即分子新状态；

S107，根据步骤S106产生的分子新状态，计算出总的碱基配对数、平均每组碱基数和假结总交叉分组数；

根据下式计算新旧RNA分子状态的能量差：

ΔE＝E_新-E_旧 (1)

其中：E_新为分子新状态的能量，E_旧为分子旧状态的能量，

E＝

总的碱基匹配数×平均每组碱基数²×

[(总分组数-总假结交叉分组数)/总分组数] (2)或

E＝总的碱基匹配数×平均每组碱基数² (3)，

式(2)适用于碱基假结交叉分组数超过设定的最大假结交叉分组数的情况下，式(3)适用于碱基假结交叉分组数没有超过设定的最大假结交叉分组数的情况下；

其中，总分组数为随机点匹配列表中三元组(i，j，k)的个数，如RNA分子已匹配的三元组有(I，j，k)，(m，n，l)，(r，s，t)，则该总分组数大小为3；假结总交叉分组数为上述三元组两两交叉的个数；最大假结交叉分组数，为三元组两两交叉的个数的上限，指上述3个三元组两两交叉的个数最大值，一般设置不超过1。

S108，判断S107所得ΔE是否大于0，若大于0，则进入步骤S110，若小于0，则进入步骤S109；

S109，根据公式(4)进行计算：

其中，Random(0,1)为用随机函数Random生成一个0到1之间的数，T为系统当前温度，若等式成立，则进入步骤S110，否则进入步骤S111；

S110，将RNA分子状态替换为新分子状态，更新RNA分子长度编码；

S111，环中最小碱基数自增，再次进入步骤S105；

S112，系统温度T以T＝T₀*0.98的方式进行一次降温，然后，再次进入S104步骤中；

S113，输出RNA分子最终状态。

优选地，步骤S102中随机点匹配列表(i，j，k)中随机数i，j需要同时满足以下关系：

i<j

j-i-k>3

最小茎区数≤k≤2/3*序列长度n

若随机数不满足上述关系，则重新生成随机数，若满足，则判断是否满足k连续匹配，若不满足k连续匹配，则重新生成随机点，满足，则添加到随机点匹配列表中。

优选地，步骤S102生成随机点匹配列表(i，j，k)后对其进行k连续匹配验证，步骤如下：

首先对单独位置上的碱基组合依照Watson-Crick碱基配对规则进行校验，先把分子序列重新编码，编码规则按照A，C，G，U依次对应0，1，2，3，根据碱基配对规则，若第i位置与第j位置为基本匹配，即A-U，G-C或U-A，C-G配对时，需要满足如下条件：

RnaSeq[i]+RnaSeq[j]＝3，

其中，

RnaSeq[i]为i位置的碱基所代表的转录数字，

RnaSeq[j]为j位置的碱基所代表的转录数字，

当匹配为G-U或U-G配对时，则需满足以下条件：

RnaSeq[i]+RnaSeq[j]＝5

当从第i位置开始到第i+k-1位置为止分别与第j位置开始到第j-k+1位置为止，均满足上述条件(如图4所示)，则随机生成的三元组(i，j，k)满足k连续匹配。

附图说明

利用附图对本发明作进一步说明，但附图中的实施例不构成对本发明的任何限制。

图1是RNA分子一级结构平面图；

图2是RNA分子二级结构平面图；

图3是本发明的用于预测RNA的假结预测方法的框架流程图；

图4是分子序列第i,j位置发生k连续匹配的示意图一；

图5是分子序列第i,j位置发生k连续匹配的示意图二；

图6是本发明步骤S106中新状态产生流程图；

图7是本发明步骤S106中无效匹配的判断流程图；

图8是本发明与采用pknotsRE算法的技术方案的Sensitivity和Specificity对照表。

具体实施方式

为使本发明更加容易理解，下面将进一步阐述本发明的具体实施方式。

如图3所示，一种基于模拟退火的核酸结构预测方法，包括以下步骤：

S100，设置最小茎区数、环中最小碱基数、初始温度值T₀、临界温度T_临、分子新状态、分子旧状态、最大假结交叉分组数；最小茎区数设置为2，由于RNA序列不能剧烈折叠，环中都需要至少间隔三个碱基，因此环中最小碱基数需默认设置为3，初始温度值T₀可设为1500。

S101，将RNA序列S表示为x₁x₂x₃…x_n，其中x_i∈{A,C,G,U},序列的长度为n,1≤i≤j≤n；将RNA序列用长度编码表示为1，2，3,…,i,…,j,…,n，RNA的分子状态由长度编码进行描述。

S102，根据待预测的RNA序列，判断RNA序列中的碱基配对情况，当RNA序列中i，j位置发生碱基配对时，将其长度编码进行交换，即此时的RNA序列长度编码为：1，2，3,…,j,…,i,…,n；根据碱基配对情况计算符合碱基配对规则的随机点匹配列表：

(i，j，k)，

其中i，j分别表示RNA分子序列的第i个位置和第j个位置，k为随机点(i，j)的连续匹配数，即RNA序列长度编码为：1，2，3,…,j，j-1,…j-k+1,…,i+k-1,…,i+1,i,…,n，如图4，5所示。

当k≥2时，该随机点组合k连续匹配可称为茎区，即序列两个以上的相邻碱基对x_i·x_j,x_i+1·x_j-1,x_i+2·x_j-2,x_i+3·x_j-3…∈S所构成封闭的结构，i表示茎区的起始位置，j表示与i匹配的碱基位置，k表示茎区长度，也就是相邻碱基对的个数。

生成随机点匹配列表(i，j，k)后对其进行k连续匹配验证，步骤如下：

首先对单独位置上的碱基组合依照Watson-Crick碱基配对规则进行校验，先把分子序列重新编码，编码规则按照A，C，G，U依次对应0，1，2，3，根据碱基配对规则，若第i位置与第j位置为基本匹配，即A-U，G-C或U-A，C-G配对时，需要满足如下条件：

RnaSeq[i]+RnaSeq[j]＝3

当匹配为G-U或U-G配对时，则需满足以下条件：

RnaSeq[i]+RnaSeq[j]＝5

当从第i位置开始到第i+k-1位置为止分别与第j位置开始到第j-k+1位置为止，均满足上述条件(如图4所示)，则随机生成的三元组(i，j，k)满足k连续匹配。

此外，随机点匹配列表(i，j，k)中随机数i，j需要同时满足以下关系：

i<j

j-i-k>3

最小茎区数≤k≤2/3*序列长度n

若随机数不满足上述关系，则重新生成随机数，若满足，则判断是否满足k连续匹配，若不满足k连续匹配，则重新生成随机点，满足，则添加到随机点匹配列表中。

S103，计算每组随机点(i，j)的连续匹配数列表，将步骤S102中计算的随机点匹配列表进行合并、删除操作。根据随机点匹配列表中的任意两个数组(i，j，k)和(m，n，l)，若i＝m，j＝n，k≠l，那么随机点组合(i，j)对应的连续匹配数列表中有k、l；若k＝l，那么直接将(m，n，l)删除即可，保留一个(i，j，k)三元组。

S104，判断当前系统是否达到临界温度T_临，若没有达到，进入步骤S105，若达到进入步骤S113；

S105，判断当前环中最小碱基数是否达到上限，上限为分子长度n-2*最小茎区数，若没有达到，则进入步骤S106，若达到，则进入步骤S112。所述的最小碱基数为形成环时未发生配对的碱基数量，如图4中，i+k～i+k+l的碱基数量，累计L+1个,后面再判断最小碱基数是否达到rna.length/2这个上限，最小碱基数的循环迭代目的是为了找出当前序列的所有的随机点满足条件的k连续匹配；

S106，根据随机点匹配列表，计算新的分子状态，即利用当前状态产生新状态。如图6所示。首先，复制当前状态，从随机点匹配列表中随机挑选一组随机点(i，j)，然后从该组随机点的k连续匹配数列表中随机挑选一个连续匹配数k，根据挑选的随机点组合和挑选的连续匹配数生成新的RNA分子长度编码，即分子新状态；

由于在之前步骤中，已经验证了三元组(i，j，k)满足k连续匹配条件，因此直接根据新生成三元组进行RNA分子长度编码。在对新的RNA分子进行编码时，要检验新状态的茎区列表的茎区与旧状态的茎区列表的茎区是否有冲突。若有冲突，则直接去掉旧冲突匹配，即把冲突位置的RNA的长度编码还原，形成单独碱基配对。然后，按照RNA序列新状态的茎区列表的茎区进行长度编码。最后再根据长度编码情况，计算出每组随机点组合的连续匹配数。若该组随机点的连续匹配数小于最小茎区数，则被判断为无效匹配，如图7所示。然后剔除无效匹配，即把该组随机点位置的长度编码还原即可。

S107，根据步骤S106产生的分子新状态，计算出总的碱基配对数、平均每组碱基数和假结总交叉分组数；

根据下式计算新旧RNA分子状态的能量差：

ΔE＝E_新-E_旧 (1)

其中：

E＝

总的碱基匹配数×平均每组碱基数²×

[(总分组数-总假结交叉分组数)/总分组数] (2)或

E＝总的碱基匹配数×平均每组碱基数² (3)，

式(2)适用于碱基假结交叉分组数超过设定的最大假结交叉分组数的情况下，式(3)适用于碱基假结交叉分组数没有超过设定的最大假结交叉分组数的情况下；

S108，判断S107所得ΔE是否大于0，若大于0，则进入步骤S110，若小于0，则进入步骤S109；

S109，根据公式(4)进行计算：

若等式成立，则进入步骤S110，否则进入步骤S111。其中，Random(0,1)为用随机函数Random生成一个0到1之间的数，公式(4)表示接受RNA分子新状态的概率，概率会随着系统温度减小而减小，在系统初期，等式成立的概率偏大，而在系统温度快达到临界值时则等式成立的概率偏小。

S110，将RNA分子状态替换为新分子状态，更新RNA分子长度编码；

S111，环中最小碱基数自增，再次进入步骤S105；

S112，系统温度T以T＝T₀*0.98的方式进行一次降温，然后，再次进入S104步骤中；

S113，输出RNA分子最终状态，即当前分子的随机点匹配列表(茎区列表)。

如上所述，即使RNA序列中包含了非嵌套结构和假结结构，本发明的能够较准确地确定RNA的二级结构。本发明的技术方案通过最小茎区数和环中最小碱基数来确定随机点匹配列表和每组随机点对应的连续匹配数，生成茎区候选区域。再利用了模拟退火的思想，使得更快的产生有效匹配。最后引入能量评价函数，来提高RNA分子假结预测的准确率，从而达到降低了时间复杂度和空间复杂度，以及提高了RNA分子假结预测的准确率。

图8示出本发明与采用pknotsRE算法的技术方案的Sensitivity和Specificity对照表。在图8中，Sensitivity＝TP/RP,Specificity＝TP/(TP+FP)，其中TP表示RNA结构中正确预测的碱基对数量，FP表示RNA结构中错误预测的碱基对数量。RP表示RNA结构中真实的碱基对的数量。

最后所应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (2)

1.一种基于模拟退火的核酸结构预测方法，其特征在于，包括以下步骤：

S100，通过用户直接输入设置最小茎区数、环中最小碱基数、初始温度值T₀、临界温度T_临、最大假结交叉分组数，然后初始化分子新状态、分子旧状态；

S101， RNA序列S表示为x₁ x₂ x₃…x_n，x代表RNA分子的四种碱基：腺嘌呤A、胞嘧啶C、鸟嘌呤G、尿嘧啶U中的任一种，即x_i∈{A，C，G，U}，序列的长度为n，1≤i≤j≤n；因此，将RNA序列的每个碱基，用该碱基所在的位置序号来替代，即表示为1，2，3，…，i，…，j，…，n，这种编码方式被称为长度编码；

S102，根据待预测的RNA序列，判断RNA序列中的碱基配对情况，当RNA序列中i，j位置发生碱基配对时，将长度编码中其对应的位置进行交换，即此时的RNA序列长度编码为：1，2，3，…，j，…，i，…，n；根据碱基配对情况计算符合碱基配对规则的随机点匹配列表：

（i，j，k），

其中i，j分别表示RNA分子序列的第i个位置和第j个位置，k为随机点（i，j）的连续匹配数，即RNA序列长度编码为：

1，2，3，…，j，j-1，…，j-k+1，…，i+k-1，…，i+1，i，…，n；

S103，计算每组随机点（i，j）的连续匹配数列表：

该连续匹配数列表的计算方法为：根据随机点匹配列表中的任意两个数组（i，j，k）和（m，n，l），若i=m，j=n，k≠l，那么随机点组合（i，j）对应的连续匹配数列表中有k、l；

S104，判断当前系统是否达到临界温度T_临，若没有达到，进入步骤S105，若达到进入步骤S113；

S105，判断当前环中最小碱基数是否达到上限，上限为分子长度n-2*最小茎区数，若没有达到，则进入步骤S106，若达到，则进入步骤S112；

S106，根据随机点匹配列表，计算新的分子状态，具体过程为首先从随机点匹配列表中随机挑选一组随机点（i，j），然后从该组随机点的k连续匹配数列表中随机挑选一个连续匹配数k，根据挑选的随机点组合和挑选的连续匹配数生成新的RNA分子长度编码，即分子新状态；

组碱基数和假结总交叉分组数；

S107，根据步骤S106产生的分子新状态，计算出总的碱基配对数、平均每组碱基数和假结总交叉分组数；

根据下式计算新旧RNA分子状态的能量差：

∆E=E_新-E_旧（1）

其中：E_新为分子新状态的能量，E_旧为分子旧状态的能量，

(2) ，当碱基假结交叉分组数超过设定 的最大假结交叉分组数的情况时；

（3），

当碱基假结交叉分组数没有超过设定的最大假结交叉分组数的情况时；

其中，总分组数为随机点匹配列表中三元组（i，j，k）的个数，假结总交叉分组数为上述三元组两两交叉的个数；

S108，判断S107所得△E是否大于0，若大于0，则进入步骤S110，若小于0，则进入步骤S109；

S109，根据公式（4）进行计算：

（4），

其中，Random（0，1）为用随机函数Random生成一个0到1之间的数，若等式成立，则进入步骤S110，否则进入步骤S111；

S110，将RNA分子状态替换为新分子状态，更新RNA分子长度编码；

S111，环中最小碱基数自增，再次进入步骤S105；

S112，系统温度T以T=T₀*0.98的方式进行一次降温，然后，再次进入S104步骤中；

S113，输出RNA分子最终状态。

2.根据权利要求1所述的基于模拟退火的核酸结构预测方法，其特征在于，步骤S102中随机点匹配列表（i，j，k）中随机数i，j需要同时满足以下关系：

i＜j

j-i-k＞3

若随机数不满足上述关系，则重新生成随机数，若满足，则判断是否满足k连续匹配，若不满足k连续匹配，则重新生成随机点，满足，则添加到随机点匹配列表中。

Images