FIT2006(第5回情報科学技術フォーラム)A-007

Approximation Algorithms for Optimal RNA SecondaryStructures Common to Multiple Sequences
Takeyuki Tamura∗
1

Tatsuya Akutsu∗

Introduction

RNA secondary structure prediction is an important problem in computational biology and thus many computational
studies have been done. This is a problem of, given an RNA sequence of length n, finding its correct secondarystructure. Usually, RNA secondary structure prediction is modeled as a free energy minimization problem [7, 9]. For
this problem, [10] and [11] proposed simple DP (dynamic programming) algorithms. The time complexities of thoseDP algorithms were O(n3 ) if we ignore the destabilizing energy due to loop regions, otherwise they were at least O(n4 ).
In a basic and simplest version, free-energy minimization of an RNA secondary structure is defined as a problem ofmaximizing the number of complementary base pairs, which is denoted by RNA0 in this paper. Even for RNA0 , onlyan O(n3 ) time simple DP algorithm had been known [7, 9].
An O(n2.776 + (1/²)O(1) ) time approximation algorithm was also shown for RNA0 in [1], which always outputs anRNA secondary structure with the score at least 1 − ² of the maximum, where ² is any positive constant number andthe score denotes the number of complementary base pairs in RNA0 . Although this algorithm can be considered as a
PTAS (polynomial time approximation scheme), it is different from usual PTAS since the problem is not NP-hard but
belongs to P. This algorithm is a combination of an approximation algorithm Aapprox and an exact algorithm Aexact ,where Aapprox is obtained by modifying the original O(n3 ) time DP algorithm for RNA0 , and Aexact is obtained by
combining Valiant’s algorithm with fast funny matrix multiplication.
In order to improve the prediction accuracy, an approach using multiple RNA sequences from the same RNA family
was proposed [6]. An O(n6 ) time exact algorithm was shown in [6] which can optimize structure and alignments when
two RNA sequences are given. Though some efforts have been done [2, 3], the worst case time complexity has not been
improved. In this paper, we show an O(n5 ) time approximation algorithm for optimizing structure and alignments of
two RNA sequences with assuming that the optimal number of base-pairs is more than O(n0.75 ). We also show thatthe problem to optimize structure and alignments for given N sequences is NP-hard and introduce a constant-factor
approximation algorithm.

2

RNA secondary structure when multiple sequences are given

Let A1 = a1,1 a1,2 . . . a1,n1 , A2 = a2,1 a2,2 . . . a2,n2 , . . ., and AN = aN,1 aN,2 . . . aN,nN be RNA sequences, wheremax{n1 , n2 , . . . , nN } = n. Thus, A1 , A2 , · · ·, AN are strings over an alphabet Σ = {a, u, g, c}. A family of pairsof indices MN = h{(1i , 1j )|1 ≤ 1i < 1j ≤ n1 , (a1,1i , a1,1j ) is a base pair}, {(2i , 2j )|1 ≤ 2i < 2j ≤ n2 , (a2,2i , a2,2j ) is abase pair} ,. . ., {(Ni , Nj )|1 ≤ Ni < Nj ≤ nN , (aN,Ni , aN,Nj ) is a base pair}i is called an N-common RNA secondary
structure if a1,1i = a2,2i = · · · = aN,Ni and a1,1j = a2,2j = · · · = aN,Nj , and no distinct pairs (xi , xj ), (xh , xk ) in MN
satisfy xi ≤ xh ≤ xj ≤ xk for all x (1 ≤ x ≤ N ). The score of MN is defined as the number of base pairs in each elementof MN (i.e., |e| for any e in MN ), and denoted by score(MN ). Then, RNA0 (N ) is defined as follows: given N RNA
sequence A1 , A2 , . . . , AN , to find an N -common RNA secondary structure M with the maximum score. In RNA0 (N ),such a structure is also called an optimal N-common RNA secondary structure, and denoted by OP T (RNA0 (N )).

3

1 − ² approximation algorithm for RNA0 (2)As mentioned above, in RNA0 (2), two sequences are given. Let (i1 , j1 ) be a pair of indices which correspond to
the leftmost and rightmost residues of the first sequence respectively. Similarly, let (i2 , j2 ) be a pair of indices whichcorrespond to the leftmost and rightmost residues of the other sequence respectively. RNA0 (2) can be solved in O(n6 )time by the following DP procedure [6]:
D(i1 + 1, j1 , i2 , j2 ) D(i1 , j1 − 1, i2 , j2 )D(i1 , j1 , i2 + 1, j2 )D(i1 , j1 , i2 , j2 ) = max
D(i1 , j1 , i2 , j2 − 1) D(i1 + 1, j1 − 1, i2 + 1, j2 − 1) + f (i1 , j1 , i2 , j2 )maxi1 <k1 <j1 ,i2 <k2 <j2 {D(i1 , k1 , i2 , k2 ) + D(k1 + 1, j1 , k2 + 1, j2 )},where f (a, u, a, u) = 1, f (u, a, u, a) = 1, f (g, c, g, c) = 1, f (c, g, c, g) = 1, otherwise f is zero.
The technique for the approximation algorithm of RNA0 [1] can be applied to RNA0 (2) with assuming OP T (RNA0 (N ))is large, where additional ideas are required for analysis of the approximation ratio. As in [1], we do not compute maxi1 <k1 <j1 ,i2 <k2 <j2 {D(i1 , k1 , i2 , k2 ) + D(k1 + 1, j1 , k2 + 1, j2 )} exactly. Instead, we compute the maximum ofD(i1 , k1 , i2 , k2 ) + D(k1 + 1, j1 , k2 + 1, j2 )} for O({nα + n1−β }2 ) values of (k1 , k2 ), where α and β (0 < α, β < 1) areappropriate constants to be determined later. We define a sequence of indices fi+1 (h) and fj+1 (h) for h1 = 0, 1, 2, . . . byfi+1 (0) = i1 + dnα e,fj−1 (0) = j1 − dnα e+++βfi1 (h1 + 1) = fi1 (h1 ) + d(fi1 (h1 ) − i1 ) e, fj−1 (h1 + 1) = fj−1 (h1 ) − d(j1 − fj−1 (h1 ))β e
Next, we define I1 (i, j) by
I1 (i1 , j1 ) ={k1 |i1 < k1 ≤ nαor j1 − nα ≤ k1 ≤ j1 } ∪ {fi+1 (h1 )|fi+1 (h1 ) ≤ (i1 + j1 )/2}∪{fj−1 (h1 )|fj−1 (h1 ) ≥ (i1 + j1 )/2}.

∗ Bioinformatics Center, Institute for Chemical Research, Kyoto University, Uji, Kyoto 611-0011, Japan.

15

FIT2006(第5回情報科学技術フォーラム)Similarly, we define a sequence of indices fi+2 (h2 ) and fj+2 (h2 ) for h2 = 0, 1, 2, . . . and I1 (i, j) by
fi+2 (0) = i2 + dnα e,fj−2 (0) = j2 − dnα e+++β
fi2 (h2 + 1) = fi2 (h2 ) + d(fi2 (h2 ) − i2 ) e, fj−2 (h2 + 1) = fj−2 (h2 ) − d(j2 − fj−2 (h2 ))β e
I2 (i2 , j2 ) ={k2 |i2 < k2 ≤ nα

or j2 − nα ≤ k2 ≤ j2 } ∪ {fi+2 (h2 )|fi+2 (h2 ) ≤ (i2 + j2 )/2}∪{fj−2 (h2 )|fj−2 (h2 ) ≥ (i2 + j2 )/2}.

Then, the approximation algorithm
 0 Aapprox (2) is expressed by the following DP procedure:
D (i1 + 1, j1 , i2 , j2 )

D0 (i , j − 1, i2 , j2 )
 0 1 1D (i1 , j1 , i2 + 1, j2 )D0 (i1 , j1 , i2 , j2 ) = max
D0 (i1 , j1 , i2 , j2 − 1)0+ f (i1 , j1 , i2 , j2 )
 D (i1 + 1, j1 − 1, i2 + 1, j2 − 1)
maxk1 ∈I1 (i1 ,j1 ),k2 ∈I2 (i1 ,j1 ) {D0 (i1 , k1 , i2 , k2 ) + D0 (k1 + 1, j1 , k2 + 1, j2 )},where f (a, u, a, u) = 1, f (u, a, u, a) = 1, f (g, c, g, c) = 1, f (c, g, c, g) = 1, otherwise f is zero.
Lemma 1 Aapprox (2) works in O(n2α+4 + n6−2β + n5+α−β ) time.
Here, we define the error of an N-common secondary structure MN to OP T (RNA0 (N )) to be score(RNA0 (N ))) −score(MN ) (note that this value must be non-negative).
Lemma 2 The error of a secondary structure MN computed by Aapprox is O(n1+αβ−α ).
Theorem 1 When OP T (RNA0 (N )) > O(n0.75 ), an N-common RNA secondary structure with the score at least 1 − ²of the maximum can be computed in O(n5 ) time, where ² is any positive constant number.
It is known that the Longest Common Subsequence problem (LCS) over an alphabet of size 2 is NP-hard [4, 5]. Usinga reduction from LCS, we have:
Theorem 2 RNA0 (N ) is NP-hard if N is not fixed.
Theorem 3 There is an O(N n) time approximation algorithm for RNA0 (N ) with the score at least 1/4 of themaximum.

4Concluding remarks

In this paper, we proposed an O(n5 ) time (1-²)-approximation algorithm for optimal RNA secondary structures common
to two sequences with assuming that the optimal score is more than O(n0.75 ). In order to delete this assumption, weshould combine Valiant’s algorithm [8] with the proposed algorithm. We also showed that the problem is NP-hard forgeneral N and introduced an O(N n) time 1/4-approximation algorithm. Improvement of this approximation ratio isleft as an open problem.

References[1] T. Akutsu, Approximation and exact algorithms for RNA secondary structure prediction and recognition ofstochastic context-free languages, Journal of Combinatorial Optimization (1999) 3:321–336[2] V. Bafna, S. Muthukrishnan, and R. Ravi, Computing similarity between RNA strings, Proc. 6th Symp.
Combinatorial Pattern Matching (1995) 1–16
[3] V. Bafna, H. Tang, and S. Zhang, Consensus folding of unaligned RNA sequences revisited, Journal ofComputational Biology (2006) 13(2):283–295[4] D. Maier , The complexity of some problems on subsequences and supersequences, Journal of the ACM (1978)
25(2):322–336[5] M. Middendorf, On finding various minimal, maximal, and consistent sequences over a binary alphabet, Theoretical Computer Science (1995) 145:317–327[6] D. Sankoff, Simultaneous solution of the RNA folding, alignment and protosequence problems, SIAM Journal
on Applied Mathematics (1985) 45(5):810–825[7] J. Setubal and J. Meidanis, Introduction to Computational Molecular Biology, PWS Publishing Company(1997)[8] L. G. Valiant, General context-free recognition in less than cubic time, Journal of Computer and System Science(1975) 10:308–315[9] M. S. Waterman, Introduction to computational biology, Chapman and Hall (1995)[10] M. S. Waterman and T. F. Smith, RNA secondary structure: A complete mathematical analysis, Mathematical
Biosciences (1978) 41:257–266[11] M. Zuker and P. Stiegler, Optimal computer folding of larger RNA sequences using thermodynamics andauxiliary information, Nucleic Acids Research (1981) 9:133–148

16