【摘要】 本研究提出了一种考虑了结构信息的同源RNA多序列比对算法，它先利用热力学方法计算出每条序列的配对概率矩阵，得到结构信息，由此构造各条序列的结构信息矢量，结合传统序列比对方法，提出优化目标函数，采用动态规划算法和渐进比对得到最后的多序列比对。试验证实了该方法的有效性。
【关键词】 多序列比对；RNA二级结构；配对概率矩阵；结构信息矢量；动态规划
Abstract：We presented a RNA sequences multi-alignment method based on structural information. Firstly， we computed base pairing probability of every sequence by thermodynamic method. Secondly， the structural information vector was constructed through gotten structure information and been pair alignment each other， as result， a guide tree was constructed. Finally， combine traditional sequence alignment， we presented the objective function and got the final multi-alignment by dynamic programming algorithm and progressive alignment with guide tree. We test validity of our method on 7 sequences of IRE through comparing with Clutal W and T-Coffee.
Key words：Multiple sequences alignment; Secondary structures; Base pairing probability; Structural information vector; Dynamic programming algorithm
1 引 言
多序列比对是生物序列分析的基础，传统的多序列比对（如ClustralW［1］、T-Coffee［2］）通常用于数据库搜索或是结构特点探测，但是对RNA分子，这些方法就不适用了，因为RNA分子的功能主要由其二级结构确定，在进化过程中RNA的结构比序列具有更强的保守性，许多RNA有关的分析研究也正是应用了这一特点，如RNA结构分析［3-5］、RNA同源搜索［6］、非编码RNA探测［7-8］和基于RNA的系统进化推断［9］。而这些RNA序列分析方法都是要求先进行准确的多序列比对，这里的准确，就是指序列比对不仅要考虑序列信息，而且要更多的考虑结构信息。
基于序列和结构信息的RNA多序列比对一般可以分为两类［10］：概率方法和非概率方法。概率方法基于上下文无关语法（SCFG），要求一个初始比对作为输入，而输出的质量对初始比对的依靠性较强。该方法被用于对RNA家族进行建模或是通过比较分析来预测二级结构，比如Cove［11］、RNACAD［12］和Pfold［4］。非概率方法，像MARNA［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］就是采用的这种算法。因为对能量规则进行了简化，对多分支环的处理是用单链碱基的自由能来模拟碱基堆积间的能量增量，所以这种算法可能会模糊两个候选结构间的差别。Mathews提出了一种动态规划算法［22］，来计算RNA二级结构的划分函数，该算法用了最新的最近邻自由能参数［23］，它充分考虑了完整的RNA二级结构的热力学参数，本研究采用该算法来计算划分函数和碱基对配对概率矩阵。
为计算划分函数，需要6个N×N矩阵和两个矢量：
V(i，j)，W(i，j)，WL(i，j)，WMB(i，j)，WMBL(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∈S2Sδ〔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，l+1+Ψ(S1)ij+Ψ(S2)kl+τ(S1(i)，S1(j);S2(k)，S2(l))(6)
初始条件为：Mi，j;k，l=｜(j-i)-(l-k)｜γ，用H表示发卡环的最小尺寸（一般为3），则有条件：j-i≤H+1或l-k≤H+1。上式的前两项分别对两序列中的空位进行计数，第三项表示不配对情况，最后一项表示将两序列分成两段进行比对。该递归空间复杂度为O(n4)，时间复杂度为O(n6)。为了降低计算的复杂度，可以限制两个碱基对(i，j)和(k，l)间的跨度Δ=｜(j-i)-(l-k)｜，对所有部分比对都用该限制。这时算法的计算复杂降为：空间复杂度为O(n3)，时间复杂度为O(n4)。计算完最优比对分值后，采用标准的回溯算法就可以将两条序列对齐。
2.3 结构信息矢量比较
为了进一步降低计算的复杂度，利用Bonhoeffer［24］提出的上、下游概率的概念，我们构造了结构信息矢量。位置 的上游概率的定义为：p＜(i)=∑j＞iPij，下游概率是p＞(i)=∑j＜iPij，根据这个定义，对序列S1可以构造上、下游结构信息矢量p＜s1、p＞s1，它们分别包括了序列S1中所有碱基的上、下游概率。同样可以定义序列Ｓ２的上、下游结构信息矢量p＜s2、p＞s2。这样可能会丢失某些具体的配对信息，但上、下游矢量仍能很好地捕捉到足够的结构信息。这样，目标函数（5）式就可以变为：
∑(ij;kl)∈S(p＜s1(i)·p＜s2(k)+p＞s1(j)·p＞s2(l)+μ)+Nopen·wopen+Next·wext+∑i∈S1，k∈S2Sδ(S1(i)，S2(k))(7)
式中，μ＜0表示对两个在结构水平碱基错配（他们中至少有一个没有构成碱基对）的罚分。剩下的部分是序列相似性分值和空位罚分，这里采用了仿射空位罚分。wopen＜0是一个新开空位罚分，Nopen是新开空位的数量。wext＜0空位扩展罚分，Next表示空位扩展的数量。
同样可以用动态规划算法来得到最好的比对分值：
Mi，j;k，l=MAX(Li，j;k，l，Ri，j;k，l，Pi，j;k，l)(8)
碱基i，j和k，l构成碱基对的情况下：
Pi，j;k，l=Mi，j-1;k，l-1+p＜s1(i)p＜s2(k)+p＞s1(j)p＞s2(l)+∑i∈S1，k∈S2Sδ(S1(i)，S2(k))(9)
碱基i，j和k，l不构成碱基对的情况下：
Pi，j;k，l=Mi，j-1;l-1+μ+∑i∈S1，k∈S2Sδ(S1(i)，S2(k))
Li，j;k，l=MAX(Li，j;k，l-1，Mi，j;k，l-1-wopen)-wext
Ri，j;k，l=MAX(Ri，j-1;k，l，Mi，j-1;k，l-wopen)-wext(10)
Pi，j;k，l表示序列S1(i，j)和序列S2(k，l)对齐的情况，Li，j;k，l表示S1(i，j)和S2(k，l)左边对齐，Ri，j;k，l表示S1(i，j)和S2(k，l)右边对齐。
此时算法的时间复杂度为O(n2)，空间复杂度为O(n3)。
2.4 多序列比对
为计算多序列比对，我们先要计算双序列比对的一致碱基配对概率矩阵，定义如下［14］：
P(S1S2)pq=
P(S1(i，j))p，qP(S2(k，l))p，q for matches
0 for gaps(11)
式中，P(S1(i，j))p，q表示，比对前（比对后的下标为p，q）碱基i，j的配对概率矩阵，P(S2(k，l))p，q的含义相同。对一致碱基配对概率矩阵同样可以构造其结构信息矢量p＜S、p＞S，S表示序列S1和序列S2的共同结构。
对于N条序列的比对，可以用渐进比对的办法完成。首先对N条序列计算所有的两两比对，并根据相似性分值，用权重对组聚类法产生一个指导树，然后沿指导树渐进地对齐N条序列。首先比对两个最高分值的序列，计算出他们的一致碱基配对概率矩阵，构造结构信息矢量，然后由指导树选出下一个序列，用其结构信息矢量和前面得到的一致结构信息矢量比较，再产生一个一致结构信息矢量，一直持续下去，直到比完所有的序列。
3 验证
为验证算法的有效性，我们在Rfam8.0数据库［25］上进行了测试。作为参考，我们用了另两种不同的比对工具：MARNA［10］和RNAlign［13］。我们从IRE(Iron response element)（Rfam［25］的索引号为RF00037）中随机选取了7条序列，IRE是一个短的、保守的茎环结构（见图1（a）中红色区域），它由IRPs（iron response proteins）约束。可以在不同的mRNA的UTRs中发现IRE，这些mRNA一般都参与到铁代谢过程。Rfam库中提供了手工比对结果和其二级结构，我们就以它为标准结果。
我们采用两个独立但互补的分值来评估比对的质量，一个是广泛采用的SPS分值（sum-of-pairs score）［26］，另一个是结构保守指数（structure conservation index，SCI）［27］。SPS定义为所有既在预测比对也在参考比对中出现的字符对的比例，值域为［0，1］，预测比对的质量越好，分值越接近1。与SPS不同，SCI分值不涉及参考比对，它给出了比对中所包含的保守的二级结构信息，由RNAalifold程序计算的分值推出。如果RNAalifold没有在比对中识别出共同的RNA结构，则SCI为0，反之，如果有一组非常保守的结构，则SCI接近1。SCI的值可以大于1，此时暗示着保守的二级结构是由互补突变或是一致突变来维系的。表1给出各个方法的SCI和SPS分值。表1 各个方法的SPS和SCI分值从表1可以看出，我们方法的比对质量不但在序列水平上（SPS分值），而且在结构保守性上（SCI值）要优于MARNA和RNAlign两个程序。
多序列比对的目的就是为了进一步的RNA结构分析，为比较不同方法的差异，我们将不同方法得到的比对结果作为输入，采用比较分析法，Mifold［28］预测了相应的二级结构。图1显示了用各种比对工具得到的不同结果，以及以此比对结果为输入的结构预测结果。通过比较可知，MARNA（图1（b））和RNAlign（图1（c））只是显示了序列间碱基的相似性，而没有考虑保守的茎环信息（图1（a）中红色区域），我们的方法（图1（d））就很好地探测到这些信息。图2给出了由此分值构造的多序列比对指导树。
(a)Manual multiple sequence alignment
AY120878 GGUCGC-GUCAACAGUGUUUGAU-CGAAC
AF266195 AUUCUUGCUUCAACAGUGUUUGAACGGAAU
D8662 GUUCUUGUUUCAACAGUGAUUGAACGGAAC
M16343 GUUCCUGCGUCAACAGUGCUUGGACGGAAC
M58040 UAUAUC-GGAGACAGUGACCUCC-AUAUG
BC001188 AUUAUC-GGGAACAGUGUUUCCC-AUAAU
AF086 UCGUUC-GUCCUCAGUGCAGGGC-AACAG
(b)MARNA
AY120878 -GGUCGCGUCAACAGUGUUUGAUCGAAC-
AF266195 AUUCUUGCUUCAACAGUGUUUGAACGGAAU
D8662 GUUCUUGUUUCAACAGUGAUUGAACGGAAC
M16343 GUUCCUGCGUCAACAGUGCUUGGACGGAAC
M58040 -UAUAUCGGAGACAGUGACCUCCAUAUG-
BC001188 -AUUAUCGGGAACAGUGUUUCCCAUAAU-
AF086786 -UCGUUCGUCCUCAGUGCAGGGCAACAG-
(c)RNAlign
AY120878 GGUC-GCGUCAACAGUGUUUG-AUC-GAAC
AF266195 AUUCUUGCUUCAACAGUGUUUG-AACGGAAU
D8662 GUUCUUGUUUCAACAGUGAUUG-AACGGAAC
M16343 GUUCCUGCGUCAACAGUGCUUG-GACGGAAC
M58040 -UAUAUCGGAGACAGUGACCU-C-CAUAUG
BC001188 -AUUAUCGGGAACAGUGUUUC-C-CAUAAU
AF086786 -UCGUUCGUCCUCAGUGCAGGGCAACAG-
(d)Our method
AY120878 GGUCGCGUC-AACAGUGC-UUGAUCGAAC
AF266195 AUUCUUGCUUCAACAGUGUUUGAACGGAAU
D8662 GUUCUUGUUUCAACAGUGAUUGAACGGGAAC
M16343 GUUCCUGCGUCAACAGUGCUUGGACGGAAC
M58040 UAUAUCGGA-GACAGUGA-CCUCCAUAUG
BC001188 AUUAUCGGG-AACAGUGU-UUCCCAUAAU
AF086786 UCGUUCGUC-CUCAGUGC-AGGGCAACAG
图1 使用不同比对工具得到的比对结果以及由Mifold预测的相应
二级结构
Fig 1 The result of using various alignment tool， and corresponding
predicted RNA structure by Mifold
图2 采用权重两两聚类法得到的、用于多序列比对的指导树
Fig 2 Guide tree using the weighted pair-group clustering methods
for multi-alignment
4 结论
为了更好地分析RNA二级结构，我们提出了基于能量参数的、考虑结构信息的多序列比方法，本方法使用了充分考虑RNA二级结构信息的Mathews动态规划算法来计算划分函数和碱基配对概率矩阵，又通过构造结构信息矢量，在保证有足够的结构信息条件下，减少了计算复杂度。试验证明了我们的方法的有效性，这为今后进行RNA二级结构分析奠定了基础。
参考文献
［1］Thompson JD， Higgins DG， Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting， position-specific gap penalties and weight matrix choice［J］. Nucleic Acids Res，1994，22(22):4673-4680.
［2］Notredame C， Higgins DG， Heringa J. TCoffee: A novel method for fast and accurate multiple sequence alignment［J］. Journal of molecular biology，2000 ，302(1):205-217.
［3］Bindewald E， Shapiro BA. RNA secondary structure prediction from sequence alignments using a network of k-nearest neighbor classifiers［J］. RNA，2006，12(3):342-352.
［4］Knudsen B， Hein J. Pfold. RNA secondary structure prediction using stochastic context-free grammars［J］. Nucleic acids research，2003，31(13):3423-3428.
［5］Pedersen JS， Meyer IM， Forsberg R，et al.J. A comparative method for finding and folding RNA secondary structures within protein-coding regions［J］. Nucleic acids research，2004，32(16):4925-4936.
［6］Griffiths-Jones S， Baterman A， Marshall M，et al.Rfam: an RNA family database［J］. Nucleic Acids Res，2003，31，439-441.
［7］Rivas E， Eddy S. Noncoding RNA gene detection using comparative sequence analysis［J］. BMC Bioinformatics，2001，2：8.
［8］Washietl S， Hofacker I，Stadler P. Fast and reliable prediction of noncoding［J］.RNAs Proc Natl Acad Sci USA，2005，102：2454-2459.
［9］Hudelot C， Gowri-Shankar V， Jow H，et al.RNA-based phylogenetic methods: application to mammalian mitochondrial RNA sequences［J］. Mol Phylogenet Evol，2003，28:241-252.
［10］Siebert S， Backofen R. MARNA: multiple alignment and consensus structure prediction of RNAs based on sequence structure comparisons［J］. Bioinformatics，2005，(16):3352-3359.
［11］Eddy SR， Durbin R. RNA sequence analysis using covariance models［J］. Nucleic Acids Res，1994，22(11):2079-2088.
［12］Brown MPS. RNA modeling using stochastic context-free grammars［M］.Santa Cruz: University of California，1999.
［13］Corpet F， Michot B. RNAlign program: alignment of RNA sequences using both primary and secondary structures［J］. Comput Appl Biosci，1994，10(4):389-399.
［14］Hofacker IL， Bernhart SH， Stadler PF. Alignment of RNA base pairing probability matrices［J］. Bioinformatics，2004，20(14):2222-2227.
［15］Sankoff D. Simultaneous solution of the RNA folding，alignment， and proto-sequence problems［J］. SIAM Journal on Applied Mathematics，1985，45:810-825.
［16］Gorodkin J， Heyer LJ， Stormo GD. Finding common sequence and structure motifs in a set of RNA sequences［A］. Proceedings/International Conference on Intelligent Systems for Molecular Biology［C］.1997.5:120-123.
［17］Mathews DH， Turner DH. Dynalign: an algorithm for finding the secondary structure common to two RNA sequences［J］. Journal of molecular biology，2002，317(2):191-203.
［18］McCaskill JS. The equilibrium partition function and base pair binding probabilities for RNA secondary structure［J］. Biopolymers，1990，29(6-7):1105-1119.
［19］Wuchty S， Fontana W， Hofacker IL， Schuster P. Complete suboptimal folding of RNA and the stability of secondary structures［J］. Biopolymers，1999，49(2):145-165.
［20］Hofacker IL， Fontana W， Stadler PF，et al.Fast folding and comparison of RNA secondary structures［J］. Monatsh Chemie，1994，125:167-188.
［21］Hofacker IL. Vienna RNA secondary structure server［J］. Nucleic Acids Res，2003，31(13):3429-3431.
［22］Mathews DH. Using an RNA secondary structure partition function to determine confidence in base pairs predicted by free energy minimization［J］. RNA，2004，10(8):1178-1790.
［23］Mathews DH， Disney MD， Childs JL，et al.Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure［A］. Proceedings of the National Academy of Sciences of the United States of America［C］. 2004.101(19):7287-7292.
［24］Bonhoeffer S， McCaskill JS， Stadler PF，et al. RNA multi-structurelandscapes. A study based on temperature dependent partition functions［J］. Eur Biophys J，1993，22(1):13-24.
［25］Griffiths-Jones S， Bateman A， Marshall M，et al.Rfam: an RNA family database［J］. Nucleic Acids Res，2003，31(1):439-441.
［26］Thompson JD， Plewniak F， Poch O. A comprehensive comparison of multiple sequence alignment programs［J］. Nucleic acids research，1999，27(13):2682-2690.
［27］Washietl S， Hofacker IL， Stadler PF. Fast and reliable prediction of noncoding RNAs［A］. Proceedings of the National Academy of Sciences of the United States of America［C］. 2005.102(7):2454-2459.
［28］Freyhult E， Moulton V， Gardner P. Predicting RNA structure using mutual information［J］. Applied bioinformatics，2005，4(1):53-59.