- Proceedings
- Open Access
Local structural alignment of RNA with affine gap model
- Thomas King-Fung Wong^{1}Email author,
- Brenda Wing-Yan Cheung^{1},
- Tak-Wah Lam^{1} and
- Siu-Ming Yiu^{1}Email author
https://doi.org/10.1186/1753-6561-5-S2-S2
© King-Fung Wong et al; licensee BioMed Central Ltd. 2011
- Published: 28 April 2011
Abstract
Background
Predicting new non-coding RNAs (ncRNAs) of a family can be done by aligning the potential candidate with a member of the family with known sequence and secondary structure. Existing tools either only consider the sequence similarity or cannot handle local alignment with gaps.
Results
In this paper, we consider the problem of finding the optimal local structural alignment between a query RNA sequence (with known secondary structure) and a target sequence (with unknown secondary structure) with the affine gap penalty model. We provide the algorithm to solve the problem.
Conclusions
Based on an experiment, we show that there are ncRNA families in which considering local structural alignment with gap penalty model can identify real hits more effectively than using global alignment or local alignment without gap penalty model.
Keywords
- Query Sequence
- Structural Alignment
- Alignment Score
- Global Alignment
- ncRNA Family
Background
A non-coding RNA (ncRNA) is a RNA molecule that does not translate into proteins. It has been shown to be involved in many biological processes [1–4]. The number of ncRNAs within the human genome was underestimated before, but recently some databases reveal over 212,000 ncRNAs [5] and more than 1,300 ncRNA families [6]. Large discoveries of ncRNAs and their families show the possibilities that ncRNAs may be as diverse as protein molecules [7]. Identifying ncRNAs is an important problem in biological study. However, it is time consuming and there is no effective method to identify ncRNAs in a laboratory, predicting ncRNAs based on known ncRNAs using comparative computational approach is one of the promising directions to identify potential candidates for further verification.
Most of the computational approaches are based on the observation that if two different ncRNA molecules are in the same family (with similar biological functions), they usually exhibit similar sequences as well as secondary structures. One common approach [8–10] is as follows. We pick an ncRNA member of a family with known sequence and secondary structure (referred as the query), scan along a genomic sequence and for each possible region (referred as the target), perform an alignment between the query and the target to obtain a similarity measure to decide if the region is a potential ncRNA candidate for that family. The similarity measure may only base on the sequence or both the sequence and secondary structure (the latter case is referred as structural alignment). Along this direction, there are some approaches [11–14] that make use of secondary structure prediction tools to predict the secondary structure to be formed by the target assuming that it is an ncRNA before performing the alignment. The accuracy may, however, depend on the accuracy of the secondary structure prediction tools.
Instead of using one member of a family, some other approaches [15] use a set of ncRNAs from the same family to train a model (e.g. covariance model). Then, using this model to scan a genomic sequence to identify potential regions that are ncRNA candidates of that family. What information (sequence similarity and/or secondary structure) to be captured from the known ncRNAs depends on how we define the model. However, in some cases, we may not have enough known members in a family to train a model. In this paper, we focus on the problem that uses one known member as the query and align it with a target sequence. We remark that there are also other computational methods that identify ncRNAs without using known members in a family. For example, some try to identify ncRNAs by considering the stability of secondary structures formed by the substrings of a given genome [16]. This method may not be very effective because a random sequence with high GC composition also allows an energetically favorable secondary structure [17]. So, the comparative approach we described in the above is still one of the most popular approaches.
We consider the following problem. Given a query sequence together with its secondary structure, we try to identify the substring in the given target sequence (with unknown secondary structure) that can align to a substring in the query sequence with the highest structural similarity score based on the affine gap model (see next section for formal definitions). We assume that the secondary structures of the ncRNAs are regular, that is, they do not have pseudoknots (no two base pairs crossing each other). This type of ncRNAs is found to be the most abundant in existing databases. We consider all possible substrings of the query sequence, even for those substrings that cover only one of the end points of some base pairs in the structure.
Our result
We propose a local structural alignment algorithm with affine gap model which assumes the secondary structure of the query is known while that of the target sequence is unknown. The time complexity of our algorithm is O(mn^{3}) which is the same as the best algorithm for global alignment for this problem where m, n are the lengths of the query and the target, respectively. We evaluated our algorithm using real data from Rfam database. According to the preliminary experiment, it shows that there are ncRNA families in which considering local structural alignment algorithm with affine gap model can distinguish real members from false hits more effectively than using global alignment or local alignment without affine gap model.
Preliminaries
An ncRNA molecule can be regarded as a sequence of four characters {A, C, G, U}, each character is referred as a base. Some of these bases may form pairs (linked up by a hydrogen bond) with some restrictions such as each base can only pair up with at most one other base and only complementary bases can form a pair (e.g. (A, U), (C, G), (U,G)). The set of base pairs formed by the molecule is referred as its secondary structure.
Formally speaking, let S = s_{1}s_{2} … s_{ m } be a length-m ncRNA sequence where s_{ i } ∈ {A, C, G, U} for 1 ≤ i ≤ m and M be the secondary structure of S. M is represented as a set of base pair positions. i.e. M = {(i, j)|1 ≤ i <j ≤ m, (s_{ i }, s_{ j }) is a base pair}. If (s_{ i }, s_{ j }) is base pair, then (s_{ i }, s_{ j }) ∈ {(A, U), (C, G), (G, C), (G, U), (U, A), (U, G)}. Let M_{ x,y } ⊆ M be the set of base pairs within the subsequence s_{ x }s_{ x }_{+1}…s_{ y }, 1 ≤ x <y ≤ m, i.e., M_{ x,y } = {(i, j) ∈ M|x ≤ i <j ≤ y}. Note that if (i, j) ∈ M and only i or j inside the region [x…y], then (i, j) ∉ M_{ x,y }. We assume that there is no two base pairs sharing the same position, i.e., for any (i_{1}, j_{1}), (i_{2}, j_{2}) ∈ M, i_{1} ≠ j_{2}, i_{2} ≠ j_{1}, and i_{1} = i_{2} if and only if j_{1} = j_{2}.
A regular structure is the structure in which there does not exist any two base pairs crossing each other. The formal definition is as follows:
Definition 1 M_{ x,y } is a regular structure if there does not exist two base pairs (i, j), (k, l) ∈ M_{ x,y } such that i <k <j <l or k <i <l <j.
Note that an empty set is also considered as a regular structure.
Problem definition
Structural alignment with affine gap model
where η(i) is the corresponding position in S according to the position i in S′; γ(u_{1},u_{2}) and δ(u_{1}, u_{2}, v_{1}, v_{2}) where u_{1}, u_{2}, v_{1}, v_{2} ∈ {A, C, G, U}, are scores for character similarity and for base pair similarity respectively; k and l is the number of gaps and the total length of all gaps; h and s is the gap starting and extending penalty.
Definition 2 An optimal global structural alignment between S and T is a structural alignment of S and T such that the alignment score is maximum.
Let S[x…y] where 1 ≤ x, y ≤ m be a substring of S with secondary structure M_{ x,y } (where S[x…y] is an empty string with empty structure if x >y). Similarly, let T[x′…y′] where 1 ≤ x′, y′ ≤ n be a substring of T (where T[x′…y′] is an empty string if x′ >y′).
Definition 3 An optimal local structural alignment between S and T is a global structural alignment between two sub stings of S and T, S[x…y] and T[x′…y′] where 1 ≤ x, y ≤ m and 1 ≤ x′, y′ ≤ n of S and T such that the alignment score between them is maximum over all possible substrings.
Given S (with known secondary structure) and T (with unknown structure), we want to compute an optimal local structural alignment with affine gap penalty between S and T.
Results and discussion
The details of the algorithm for solving the problem will be given in Method Section. In this section, we evaluate the resulting algorithm and show that considering local structural alignment with affine gap model can improve the effectiveness of locating ncRNAs for the families in which members may have variable size of hairpins, loops or stems when compared to using global alignment [10], local alignment without gap penalty model and Gotohscan [8]. Note that the differences in size of hairpins, loops or stems represent gaps in the corresponding sequences.
To test the algorithm, we selected around twenty ncRNA families in which the members have variable sizes of hairpins, loops or stems. We construct our testing cases based on real ncRNAs as follows. For each family, we first select a seed member (i.e. In Rfam database, there is a set of reliable members which are regarded as seed members) as the query sequence Q. To demonstrate the effectiveness of the affine gap model, we select the longest seed member as this query sequence. We then created a long random sequence with even distribution of four characters {A, C, G, T} to simulate a long genome. The length of this long random sequence is around ten times of the total length of all the seed members of the family. Finally, we embedded the full members (i.e. all the members including the seed members) of the family (except the one chosen as query) into this long random sequence in arbitrary positions. If there are more than 100 members, then we randomly picked 100 of them. This resulting sequence is our T.
The details of the ncRNA families used in the experiments.
Family | Query Sequence ID | Length | Number of members embedded |
---|---|---|---|
RF00014 | CP000468.1/2032552-2032638 | 87 | 96 |
RF00021 | CP000851.1/113395-113522 | 128 | 100 |
RF00022 | AAND01000021.1/495-707 | 213 | 100 |
RF00027 | AAPE01289140.1/8905-8994 | 90 | 100 |
RF00032 | S49118.1/1081-1106 | 26 | 100 |
RF00033 | Y15844.1/450-543 | 94 | 100 |
RF00034 | BX571867.1/288515-288628 | 114 | 100 |
RF00038 | AJ132964.1/66-198 | 133 | 100 |
RF00039 | AF370716.1/3603-3656 | 54 | 100 |
RF00042 | X55895.1/474-565 | 92 | 100 |
RF00043 | Z47410.1/1220-1294 | 75 | 21 |
RF00044 | M11813.1/4883-5126 | 244 | 8 |
RF00046 | AY013245.2/62208-62303 | 96 | 76 |
RF00048 | AF504534.1/666-726 | 61 | 100 |
RF00386 | AF363455.1/1-122 | 122 | 100 |
RF00643 | AASG02000279.1/67999-67862 | 138 | 100 |
RF00661 | AC154049.1/4734-4855 | 122 | 100 |
RF01051 | AE014299.1/1112481-1112574 | 94 | 100 |
Summary of comparison on results between global alignment, local alignment without gap penalty and local alignment with affine gap penalty when using the smallest threshold such that there is no false positive.
Family | Number of members | Number of misses | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|
Gotohscan [8] | % | Global [10] | % | Local | % | Local with affine gap | % | |||||
RF00014 | 96 | 2 | 2.1% | 0 | 0% | 0 | 0% | 0 | 0% | |||
RF00021 | 100 | 10 | 10% | 5 | 5% | 5 | 5% | 2 | 2% | |||
RF00022 | 100 | 59 | 59% | 20 | 20% | 19 | 19% | 4 | 4% | |||
RF00027 | 100 | 100 | 100% | 15 | 15% | 9 | 9% | 2 | 2% | |||
RF00032 | 100 | 59 | 59% | 4 | 4% | 1 | 1% | 0 | 0% | |||
RF00033 | 100 | 29 | 29% | 27 | 27% | 27 | 27% | 25 | 25% | |||
RF00034 | 100 | 71 | 71% | 11 | 11% | 22 | 22% | 7 | 7% | |||
RF00038 | 100 | 88 | 88% | 0 | 0% | 0 | 0% | 0 | 0% | |||
RF00039 | 100 | 100 | 100% | 1 | 1% | 1 | 1% | 1 | 1% | |||
RF00042 | 100 | 10 | 10% | 0 | 0% | 0 | 0% | 0 | 0% | |||
RF00043 | 21 | 3 | 14.3% | 0 | 0% | 0 | 0% | 0 | 0% | |||
RF00044 | 8 | 1 | 12.5% | 0 | 0% | 0 | 0% | 0 | 0% | |||
RF00046 | 76 | 9 | 11.8% | 2 | 2.6% | 1 | 1.3% | 0 | 0% | |||
RF00048 | 100 | 17 | 17% | 0 | 0% | 0 | 0% | 0 | 0% | |||
RF00386 | 100 | 88 | 88% | 63 | 63% | 62 | 62% | 6 | 6% | |||
RF00643 | 100 | 98 | 98% | 4 | 4% | 13 | 13% | 0 | 0% | |||
RF00661 | 100 | 100 | 100% | 87 | 87% | 77 | 77% | 30 | 30% | |||
RF01051 | 100 | 100 | 100% | 91 | 91% | 85 | 85% | 52 | 52% | |||
average | 53.9% | 18.4% | 17.9% | 7.2% |
Summary of comparison on results between global alignment, local alignment without gap penalty and local alignment with affine gap penalty when setting the threshold which allows 5% or 10% of false positives.
Family | Number of members | Number of misses | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
False positive rate=5% | False positive rate=10% | |||||||||
Gotohscan | Global | Local | Local with affine gap | Gotohscan | Global | Local | Local with affine gap | |||
RF00014 | 96 | 2 | 0 | 0 | 0 | 2 | 0 | 0 | 0 | |
RF00021 | 100 | 10 | 1 | 1 | 1 | 10 | 1 | 1 | 1 | |
RF00022 | 100 | 51 | 9 | 5 | 2 | 35 | 4 | 4 | 2 | |
RF00027 | 100 | 100 | 3 | 5 | 0 | 100 | 2 | 2 | 0 | |
RF00032 | 100 | 59 | 0 | 0 | 0 | 37 | 0 | 0 | 0 | |
RF00033 | 100 | 27 | 1 | 25 | 24 | 26 | 1 | 1 | 24 | |
RF00034 | 100 | 71 | 1 | 0 | 0 | 71 | 1 | 0 | 0 | |
RF00038 | 100 | 88 | 0 | 0 | 0 | 88 | 0 | 0 | 0 | |
RF00039 | 100 | 100 | 0 | 0 | 0 | 100 | 0 | 0 | 0 | |
RF00042 | 100 | 10 | 0 | 0 | 0 | 10 | 0 | 0 | 0 | |
RF00043 | 21 | 3 | 0 | 0 | 0 | 3 | 0 | 0 | 0 | |
RF00044 | 8 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | |
RF00046 | 76 | 9 | 0 | 0 | 0 | 9 | 0 | 0 | 0 | |
RF00048 | 100 | 11 | 0 | 0 | 0 | 11 | 0 | 0 | 0 | |
RF00386 | 100 | 88 | 58 | 56 | 1 | 88 | 48 | 38 | 1 | |
RF00643 | 100 | 98 | 1 | 4 | 0 | 98 | 0 | 2 | 0 | |
RF00661 | 100 | 100 | 87 | 66 | 23 | 100 | 81 | 52 | 14 | |
RF01051 | 100 | 100 | 79 | 85 | 47 | 100 | 79 | 81 | 39 |
Summary of the area (normalized) under ROC curve for false positive rate ≤ 10%
Family | Area (normalized) under ROC curve | |||
---|---|---|---|---|
Gotohscan | Global | Local | Local with affine gap | |
RF00014 | 0.98 | 1.0 | 1.0 | 1.0 |
RF00021 | 0.9 | 0.99 | 0.99 | 0.99 |
RF00022 | 0.53 | 0.92 | 0.93 | 0.98 |
RF00027 | 0.0 | 0.96 | 0.96 | 1.0 |
RF00032 | 0.61 | 0.99 | 1.0 | 1.0 |
RF00033 | 0.73 | 0.93 | 0.79 | 0.76 |
RF00034 | 0.29 | 0.98 | 0.99 | 0.99 |
RF00038 | 0.12 | 1.0 | 1.0 | 1.0 |
RF00039 | 0.0 | 1.0 | 1.0 | 1.0 |
RF00042 | 0.9 | 1.0 | 1.0 | 1.0 |
RF00043 | 0.86 | 1.0 | 1.0 | 1.0 |
RF00044 | 0.88 | 1.0 | 1.0 | 1.0 |
RF00046 | 0.88 | 1.0 | 1.0 | 1.0 |
RF00048 | 0.89 | 1.0 | 1.0 | 1.0 |
RF00386 | 0.12 | 0.42 | 0.49 | 0.98 |
RF00643 | 0.02 | 0.99 | 0.96 | 1.0 |
RF00661 | 0.0 | 0.14 | 0.36 | 0.79 |
RF01051 | 0.0 | 0.18 | 0.17 | 0.56 |
Our program take around 15 seconds for performing local structural alignment with affine gap model between query and target of around 150 bases long, and around 30 seconds for 200 bases long. We tested the program on a machine with 2.4GHz dual-core CPU and 8G memory.
Conclusions
In the paper, we provided an algorithm to handle local structural alignment with affine gap model of RNA with regular structure that compute the optimal alignment. Our experiments show that the solution is effective for some ncRNA families in which members may have varying sizes on hairpins, loops or stems (contributing to large gaps) when compared to using only global alignment or local alignment without gap model. And also we have not yet studied different types of gap penalty model and the effect of setting different gap penalty parameters. Other interesting directions include speeding up the algorithm and considering other more complicated structures (e.g. the structures with pseudoknots). In the mean time, we have completed the algorithm of computing local structural alignment for simple pseudoknots structure, and we are now in the progress of performing experiments.
Methods
We develop a dynamic programming algorithm to solve the problem. Before we describe the method, we would like to define some variations of alignments which will be used in our algorithm. Let S[1…m] be the query sequence with known structure M and T[1…n] be the target sequence with unknown structure.
Definition 4 Optimal prefix-global structural alignment between S[1…m] and T[1…n] is to find a prefix S[1…y] where 0 ≤ y ≤ m (i.e. S is an empty string when y = 0) such that the score of the optimal global structural alignment between the prefix S[1…y] and T[1…n] is maximum.
Definition 5 Optimal suffix-global structural alignment between S[1…m] and T[1…n] is to find S[x…m] where 1 ≤ x ≤ m + 1 (i.e. S is an empty string when x = m + 1) such that the score of the optimal global structural alignment between the suffix S[x…m] and T[1…n] is maximum.
Definition 6 Optimal semi-global structural alignment between S[1…m] and T[1…n] is to find a substring S[x…y] where 1 ≤ x, y ≤ m such that the score of the optimal global structural alignment between the substring S[x…y] and T[1…n] is maximum.
Let the affine gap model be h + sL, where h is the gap opening penalty, s represents a gap extension penalty, and L denotes the length of gap. Our method consists of two steps. In the first step, we compute the optimal semi-global structural alignment between S and all possible substrings of T. In the second step, we obtain the optimal local structural alignment between S and T resulted in the first step. Define A(p, q, e, f) to be the score of the optimal semi-global structural alignment between S[p…q] and T [e…f]. The score of the optimal local structural alignment between S and T can be obtained from the entry max_{ e }_{≤}_{ f }_{+1}A(1, m, e, f). We first show how to compute A, then show how to use the structure of S to guide the computation of A without considering all possible combinations of p, q.
When considering any substring S′ = S[x′…y′] of S[x…y], there are four possible cases: (1) S′ is equal to S (i.e. x′ = x, y′ = y); (2) S′ is a proper prefix in S (i.e. x′ = x, y′ <y); (3) S′ is a proper suffix in S (i.e. x′ >x, y′ = y); (4) S′ is a substring of S[x + 1…y – 1] (i.e. x′ >x, y′ <y); Therefore, we can consider each case one by one when computing the value of A.
Define A_{1}(p, q, e, f) to be the score of the optimal global structural alignment between S[p…q] and T[e…f]. Define A_{2}(p, q, e, f) to be the score of the optimal prefix-global structural alignment between S[p…q – 1] and T[e…f]. Define A_{3}(p, q, e, f) to be the score of the optimal suffix-global structural alignment between S[p + 1…q] and T[e…f]. Define A_{4}(p, q, e, f) to be the score of the optimal semi-global structural alignment between S[p + 1…q – 1] and T[e…f].
The value of A(p, q, e, f) can be computed recursively and it is the maximum value of four cases: (1) when S′ = S[p, q] (i.e. A_{1}(p, q, e, f)); (2) when S′ is a proper prefix of S[p, q] (i.e. A_{2}(p, q, e, f)); (3) when S′ is a proper suffix of S[p, q] (i.e. A_{3}(p, q, e, f); (4) when S′ is a substring of S[p + 1, q – 1] (i.e. A_{4}(p, q, e, f); Lemma 1 summarizes these cases.
The following subsections describe how to compute A_{1},A_{2},A_{3},A_{4}.
Calculation of A_{1}
When considering the optimal global structural alignment (with affine gap model) between S[p…q] and T [e…f], there are nine possible cases: (1) S[p] is aligned with T[e] and S[q] with T[f]; (2) S[p] with T[e] and S[q] with space;(3) S[p] with T[e] and T[f]withspace; (4) S[p] with space and S[q] with T[f]; (5) S[p] with space and S[q] with space; (6) S[p] with space and T[f] with space; (7) T[e] with space and S[q] with T[f]; (8) T[e] with space and S[q] with space; (9) T[e] with space and T[f] with space. Hence, we can consider each case one by one when computing the value of A_{1}.
Define A_{1}_{ x }(p, q, e, f), where 1 ≤ x ≤ 9, to be the score of the optimal global structural alignment between S[p…q] and T[e…f] where the above case x is satisfied. (i.e. if x = 1, then S[p] is aligned with T[e] and S[q] with T[f]).
The value of A_{1}(p, q, e, f) can be computed recursively and it is the maximum value of nine cases. Lemma 2 summarizes these cases.
We will describe the calculation of A_{12}. Similar skill can be applied for the others (i.e. A_{11}, A_{13}, … , A_{19}).
Calculation of A_{12}
A_{12}(p, q, e, f) is the score of the optimal global structural alignment between S[p…q] and T[e…f], which aligns S[p] with T[e] and S[q] with space. There are three situations and we need to consider them one by one. Note that according to the affine gap model, the penalty of a first space in a gap (i.e. which is h + s) is different from the penalty of the other space in a gap (i.e. which is s). Situation I: when (p, q) is a base pair - aligning the base pair S[p] with T[e] and S[q] with space. Considering the alignment between S[p + 1…q – 1] and T[e + 1…f], if S[q – 1] is aligned with space (i.e. case 2, case 5 and case 8), then a penalty s should be considered. Otherwise (i.e. for the other six cases), a penalty h + s should be considered. Situation II: when ∃q′ where p <q′ <q such that (p, q′) is a base pair - we need to find k ∈ [e – 1, f] such that the sum of the alignment score between S[p, q′] and T[e, k], and that between S[q′ + 1, q] and T[k + 1, f] is maximum. Since S[p] is aligned with T[e] and S[q] with space, the alignment between S[p,q′] and T[e, k] should satisfy the case 1, case 2 and case 3 (i.e. S[p] is aligned with T[e]). Similarly, the alignment between S[q′ + 1, q] and T[k + 1, f] should satisfy the case 2, case 5 and case 8 (i.e. S[q] is aligned with space). Situation III: when p does not form base pair with any base q′ ∈ [p, q] - we align base S[p] with T[e]. Then the alignment between S[p + 1…q] and T[e+ 1…f] should satisfy the case 2, case 5 and case 8 (i.e. S[q] is aligned with space). Lemma 3 summarizes these situations:
Calculation of A_{2}
When considering the optimal prefix-global structural alignment (with affine gap model) between S[p…q] and T[e…f], there are four possible cases: (1) S[p] is aligned with T[e]; (2) S[p] with space; (3) T[f] with space; and (4) an empty string of S with T.
Define A_{2x}(p, q, e, f), where 1 ≤ x ≤ 3, to be the score of the optimal prefix-global structural alignment between S[p…q] and T[e…f] where the above case x is satisfied. (i.e. if x = 1, then S[p] is aligned with T[e]). Note that we do not need to define function for the case 4 because the corresponding score is – h – s(f – e + 1). The value of A_{2}(p, q, e, f) can be computed recursively and it is the maximum value of four cases. Lemma 4 summarizes these cases.
Lemma 4
A_{2}(p, q, e, f) = max{A_{21}[ p, q, e, f], A_{22}[p, q, e, f], A_{23}[p, q, e, f], – h – s(f – e + 1)}
We will describe the calculation of A_{22}. Similar skill can be applied to calculate A_{21} and A_{23}.
Calculation of A_{22}
The following lemma lists out the computation of A_{22}.
A_{22}(p, q, e, f) is the score of the optimal prefix-global structural alignment between S[p…q – 1] and T[e…f], where S[p] is aligned with space. Similar to A_{12}, there are also the same three situations. Situation I: when (p, q) is a base pair - aligning the base pair S[p] with space. Since a prefix of S[p…q – 1] is considered, there are two possibilities: a prefix of S[p + 1…q – 1] is aligned with T[e…f] (i.e. semi-global alignment), or the whole sequence S[p+ 1…q – 1] is aligned with T[e…f] (i.e. global alignment). Situation II: when ∃q′ where p <q′ <q such that (p, q′) is a base pair - we need to find k ∈ [e – 1, f] such that the sum of the alignment score between S[p, q′] and T[e, k], and that between S[q′ + 1, q] and T[k + 1, f] is maximum. Since a prefix of S[p…q – 1] is considered, there are two possibilities: (1) the whole sequence S[p, q′] is aligned with T[e, k] (i.e. global alignment) and a prefix of S[q′ + 1, q] is aligned with T[k + 1, f] (i.e. semi-global); (2) a prefix of S[p, q′] is aligned with T[e, k] (i.e. semi-global) only. Situation III: when p does not form base pair with any base q′ ∈ [p, q] - we align base S[p] with space. For each possibility of situation I & III, there are also two conditions: if S[p + 1] is aligned with T[e] or T[e] is aligned with space, the penalty score h + s should be considered. Otherwise, if S[p + 1] is aligned with space, then the penalty score s should be considered. The lemma 5 summarizes these cases.
The calculations for A_{3} and A_{4} are similar. In the following subsection, we will describe the time complexity of the algorithm.
Time complexity
To fill the dynamic programming table, not all entries for all possible subrange of S needs to be filled. According to the design of the dynamic programming, there are three cases:
Case 1: if (p, q) ∈ M_{ p,q }, then all the entries for S[p, q] of all tables (i.e. A, A_{1}, A_{2}, A_{3}, A_{4}, A_{11}, …, etc.) can be computed from the entries for S[p – 1, q + 1].
Case 2: if ∃q′ <q s.t. (p, q′) ∈ M_{ p,q }, then all the entries for S[p, q] of all tables can be computed from the entries for S[p, q′] and S[q′ + 1, q].
Case 3: if ∄q′ s.t. (p, q′) ∈ M_{ p,q }, then all the entries for S[p, q] of all tables can be computed from the entries for S[p + 1, q].
We only need to fill in the entries for all the tables provided (p, q) can be obtained from (1, m) by applying ζ function repeatedly. Considering the ζ function, each time the total size of the subregions outputted cannot be greater than the size of the input region and each of the subregions outputted is smaller than the input region. Therefore, in total there are only O(m) such (p, q) values. Also, there are O(n^{2}) values of different (e, f) values, and for each entry, it takes O(n) because of the consideration of e – 1 ≤ k ≤ f in the case that ∃q′ <q s.t. (p, q′) ∈ M_{ p,q }. After finishing the calculation of values A(1, m, e, f) for all 1 ≤ e, f ≤ n, the final answer (i.e. max_{e≤f+1}{A(1, m, e, f)}) can be computed in O(n^{2}) time. Therefore the total time complexity = O(mn^{3}) + O(n^{2}) = O(mn^{3}).
Theorem 1 For any sequence S[1..m] with regular structure and any sequence T[1…n] with unknown structure, the optimal local alignment score between S[1..m] and T[1..n] can be computed in O(mn^{3}).
Declarations
Acknowledgements
The project is partially supported by the Seed Funding Programme for Basic Research (Project number: 200911159065) of the University of Hong Kong.
This article has been published as part of BMC Proceedings Volume 5 Supplement 2, 2011: Proceedings of the 6th International Symposium on Bioinformatics Research and Applications (ISBRA'10). The full contents of the supplement are available online at http://www.biomedcentral.com/1753-6561/5?issue=S2.
Authors’ Affiliations
References
- Frank DN, Pace NR: Ribonuclease P: unity and diversity in a tRNA processing ribozyme. Annu. Rev. Biochem. 1998, 67: 153-180. 10.1146/annurev.biochem.67.1.153.View ArticlePubMedGoogle Scholar
- Nguyen VT, Kiss T, Michels AA, Bensaude O: 7SK small nuclear RNA blinds to and inhibits the activity of CDK9/cyclin T complexes. Nature. 2001, 414: 322-325. 10.1038/35104581.View ArticlePubMedGoogle Scholar
- Wadler CS, Vanderpool CK: A dual function for a bacterial small RNA: SgrS performs base pairing-dependent regulation and encodes a functional polypeptide. Proc Natl Acad Sci USA. 2007, 104 (51): 20454-9. 10.1073/pnas.0708102104.PubMed CentralView ArticlePubMedGoogle Scholar
- Yang Z, Zhu Q, Luo K, Zhou Q: The 7SK small nuclear RNA inhibits the CDK9/cyclin T1 kinase to control transcription. Nature. 2001, 414: 317-322. 10.1038/35104575.View ArticlePubMedGoogle Scholar
- Liu C, Bai B, Skogerbo G, Cai L, Deng W, Zhang Y, Bu D, Zhao Y, Chen R: NONCODE: an integrated knowledge database of non-coding RNAs. NAR. 2005, 33 (Database issue): D112-D115. 10.1093/nar/gki041.PubMed CentralView ArticlePubMedGoogle Scholar
- Griffiths-Jones S, Bateman A, Marshall M, Khann A, Eddy SR: Rfam: an RNA family database. NAR. 2003, 31: 439-441. 10.1093/nar/gkg006. [Http://www.sanger.ac.uk/Software/Rfam/]PubMed CentralView ArticlePubMedGoogle Scholar
- Eddy SR: Non-coding RNA genes and the modern RNA world. Nat Rev Genet. 2001, 2 (12): 919-929. 10.1038/35103511.View ArticlePubMedGoogle Scholar
- Hertel J, de Jong D, Marz M, Rose D, Tafer H, Tanzer A, Schierwater B, Stadler PF: Non-coding RNA annotation of the genome of Trichoplax adhaerens. Nucleic Acids Research. 2009, 37 (5): 1602-1615. 10.1093/nar/gkn1084.PubMed CentralView ArticlePubMedGoogle Scholar
- Klein RJ, Eddy SR: RSEARCH: Finding homologs of single structured RNA sequences. BMC Bioinformatics. 2003, 4: 44-10.1186/1471-2105-4-44.PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang S, Haas B, Eskin E, Bafna V: Searching genomes for noncoding RNA using FastR. IEEE/ACM TCBB. 2005, 2: 4-Google Scholar
- Tabei Y, Asai K: A local multiple alignment method for detection of non-coding RNA sequences. Bioinformatics. 2009, [Doi:10.1093/bioinformatics/btp261]Google Scholar
- Will S, Reiche K, Hofacker IL, Stadler PF, Backofen R: Inferring noncoding RNA families and classes by means of genome-scale structure-based clustering. PLOS Computational Biology. 2007, 3 (4): e65-10.1371/journal.pcbi.0030065.PubMed CentralView ArticlePubMedGoogle Scholar
- Jiang T, Lin G, Ma B, Zhang K: A general edit distance between RNA structures. Journal of Computational Biology. 2002, 9 (2): 371-388. 10.1089/10665270252935511.View ArticlePubMedGoogle Scholar
- Lin GH, Chen ZZ, Jiang T, Wen J: The longest common subsequence problem for sequences with nested arc annotations. Journal of Computer and System Sciences. 2002, 65 (3): 465-480. 10.1016/S0022-0000(02)00004-1.View ArticleGoogle Scholar
- Nawrocki EP, Eddy SR: Query-Dependent Banding (QDB) for faster RNA similarity searches. PLoS Comput Biol. 2007, 3 (3): e56-10.1371/journal.pcbi.0030056.PubMed CentralView ArticlePubMedGoogle Scholar
- Le S, Chen J, Maizel J: Efficient searches for unusual folding regions in RNA sequences. Structure and Methods: Human Genome Initiative and DNA Recombination. 1990, Adenine Pr, 1: 127-130.Google Scholar
- Rivas E, Eddy SR: Secondary structure alone is generally not statistically significant for the detection of noncoding RNAs. Bioinformatics. 2000, 16 (7): 583-605. 10.1093/bioinformatics/16.7.583.View ArticlePubMedGoogle Scholar
- Mosig A, Zhu L, Stadler PF: Customized strategies for discovering distant ncRNA homologs. Briefings in Functional Genomics and Proteomics. 2009, [To appear]Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.