RNA 二级结构预测方法

RNA 二级结构预测是指以生物信息学理论为基础，输入 RNA 一级序列，计算预测 RNA 二级结构的过程。RNA 的一级序列主要由 A，G，C，U 四类碱基组成。RNA 二级结构指由非相邻碱基相互作用折叠形成的平面结构。
最近的研究表明，RNA 的二级结构对于调控来说具有非常重要的作用。例如，RNA 的不同结构模式与编码区域、剪接点和多腺苷酸化点相关；动物细胞的活动过程和酵母受 RNA 二级结构的影响明显不同；又如，RNA 二级结构在蓝藻多变鱼腥藻的 nifH1 转录稳定性中发挥了重要作用。同时，有证据表明二级结构可以极大地影响三级结构。所以确定 RNA 二级结构，对分析和理解基因遗传信息的传递机制至关重要。

下面将描述一些传统的 RNA 二级结构预测算法。

1、比较序列分析方法

比较序列分析法也是常用的 RNA 二级结构预测算法之一。其工作原理其实比较简单，它以 RNA 序列中互补碱基间的共变联配（英文名称为 covariant-alignment）活动为基础；以已知的 RNA 序列的数据为依据标准，以查找被测算 RNA 序列中的高近似度序列为手段；以一定的相关数学模型为依托，共同研究推算所给 RNA 序列的二级结构。
比较序列分析适用于多条 RNA 序列的结构预测问题。在生物实验过程中，常常需要同时处理一组或几组同源的 RNA 序列。通常认为序列的保守性要小于结构的保守性。因此即使这些同源 RNA序列的长度、碱基排序不同，结构却十分相似，在生物体内也发挥出相似的功能，例如，tRNA 分子的序列虽不一致，但其二级结构都是三叶草形的，三级结构是倒 L 形。
由于结构的保守性，比较序列的方法有着较高的预测精度。比较序列分析法主要由多序列比对与序列结构预测两个部分组成，根据序列比对和结构预测顺序的不同，比较序列分析的 RNA 二级结构预测算法可以分为如下三类:

（1）先比对后预测，首先使用多序列比对工具来构造出 RNA 序列比对结果，然后通过检测碱基突变等信息来揭示这些 RNA 序列共有的一致保守结构。然而，使用该类算法会得到多个可能的预测结果，且无法保证其中是否包含真实的二级结构。
（2）先预测后比对，首先对集合中的每个序列分别进行结构预测，然后对这些结果进行比较分析，最后获得这些 RNA 序列的保守二级结构。这种预测算法过于依赖结构比对的结果，同时准确性也无法得到很好的保证。
（3）序列比对与结构预测同时进行，算法利用多序列比对方法和最大碱基配对算法递归的完成序列比对和 RNA 结构预测等操作。由于序列比对与结构预测需要同时进行，其空间复杂度和时间时间复杂度明显增加。因此，这类算法在预测的序列长度不能太长，否则花费的成本太高，序列长度较小时可以尝试使用。

比较序列分析方法是目前测算分子结构准确度最高的方法之一，同时，比较序列分析方法在测算假结结构等问题时也有比较突出的结果，然后目前正在使用的其它方法中有的还仍不能对假结进行测算。但比较序列分析方法也有不足的地方，就是对于较少、甚至是一条序列，以及同源性低的序列不适合使用，测算结果较差。

2、动态规划法

动态规划是一种寻求最优决策的数学方法，其核心思想是将一个多阶段决策问题分解为一组子阶段决策问题，利用各个子阶段的关系对决策问题逐一求解。目前，动态规划在管理、计算机、生信等领域已有广泛的应用。
动态规划类似于分治法，都是通过将原问题分解成子问题，然后对子问题的解组合来求解出原始问题解的算法。但两者的差异是，分治法将原始问题分解为多个不重叠的子问题，分别求出这些子问题的解，然后结合起来作为原问题的解；动态规划则适用于子问题相互重叠的情况，即两个不同的子问题间具有相同的子子问题。动态规划对于分解后每个子子问题仅求解一次，将子子问题的结果保存在表格中，从而不必在每次求解这些子子问题时重复计算。减少了很多额外的计算工作。因此，动态规划本质上是一种以计算空间换取计算时间的技术，在运算的过程中它需要存储各种状态，因此其空间复杂度要大于其他算法。

2.1 最大碱基配对算法

1978 年 Nussniov 提出了最大碱基配对算法。最大碱基配对算法基于在碱基互补配对的过程中碱基间的氢键能让两个碱基比较紧密地结合在一起的基本假设，RNA 结构中配对的碱基对越多，连接的氢键越多，结构就越稳定。最大碱基配对算法认为 RNA单链折叠成碱基对数目最多时的状态即为该 RNA 的真实二级结构。简单说法为，RNA 单链的自我折叠使其碱基尽量达到最大互补配对时，就构成了 RNA 的二级结构。算法在求解 RNA序列的最大碱基配对数时使用了动态规划的思想，将整个问题分解成多个小的子问题去解决。
Nussinov 算法是最早使用动态规划进行 RNA 二级结构预测的算法，但算法的预测精度却不高，这是由于 Nussinov 算法只是简单的计算了碱基对数量对结构的影响，没有考虑到各种构件的能量的不同、碱基对类型的不同等其他因素也对序列的折叠起到的作用。另外由于最大碱基配对算法所预测的各个碱基对相对独立，没有考虑到连续碱基对可以形成茎区这一更为稳定的结构，因此预测出二级结构中各个碱基对是不连续的，不能够形成稳定的茎区。

2.2最小自由能算法

最小自由能算法(Minimum Free Energy,MFE)是由 Zuker 所提出的，与 Nussniov 算法简单计算碱基对最大配对数不同的是：Zuker 将从实验得出的能量数据引入 RNA 二级结构预测的计算中。RNA 分子中碱基对以及核苷酸之间均是靠氢键进行连接的，通过破坏 RNA 二级结构中不同子结构中的氢键所需要的能量值可以得到这些子结构在组合时所需要消耗的能量，这个能量被称为自由能。
RNA 分子结构之所以能够稳定，其原因就是碱基配对造成了能量的降低。因此，最小自由能算法认为，自由能趋向最小[30]的时候为 RNA 真实的二级结构。这是因为，自由能趋向最小的时候 RNA 分子会通过构象调整从而达到某种热力学上的平衡进而形成了 RNA 分子最稳定的状态。最小自由能（Minimum Free Energy，MFE）算法的中心思想[25]为，RNA 的二级结构在最优化组合时，通过分子动力学方法获取的各种茎环结构的总能量最小，结构也最稳定。

经实验验证，形成碱基互补配对的茎区结构的自由能为负值，而凸环、内环等未形成碱基互补配对的单链结构的自由能为正值。因此，在 RNA 的结构中，茎区结构越长，数目却多，整体结构的自由能也越小。而研究人员普遍认为，当生物分子具有的能量越小时，其结构也越稳定的。因此，最小自由能算法的核心就是将预测出来的自由能最小的结构作为 RNA 的真实结构。
自由能并不是和碱基对的数目成简单的线性关系，否则最小自由能算法就退化为最大碱基配对算法。自由能不仅和碱基对的个数、类型相关，而且茎区、环状结构等不同基本构件的自由能也有较大的差异。在此基础上，Zuker 等人将 RNA 二级结构划分为茎区、凸环和内环等不同的基本构件，通过实验的方法得到每种构件各自的自由能参数。Zuker 等人在计算 RNA 二级结构整体自由能大小时假定组成整体结构的各个构件相互独立，整体的自由能即为各个构件自由能之和，利用动态规划的方法对 RNA二级结构中的各个构件进行组合，可以得到整体自由能最小的二级结构。由于 Zuker算法需要计算出每个基本构件的自由能数值，因此最小自由能算法在时间复杂度上要高于最大碱基配对算法。
最小自由能算法所预测的结构中各个配对碱基对是嵌套或顺次的，即不能预测结构中存在的假结，Rivas 和 Eddy 在最小自由能算法的基础上，使用隔空矩阵来计算带假结的 RNA 序列的二级结构。该算法不仅计算了茎区、凸环、内环等基本构件的自由能，还计算了产生假结时产生的自由能数值。该方法虽然能够预测出假结，但也提高了算法整体的时间复杂度和空间复杂度。

3、组合优化方法

RNA 二级结构预测组合优化的方法主要思想是，根据碱基配对可以构成各式各样不同种类的茎区，茎区组合的种类繁多茎区数量大，因此该问题灵活性很大，也容易产生错误。从数学的角度来看，只要知道了这些茎区怎样进行选择以及如何组合，就能够有效地预测出 RNA 的二级结构。该类方法有两种最具代表性的算法：一种是螺旋区堆积算法（helical regions stacking methods）这种算法的基本观点是代表了一类以统一的结构单元来对待 RNA 的二级结构组成。另外一种常用方法是最大权重匹配算法（MWM,maximum weighted matching algorithm）这种算法其实是一种图论的算法，遵循最大权重匹配原则。
使用组合优化这种类型的方法来预测 RNA 二级结构与传统的预测方法的差别非常大，它们的基本思想都非常新颖独特。目前来看这类方法存在着一个共同的缺点：无法保证能够靠近全局最优解，更无法保证找到全局最优解；对于预测出来的结果也没有判断优劣的手段。RNA 二级结构预测中茎区的选择和组合是非线性问题，属于离散问题，所以利用传统的组合优化的方法并不容易解决，虽然已经提出了例如螺旋区堆积算法、最大权重匹配算法等方法，而且可以预测假结，但是在实际的应用中对假结误差较大，其预测的能力也受到了诸多方面的限制。综上，该方法中的各算法往往被作为一种参考，或者与其它算法结合在一起来使用。

3.1 螺旋区堆积算法

早在 1975 年，PipaS 和 Mc Mahon[37]就首先提出了螺旋区组合算法来预测RNA 二级结构。螺旋区堆积算法，是将 RNA 的二级结构看作是一组特定的碱基配对集合，这里的螺旋区其实就是我们平常所说的茎区。已知一条 RNA 序列，倘若明确了其碱基配对关系，即确定了其二级结构。该算法的主要操作方法就是尽最大可能列出所有螺旋区的所有组合情况，来达到找到其最小自由能的值为最小的组合结构，这个结构就是最稳定的结构状态。当然，在选择螺旋区的时候，以及将螺旋区组合的时候还要遵循一定的规则。
具体来说，螺旋区堆积算法预测 RNA 二级结构可以分为如下三个步骤来完成：
（1）找出所有可能螺旋区。针对所要预测的 RNA 序列，构建两个矩阵“A”、“B”，“A”为邻接矩阵，“B”是“A”的相容矩阵，通过“A”矩阵找出所有稳定的螺旋区，根据结果建立起螺旋区列表。因为“B”矩阵为“A”的相容矩阵，因此就可以利用这两个矩阵来判断两个螺旋区能不能同时存在于一个 RNA 的二级结构中。当两个螺旋区当中没有共享的碱基，则两个螺旋区是相容的，否则就是不相容的。
（2）把所有已经已知相容的螺旋区从全部的螺旋区中挑选出来，统一组成一个结构形成螺旋区的序列，要注意的是其螺旋区的数量不能少于三段。
（3）对于所预测出的 RNA 二级结构进行评价。评价的标准依然还是常用的最小自由能方法。求解自由能的方法是利用碱基能量函数求解，把所有结构的自由能都用这种方法求解出来，然后从大到小做一个排序，其中自由能最小的结构，就是我们要找的 RNA 二级结构的最优解。

3.2 最大权重匹配算法

最大权重匹配算法是 RNA 结构预测模型上新兴起的一种算法。它在预测假结合与三级结构的相互作用方面有着重要作用。
最大权重匹配算法是一种图论算法，算法中将 RNA 序列比拟成一个环形，利用了非二分图的最大权重匹配算法。算法的主要思想为，RNA 序列中的每一个碱基代表一个顶点，由这些顶点组成一个环的形状，并根据 Watson-Crick 配对规则和与 G-U 配对规则，在各碱基对中间，用一条线把已构成配对的碱基连接起来；当每一条线都被赋予权重值后，就可以在刻画出所能够出现所有的，真实 RNA 折叠结构中碱基配对的可能性，即每一种匹配线组合都是一个 RNA 的可能存在的结构。在整个匹配的过程当中，不允许出现一对多的碱基配对关系，只可以在两两碱基之间形成一对一配对。在其折叠图中 RNA 序列最优的二级结构就是总权重值最大的那个匹配组合。先寻找到含有一条线的最大总权重值，随后再寻找含有两条线的最大权重的匹配，通过不断来增加总权值一直继续扩展操作下去。
最大权重匹配算法最突出优点是：可以预测多种类型的 RNA 结构，例如假结、三联碱基配对等，并且可以获得全局最优的结果，其算法简单、思路清晰。但是，最大权重匹配算法的缺点也同样在于此处，它一直追寻最后匹配的最大总权值，导致了许多假性的碱基配对（spurious base pairs）也被当成解预测出来。

4 启发式方法

启发式算法是相对于最优算法提出的，其基本思想为在可接受的计算费用范围内搜索最优解，启发式算法不能保证所得到的解是否可行以及是否是最优解，甚至许多数情况下，其所得到解与最优解之间的近似程度都无法表示出来[42]。通过概念可以看出，使用启发式算法来预测 RNA 二级结构我们也无法去判断得到的结果是否是最优解或者说离最优解是不是最近。正因为如上所述，这种方法在使用的时候要考虑与其他的方法结合在一起，以弥补其不足，或者是用其得到的解作为一种参考。

4.1遗传算法预测

使用遗传算法对 RNA 二级结构进行预测，寻找出所有可能存在的茎区是第一步要进行的工作，将找到的茎区形成茎区池，然后进行初始结构（每一个构建的结构相当于遗传算法中的一个染色体），形成具有多个二级结构的初始结构群体。
在构建结构的过程中要注意所添加的茎区彼此之间的相容性。采用自由能作为构造群体的适应度标准一代一代的进化下去，能量越小其结构越好。 突变操作就是对结构中的茎区替换，把某些茎区从父代结构中移出，然后在茎区池中重新选择一个其它的茎区添加回去，可以使解搜索有机会跳出局部最小。按照轮盘旋转法自由能越大的茎区越有可能被移出。茎区的加入则根据随机性选择，随机性选择出来的茎区要通过 Metropolis 准则来判断，不是直接就可以添加的，如果不符合准则则需重新随机挑选，直到符合了原则，才能够替换成功。在替换的过程当中要注意保证茎区之间的相容性。
交换为把作为父本的两个结构（注意两个副本不可以是同一个结构）从一个或者多个位点进行交换。把两个父本的茎区合并成一个新的茎区池，根据池内茎区各自的自由能分支按照轮盘旋转的方法来选，构造一个新的结构。适应度同样也作为选择两个父本结构的准则依据，其值如果越高那么意味着被选取的概率就越大。这样操作可以使子代继承父代中结构较稳定的结构单元。
如果初始群体当中有 N 个个体，交换和突变操作之后会分别又产生 N 个新个体，候选个体的总数量变为3N 。每一代群体都是上一代群体遵循适者生存，优胜劣汰的原则所挑选出来的，这样一代一代的循环，就可以保证适应性强的个体有更高的概率和更多的后代被保留下来。由于遗传算法是从一个初始群体开始计算，而不是个体，导致了具有隐含并行搜索的特性，也是因为如此遗传算法容易过早收敛陷入局部最小，所以需要重点考虑解决这一问题。
遗传算法与其它经典算法具有很多优点：
（1）遗传算法以适应度作为搜索信息，不需要其它信息，而且不容易落入局部最小点。
（2）遗传算法可在解空间内平行搜索，而不是点对点搜索，可以有效地搜索空间中的许多区域，正因如此也是对局部最小值不敏感的原因。
（3）遗传算法代表了潜在的解，不是解本身，因此不需要完全理解问题模型，只需要去评估试验解的正确性即可。

4.2 遗传模拟退火算法

模拟退火算法（Simulated annealing，SA）最早的思想是由 Metropolis[48]等人在 1953 年提出的，基本思想是模拟统计物理中固体物质的结晶过程和一般的组合优化问题之间所具有的相似性质，基本原理是模拟热力学中粒子系统的降温过程求解规划问题的极值。模拟退火算法是一种非常通用的优化算法，整个模拟物理退火的过程总共有三个基本组成部分分别为：加温的过程、等温的过程和冷却的过程。将需要组合优化的待解决问题看成是一个物理系统，将目标函数比拟为物理系统的能量。把寻求组合优化问题最优解模拟成为物理系统逐步降温并且达到最低能量状态的退火过程。
将模拟退火算法和遗传算法结合来预测 RNA 二级结构可以互相弥补两种算法中的缺点和不足取得更好的效果。 在使用遗传模拟退火预测 RNA 二级结构过程与遗传算法相比在中最大的区别是，按照 Metropolis 准则，使用模拟退火思想来选择或淘汰茎区。在种群中，个体的之间所发生的变异和交叉操作与 Boltzmann 分布是非常接近的，而Boltzmann 分布却很容易陷入局部最小，因此可以通过降低温度遗传算法来对此加以改善。
退火算法具有更好地收敛到全局最优的能力。遗传算法因为有其并行搜索的特性所以在对搜索空间的搜索过程中具有广泛灵活的优点。把这两种算法结合在一起退火算法由于其具有收敛性好，稳定性好的优点补充了遗传算法易陷入局部最小的缺点，而遗传算法对搜索空间的灵活广泛性是退火算法提高了向最优解求解的效率。利用这两种算法结合起来一起应用于 RNA 二级结构预测，为 RNA 二级结构预测提供了除了最小自由能方法之外的又一种比较理想的方法，而且也进一步改善了单独使用遗传算法、或退火算法去预测 RNA 二级结构的效果，多种方法结合从不同角度预测 RNA 二级结构也是以后发展的必然。

5、常用软件的比较

参考文献：
RNA二级结构预测算法的研究_邢翀
基于卷积神经网络的RNA二级结构预测方法研究_张春鹤