基于结构信息的RNA多序列比对.docVIP

　　基于结构信息的RNA多序列比对 【摘要】 本研究提出了一种考虑了结构信息的同源RNA多序列比对算法，它先利用热力学方法计算出每条序列的配对概率矩阵，得到结构信息，由此构造各条序列的结构信息矢量，结合传统序列比对方法，提出优化目标函数，采用动态规划算法和渐进比对得到最后的多序列比对。试验证实了该方法的有效性。 【关键词】 多序列比对；RNA二级结构；配对概率矩阵；结构信息矢量；动态规划 Abstract：ARNA［10］，RNAlign［13］，PMmulti［14］，这种方法先进行双序列比对，然后渐进的完成多序列比对。我们提出的方法属于后者。 2 算法 Sankoff［15］首先提出同时进行序列比对和结构预测，但是该算法的时间复杂度为O(N6)，空间复杂度为O(N4)，其中N为序列长度。已有的几个采用此方法的程序都使用了不同的限制，比如，Foldalign［16］利用了核心比对和贪婪算法，而Dynalign［17］则是通过限制两个序列间的最大距离来减少复杂度。 我们采用类似Sankoff算法的思路，但不是为了同时进行序列比对和结构预测，只是为了得到考虑结构信息的多序列比对。基本步骤是：首先，对每条序列，分别计算出其碱基配对概率矩阵，然后将这些矩阵变换成易于比较的结构信息矢量，通过两两比对这些矢量，构造出一个比对指导树，最后根据比对指导树，渐进的得到多序列比对。 2.1 碱基对配对概率矩阵 为了得到配对概率矩阵，首先要进行划分函数的计算，McCaskill［18］给出了RNA二级结构的划分函数的概念。RNA二级结构的划分函数Q定义为： Q=∑Se-△G(S)/RT(1) 式中，ΔG是结构的Gibbs自由能变化量，R是气体常数，T是绝对温度，S 是所有可能二级结构的集合。McCaskill提出了一种动态规划算法来确定二级结构形成中的划分函数，该算法给出了序列中每个可能碱基对的配对概率，用一个概率点图显示，程序RNAfold［20-21］就是采用的这种算法。因为对能量规则进行了简化，对多分支环的处理是用单链碱基的自由能来模拟碱基堆积间的能量增量，所以这种算法可能会模糊两个候选结构间的差别。MatheB(i,j)，BL(i,j)，Wcoas(i,j)，W5(i)，W3(i) 以上矩阵分别对应不同的分支环情况，W5(i)是从序列的5’端到包括碱基的这一段序列的划分函数。W3(i)是从碱基i（包括碱基i）到3’端的这段序列的划分函数。 由以上定义可得，碱基i到j配对的概率Pij为： Pi,j=∑ke-△G(k)/RTQ=1Q∑ke-△G(k)/RT= V(i,j)×V(j,i+N)W5(N)(2) 2.2 概率矩阵比较 有了碱基配对概率矩阵，就可以进行比较了，首先考虑两条序列的情况。文献［14］给出了比较两个概率矩阵的方法，我们用了其中相关定义。给定两条序列S1和S2，通过（2）式计算出其碱基配对概率矩阵P(S1)和P(S2)，这里包含了结构信息，可以通过比较这两个矩阵来计算他们之间的相似性。换而言之就是要找出他们之间共有的结构，并用权重来表示他们之间的相似程度。我们要找出序列S1中和序列S2中所有相似的碱基对（序列S1中碱基对为(i,j)，序列S2中的碱基对为(k,l)），并满足下面的 条件： ∑matches(ij;kl)〔Ψ(S1)ij+Ψ(S2)kl〕→max(3) 式中Ψ(S1)ij是序列S1中碱基i,j配对的权重，可以用配对概率来计算： Ψij=log(Pij/Pmin)(4) pmin为有意义的最小配对概率。对一般化的序列比对情况，我们考虑序列S1和S2有Ngap个空位，同时考虑共同结构S时，目标函数为［14］： ∑(ij;kl)∈S(Ψ(S1)ij+Ψ(S2)kl+τ(S1(i)，S1(j)；S2 (k),S2(l)))+γNgap+∑i∈S1，k∈S2Sδ〔S1(i),S2(k)〕(5) 要使它为最大。式中γ＜0是空位罚分。分值δ(S1(i),S2(k))和τ(S1(i),S1(j);S2(k),S2(l))分别表示未配对和配对碱基的替代情况。在最简单的情况下，我们不考虑序列部件，设σ=τ=0。 用Mi,j;k,l表示子序列S1(i……j)和S2(k……l)间最好的比对分值，用M′i,j;k,l表示碱基i,j和k,l构成碱基对的情况下最好的比对分值。这样我们就可以用动态递归来计算比对分值了，递归公式如下［14］： Mi,j;k,l=MAXMi+1,j;k,l+γ Mi，j;k+,l+γ Mi+1，j;k+1,l+δ(S1(i),S2(k)) MAXh≤j,q≤l｛M′i,h;k,q+Mh+1,j;q+1,l｝ M′i,j;k,l=Mi+1,j+1;k+1,