计算生物学中相关问题算法及思考
1 背景介绍计算生物学是指研究利用数据分析、数据建模和计算机技术等来进行生物学数据的分析、生物学结构的预测等方法的学科。随着生物学研究的不断深入，其需要处理的数据也呈很快的速度增长。就其学科涉及的内容而言，涉及数学、化学、物理、生物、统计、计算机等传统学科。就计算机方面所要研究的内容来讲，就是对通过各种方法获得的数据进行处理，或设计一些对一些生物分子、生物组织结构的形态规律的算法，以期利用计算机的高速计算性能，帮助分析数据、预测数据。
本文主要讨论了RNA 结构预测的算法，并在后面提出了自己的想法。
2 关于RNA 结构
一级结构是指RNA 单链中由四种核某酸排列而成的有限线性序列。一级结构定义：RNA 一级结构是取自字母表Σ ={A，U，G，C} 的字符串R。即R=r1r2r3…rn，其中ri∈Σ。二级结构是指RNA 序列通过自身回折而形成部分碱基配对和单链交替出现的茎环结构。不同核酸之间的配对规则：C-G(G-C)，A-U(U-A)，G-U(U-G).二级结构定义：用(i,j) 表示ri 和rj 构成基对，1 ≤ i ≤ j ≤ n。设S 为其一级结构的集合，称S 为RNA 二级结构，如果S 满足：(1) 对于任意(i1,j1) 和(i2,j2)，要么i1=i2 且j1=j=, 要么i1 ≠ i2，i1 ≠ j2，j1 ≠ i2 且j1 ≠ j2。(2) 如果h<i<j<k，则S 不能同时包含(h,j) 和(i,k)。(3) 如果S 包含(i,j) 则|i-j| ≥ 4。三级结构是由各二级结构单元之间相互作用并在空间中形成稳定的地位和取向而构成的。
3 其他一些算法的介绍
3.1 随机文法
在随机文法方法中，把RNA 序列看成具有一定语法规则的语句，通过这些语法规则来解析RNA 序列中存在的碱基配对关系，也就是它的语义，从而得到该序列的二级结构。该方法是一种基于已有序列的先验知识的方法。该方法可以通过扩展随机上下文无关文法来预测RNA 二级结构中假结的存在。
该方法需要解决的三个问题：(1) 对一个未知序列与用一套样本序列训练生成的模型进行比对，判定该序列在模型中走哪条推导路径；(2) 在训练好的模型中计算生成一条给定序列的概率；(3) 用一套样本序列( 或结构) 对随机文法模型的最优概率参数进行估计。
3.2 基于堆积能量和协变信息的预测方法( 念心算法)
给定一个RNA 分子的碱基序列R，基于能量的动态划分算法的思想是：首先找出整条链中最稳定的结构，然后，在保留这个结构的前提下逐步找出余下部分的最稳定结构．具体说来，算法分为两部分，一部分是划分算法，另一部分是寻找最小能量茎区的算法．通俗来讲，RNA 二级结构是由一个个茎区组成的，如果所有茎区确定了，那其二级结构就确定了。该方法通过 迭代得到二级结构的茎区集，每 次迭代寻求具有最低能量的最大茎区，存储，并标记使其不参与以后的迭代。迭代结束后得到茎区集，从而确定RNA 二级结构。
4 GTfold 算法介绍
关于能量定义：考虑到茎区对分子的结构起着稳定作用。而环会降低结构的稳定性，Estack 为负值，E-hairpin、E-bulge、E-interml 和E-multi 为正值，原始一级序列的能量值为E-o=0。
该算法求出可能的RNA 二级结构中最有的最小能量，然后回溯得到最小能量。
5 算法的改进
a. 并行化计算
在使用原始算法的情况下，注意到VBI，VM，WM 填充过和中，填充每个元素需要的基础要么是其他表提供的，要么是自己表内左下方的元素，这说明整个工作是可分割的，即在每一个部分都可分割成相互不依赖的部分。一种简单的分割方法如下图所示，方格中数字表示当前计算的步骤。
每个数字对应几个区域，表明在该步可用几个CPU 并行计算。这样有一个缺点会导致在后面的计算中能用的CPU 越来越少。一个解决方法是重新划分工作份额。下图表示了一种可能的划分，但并不一定是最优。划分发生在第3 步后，用第4 步将工作重新分成了6 份。
具体在何时重新分配，将是一个分配所需工作量和实际计算工作量之间的平衡。
b. 内部环加速算法
由上面的分析，在不舍弃运算精度的情况下，内部环运算需O(n4) 的时间和O(n2) 的空间。
我们重新考虑内部环计算时的情况。实际上，在计算VBI(i,j) 时，我们考虑了从i+1 到j-1 中任一对碱基的配对情况，即利用了V 矩阵中(i,j) 左下方的所有点。利用结论△，我们可以首先考虑内环两条单链长度都小于等于一常量c 的情况，由这里算出的最小值去更新所有满足VBI 中的(i’,j’) 的值，其中1《i’<i 且N>j’>j。这种算法的效率比原始算法要高出不少，因为更新只发生在环中两条单链长度恰为c 或一边为c 另一边为c+1 时。
6 一些其它的想法
诚然，GTfold 在效率和预测准确性上有很好的表现。但在哪些不足导致它的预测与实际结构有一定偏差？一种可能的是解释是：在实际条件下，RNA 并不能总达到最低的能量，而是处在一个极低能量的结构上。但是，从极低到最大低的变换过程可能需要经过一个比较高能量的状态，这导致了在实际中这种变换很难完成。一方面，如果RNA 分子的没有从外界获取能量，那么它就不能跳跃着变化。（隧穿效应是微观世界的，RNA 显然太大了。）虽然有的结构能量不是最低，但实际上可能更接近实际结构。另一方面，如果我们假设RNA 可以吸收能量而处于一个比较高的不稳定状态，然后跃回最其它的极低能量状态。
我们做出三种可能的假设：1、RNA 分子是转录完再折叠的。2、RNA 分子是边转录边折叠的。3、折衷地，RNA 每转录k 长度重新进行一次折叠。估计第一种情况应该简单一些。不论哪种情况，算法基本思想应是一致的。
我们从另一个角度出发，不再考虑分子的最低能量，而考虑分子从基本状态走哪条路径，其自由能下降最快。这样找到的结果未必是最低自由能，但却是实际情况下最可能发生的结果。
如果我们考虑一个阈值E，它是RNA 可能偏移当前极低能量状态的最大能量值，那么，如果状态B 的能量小于状态A，但从A 到B 有一条路径（路径上每一点对应一种RNA 二级结构，能量变化来自打开碱基再重新配对，相邻点间只需进行“一次”结构变化即可到达），满足该路径上任两点间的差不大于E，则我们认为A 是可以到达B 的，应该将B 作为最可能结构。
另外，查阅相关资料，假结对RNA 的功能有很大的影响。
GTfold 如果能将假结的相关算法巧妙加入其中而不带来复杂度上太大的变化，将是一件很值得期待的成果。这一点，我们现在还没有多少可行的想法。__
谢谢分享！:D
页:
[1]