[发明专利]一种基于并行改进人工蚁群算法的RNA二级结构预测方法

权利要求书

1.一种基于并行改进人工蚁群的RNA二级结构预测方法，其特征在于所述方法具体过程按如下步骤进行：

步骤一：构建RNA二级结构茎区池；

步骤二：基于并行人工蚁群进行计算最小自由能并构建相应的二级结构，避免陷入局部最优以获得全局最优解；

步骤三：对构建好的二级结构进行可视化。

2.根据权利要求1所述一种基于并行改进人工蚁群的RNA二级结构预测方法，其特征在于：所诉步骤一中根据RNA序列S和茎区最小长度n，构建茎区最小长度构建茎区池：

1)根据序列S的长度N，构建N*N大小的茎区矩阵，并将序列按次序放在首行的上一行，首列的前一列，然后在下三角矩阵中，根据碱基对匹配的规则，将矩阵中匹配碱基的对应位置设为1，将矩阵中不匹配碱基的对应位置设为0，碱基包括A(腺嘌呤，adenine)、G(鸟嘌呤，guanine)、C(胞嘧啶，cytosine)和U(尿嘧啶，uracil)，根据Watson-Crick和GU摇摆碱基对配对法有六种配对方法：“AU”、“UA”、“CG”、“GC”、“GU”、“UG”，以RNA序列CGCCCAGCGAAAUGCAAAGUC为例。

2)茎区矩阵建立好之后，在茎区矩阵中寻找连续的碱基对，其过程如下：

a)从矩阵的左上角第i(0≤i≤N)行第j(0≤j≤N)开始，向矩阵右上(i+1，j-1)的方向搜索连续的碱基对，并在首个连续处停下。

b)计算连续碱基的长度是否大于茎区最小长度n，如果连续碱基的长度大于n则将此段连续的碱基保存起来，反之，则不保存。

c)继续向茎区矩阵的右上方寻找连续的碱基对，如果遇到匹配的碱基则执行b)。

d)如个寻找的位置到了矩阵的边缘，判断的条件为i＝N或j＝N，就将i和j重置为原来的位置，并将i设置为i+1，以RNA序列CGCCCAGCGAAAUGCAAAGUC，n＝3为例，其茎区矩阵中连续如图1中所圈出碱基。

3)完成茎区的寻找后，需要将寻找到的茎区作为茎区池保存，并作为下一步的输入。

3.根据权利要求1所述一种基于并行改进人工蚁群的RNA二级结构预测方法，其特征在于：所诉步骤二中一种基于并行改进人工蚁群的RNA二级结构预测方法，避免陷入局部最优以获得全局最优解；具体过程为：

1)从设置人工蚁群算法的参数，控制算法的参数包括种群大小colony_size、最大迭代次数steps，信息素蒸发的常数ρ，信息素初始值initial_pheromone，信息素的权重α，启发函数的权重β，启发函数η(i)，线程数量m，信息素矩阵初始化全部为initial_pheromone；信息素矩阵的大小为与茎区池中茎区数量一致，将每只人工蚂蚁分配到不同的线程上，每只人工蚂蚁内部通过存储一个bitmap结构，判断对应位置的茎区是否被访问，初始化时，将所有的bitmap设为false，当访问过某一个节点后设为true，设置为true之后不再重复访问。

2)当判断两个茎区是否存在冲突和假结，对茎区A和茎区B而言冲突的判读如公式1所示

对于假结的判断如公式2所示

i＜i′＜j＜j′(i，i′∈Stem_A and j，j′∈Stem_B) (2)

3)判断茎区i是否与解决方案S冲突的三个条件

1.茎区i与解决方案S中每一个茎区都不冲突。

2.茎区i与解决方案S中每一个茎区都不构成假结(如果有需要)。

3.茎区i没有被蚂蚁k访问过。

4)将所有符合(2.3)要求的茎区构造集合N_k。

5)茎区选择算法的公式如公式3所示：

其中p_k(i)表示第k个蚂蚁选择第i个茎区的概率，τ是表示信息素，η表示启发函数。

6)每个节点的信息素更新算法为公式4所示

其中表示新的所有人工蚂蚁经过的结点的信息素，Δτ^k表示单个节点新的信息素。

7)每个节点信息素的挥发遵从公式5：

τ(i)＝(1-ρ)·τ(i) (5)

每次信息素都会在原有的基础上按照一定比例挥发。

8)在RNA序列的实际折叠过程中长度更长的序列更容易被选中，启发函数的定义见公式6

其中i表示某个茎区，i.length表示茎区序列的长度，paired_length表示所有匹配碱基的长度。

9)在公式3为每个茎区生成了不同的概率，再通过轮盘赌算法返回通过概率计算得来的茎区。

10)将计算得来的茎区添加到解决方案S中。

11)计算解决方案S的自由能，自由能的计算如公式7所示：

E＝E_helices+E_hairpin+E_bulge+E_intermal+E_multibranch+E_pknots (7)

RNA二级结构的自由能根据每一个部分不同计算得来，所以首先需要使用区分算法区分不同的结构。

12)如果单链末梢两端碱基与同一个茎区a相邻，则这条单链s与茎区a属于同一个发夹环子结构，s属于发夹环子结构的环状单链，茎区a属于发夹环的茎区，如公式8所示：

s_1-1，s_n+1∈StemA (8)

13)如果单链s的末梢两端所对应的碱基b₁，b₂，都是碱基对，且b₁和b₂所配对的碱基相邻一个单元，那么就认为单链s与b₁ b₂所在茎区形成一个凸环子结构，如公式9所示。

b₁₊₁＝b₂(b₁＝pair(s₁)，b₂＝pair(s_n)) (9)

14)对于一条单链s如果s末梢两端的碱基b₁、b₂是有配对，与b₁、b₂配对的碱基分别记做b′₁和b′₂，如果碱基b′₁、b′₂之间不存在其它的碱基对，则称b′₁和b′₂之间单链为s′，内环子结构包含单链s、s′、b₁和b₂所在茎区，如公式10所示。

b₁＝pair(b′₁)，b₂＝pair(b′₂)(b₁，b₂∈s，b′₁，b′₂∈s′) (10)

15)从单链角度出发给出多分支环的定义，对于一条单链s₁来说，如果这条单链s₁的末梢两端所对应的碱基b₁，b₂，都是碱基对，并且任意一个碱基对，假设b₁，令b₁的配对碱基为b′₁、b′₁必然与一条新的的单链s₂、s₂的长度length(s₂)≥0，单链s₂的末梢两端的碱基为b′₁和b₃，b₃又有配对碱基b′₃，像这样继续下去，不断出现新的配对碱基和新的单链，直到有新的碱基b_n有配对碱基b′_n，并且b′_n＝b₂。

16)对于假结的判断较为复杂，参照(2.2)中的方法。

17)螺旋区自由能的计算方法为公式11：

E_helices＝E_initiation+E_AUpenalty+E_symmetry+E_stack (11)

其中内分子的初始自由能一般为常数E_initiation，每个以AU碱基对作为结尾的惩罚项的自由能为常数E_AUpenalty，螺旋茎区的自对称自由能为常数E_symmetry，堆积的螺旋茎区自由能E_stack。

18)发夹茎区自由能的计算方法为公式12：

E_hairpin＝E_initiation(n)+E_mismatch+E_penalty(C) (12)

其中初始化自由能E_initiation(n)与内环的单链长度有关，而E_mismatch是指螺旋茎区和发夹单链之间第一个不配的碱基的能量，E_penalty(C)是惩罚项。

19)凸环自由能的计算方法为公式13、公式14：

E_bulge(n＝1)＝E_initiation(n)+E_penalty(C)-RTln(n) (13)

E_bulge(n≥1)＝E_initiation(n) (14)

计算凸环的能量需要分为两个部分，如果凸环单链长度唯一，则使用公式13，否则使用公式14，E_penalty(C)为惩罚项，RT是固定系数，E_initiation(n)与单链长度有关。

20)如内环来说，计算公式为公式15

E_intermal＝E_initiation(n)+E_asym*|n₁-n₂|

+E_mismatch(mismatch1)+E_mismatch2(mismatch2) (15)

其中E_asym*|n₁-n₂|为非对称项乘以内环中两条单链的差的绝对值，E_mismatch(mismatch1)与E_mismatch2(mismatch2)分别代表茎区与两条单链临界处第一个不匹配的碱基的能量值。

21)多分支环来说，计算公式为公式16

E_multibranch＝a+b×[averageasymmetry]+c×[numberofbranchinghelices](16)

其中a，b，c为系数，averagea symmetry代表每个螺旋茎区两侧单链碱基数量差值的平均值，number of branchinghelices是指多分支环中所有单链碱基的数量和。

22)对于假结来说，如计算公式17，18，19所示

E_pknots＝E_initiation+E_pairedbase+E_unpairedbase (17)

E_pairedbase＝bandpenalty+stacked_energy*paired_penalty (18)

E_unpairedbase＝unpairedpenalty*n (19)

其中paired_pena，band_penalty都是惩罚项，unpairedpenalty为为匹配的惩罚项。

23)使用公式7计算完自由能E之后使用公式20判断解决方案S是否接受这次添加：

S_old＝min_energy(S_old，S_new) (20)

min_energy会返回两个解决方案S中的更小者。

24)如果迭代次数达到steps，则终止算法，输出最小自由能的结果，否则回到2)继续循环循环。

4.根据权利要求1所述一种基于改进人工蜂群算法的路径规划方法，其特征在于：所诉步骤三中将解决方案S进行可视化；具体过程为：

1)将(2.22)求解出最小自由能的解决方案S转换为点括号的形式，匹配的碱基中位于前面的碱基用”(”表示，匹配的碱基中位于后面的碱基用”)”表示，没有被匹配的碱基用”.”表示。

2)将上一步生成的点括号表达式绘制成图形输出。