Bioinformatics Advance Access originally published online on April 25, 2007
Bioinformatics 2007 23(13):1588-1598; doi:10.1093/bioinformatics/btm146
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Murlet: a practical multiple alignment tool for structural RNA sequences
1Computational Biology Research Center, National Institute of Advanced Industrial Science and Technology (AIST), 2-42 Aomi, Koto-ku, Tokyo 135-0064, 2Graduate School of Information Sciences, Nara Institute of Science and Technology, 8916-5 Takayama-cho, Ikoma, Nara 630-0192 and 3Department of Computational Biology, Faculty of Frontier Science, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8561, Japan
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT 1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSION5 CONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
Motivation: Structural RNA genes exhibit unique evolutionary patterns that are designed to conserve their secondary structures; these patterns should be taken into account while constructing accurate multiple alignments of RNA genes. The Sankoff algorithm is a natural alignment algorithm that includes the effect of base-pair covariation in the alignment model. However, the extremely high computational cost of the Sankoff algorithm precludes its application to most RNA sequences.
Results: We propose an efficient algorithm for the multiple alignment of structural RNA sequences. Our algorithm is a variant of the Sankoff algorithm, and it uses an efficient scoring system that reduces the time and space requirements considerably without compromising on the alignment quality. First, our algorithm computes the match probability matrix that measures the alignability of each position pair between sequences as well as the base pairing probability matrix for each sequence. These probabilities are then combined to score the alignment using the Sankoff algorithm. By itself, our algorithm does not predict the consensus secondary structure of the alignment but uses external programs for the prediction. We demonstrate that both the alignment quality and the accuracy of the consensus secondary structure prediction from our alignment are the highest among the other programs examined. We also demonstrate that our algorithm can align relatively long RNA sequences such as the eukaryotic-type signal recognition particle RNA that is 
300 nt in length; multiple alignment of such sequences has not been possible by using other Sankoff-based algorithms. The algorithm is implemented in the software named ‘Murlet’.
Availability: The C++ source code of the Murlet software and the test dataset used in this study are available at http://www.ncrna.org/papers/Murlet/
Contact: kiryu-h{at}aist.go.jp
Supplementary information: Supplementary data are available at Bioinformatics online.
|
1 INTRODUCTION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSION5 CONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
Recent studies have revealed that a substantial number of RNA transcripts do not code protein sequences in higher eukaryotic cells (Carninci et al., 2005; Dunham et al., 2004; Okazaki et al., 2002), and the question of whether such transcripts have any functional roles in cellular processes has attracted considerable interest. The existence of conserved secondary structures among phylogenetic relatives indicates the functional importance of such transcripts; therefore, it would be extremely interesting to detect conserved secondary structures from multiple alignments of genomic sequences. The evolutionary process of a structural RNA gene has a unique characteristic that the substitutions of distant bases are correlated in order to conserve their stem structures; hence, multiple alignment methods should account for such substitution patterns to enable accurate detection of the conserved structures. The Sankoff algorithm (Sankoff, 1985) is an alignment algorithm that naturally includes the effect of base-pair covariation in the alignment model. However, it is not practical to use the original version of the Sankoff algorithm due to its prohibitive computational cost. Hence, there have been intensive studies that have investigated practical variations of the Sankoff algorithm in recent years (Dowell and Eddy, 2006; Gorodkin et al., 1997; Havgaard et al., 2005; Hofacker et al., 2004; Holmes, 2005; Mathews and Turner, 2002; Uzilov et al., 2006). The algorithms proposed in these studies can be broadly categorized into two groups depending on how the secondary structures are scored in the algorithm.
In the first group, the algorithms score the structures using the free energy parameters collected by the Turner group (Mathews et al., 1999). These algorithms have the advantage of relatively accurate structure predictions. However, it is difficult for these algorithms to combine the structure energy with the homology information consistently. This group comprises the pairwise alignment programs Dynalign (Mathews and Turner, 2002; Uzilov et al., 2006) and Foldalign (Havgaard et al., 2005), and the multiple alignment program PMMulti (Hofacker et al., 2004).
In the second group, the algorithms score the structures as a part of the probabilistic model called the pair stochastic context-free grammar (PSCFG). These algorithms have the advantage that the parameters that score both the alignments and structures are determined in a unified manner. However, these algorithms have a potential disadvantage that the accuracies of the structure models might be only modest when compared with those in the first group; this is due to the limitations of PSCFG. The second group comprises the pairwise alignment program Consan (Dowell and Eddy, 2006) and the multiple alignment program Stemloc (Holmes, 2005).
These algorithms provide a variety of methods to reduce the enormous computation costs. Dynalign restricts the dynamic programming (DP) region to a narrow band such that only similar positions of sequences are compared. Foldalign and PMMulti limit the difference of subsequence lengths that are compared with each other. Stemloc implements a general method for combining the constraints in the structure space and those in the alignment space, which are computed using a Waterman–Eggert style suboptimal alignment algorithm (Waterman and Eggert, 1987). Consan constrains the DP region by anchoring the points in the DP matrix that have very high posterior probabilities of alignment according to the computations by the pair hidden Markov model (PHMM).
However, the computational costs remain relatively high despite these approximations, and it is impractical to use these programs for aligning sequences that are longer than 200 bases (as shown in Fig. 4). Therefore, several studies have sought for algorithms that can circumvent the Sankoff algorithm for fast computation of common secondary structures. For example, the SCARNA program (Tabei et al., 2006) aligns a pair of stem candidate sets that are extracted from the base pairing probability matrices of sequences. The RNAcast program predicts secondary structures for unaligned sequences that have a common topology or consensus shape (Reeder and Giegerich, 2005), and the RNAmine algorithm (Hamada et al., 2006) provides comprehensive list of the frequent stem motif patterns from unaligned sequences.
In this article, we propose a practical method based on the Sankoff algorithm for aligning multiple RNA sequences. We show that both the alignment quality and the accuracy of the consensus structure prediction from our alignment are the highest among the existing alignment softwares. Additionally, we show that our algorithm can align relatively long RNA sequences that have not been computable by other Sankoff-based multiple alignment algorithms. The algorithm is implemented in the software ‘Murlet’.
|
2 SYSTEMS AND METHODS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSION5 CONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
2.1 The model
First, we describe our algorithm for a pairwise sequence alignment. Our heuristic score system for the Sankoff algorithm is derived on the basis of two principles.
The first principle is the extensive preprocessing before applying the Sankoff algorithm. In general, the alignment of structural RNA sequences requires simultaneous consideration of complex information, such as base substitution score, gap insertion cost, stacking energy and various loop energies. If all these elements are included in the Sankoff model, the computation time would become unmanageably slow. Therefore, we used the match probability 
and the base-pairing probability 
to score the alignments and structures.
The match probability 
is the posterior probability that sequence positions i and j will be matched in an alignment. The match probability is calculated by using the standard PHMM (Durbin et al., 1998), as shown in Figure 1.
|
|
is the posterior probability of an alignment path 
given sequences x and y. 
is the joint probability of generating the alignment path , and it is estimated by the product of the transition and emission probabilities of the PHMM model. 
is the set of alignment paths that pass through the point 
in the DP matrix as the match state. The sum of the denominator in the second line is across all the possible alignment paths. 
is calculated using the forward and backward algorithms. The computation of 
requires 
time and 
memory.
|
. M indicates the match state, and I and D indicate the insertion and deletion states, respectively.The base-pairing probability
is the probability that the pair positions i and k in the sequence forms a base pair, and it is calculated by using the McCaskill algorithm (McCaskill, 1990). |
|
denotes a secondary structure candidate of sequence x; 
, the secondary structure free energy that is computed using the energy parameters collected by the Turner group (Mathews et al., 1999); R, the gas constant; T, the temperature; Z(x), the partition function and 
, the set of all the secondary structures that have a base pair between i and k. We let 
denote the loop probability at position i.
|
| (1) |
The computation of 
requires 
time and 
memory.
Both 
and 
can be computed by much faster algorithms than the Sankoff algorithm and compactly represent complex information such as sequence homology and structure contexts. This enables us to keep the Sankoff model very simple. Since the quantities 
and 
do not include the effects of base-pair substitution, we also apply the base-pair substitution matrix 
to score the events of the base-pair substitution.
The second principle is the maximal expected accuracy (MEA) principle. Recent studies have shown that the accuracy of the sequence alignment and the secondary structure predictions based on the principle of the maximization of expected accuracy (Holmes and Durbin, 1998; Miyazawa, 1995) perform better than those made by the conventional maximal likelihood algorithms (Do et al., 2006; Knudsen and Hein, 2003; Pedersen et al., 2006). A simple application of the MEA principle to the Sankoff algorithm that has a probabilistic scoring model such as PSCFG can be defined by the following expected accuracy z:
|
|
is the posterior probability that the bases xi and yj are aligned as unpaired bases, and 
is the posterior probability that the base pairs 
and 
aligned with each other as stem forming pairs. The sum of the first term is taken across the unpaired base matches in the candidate alignment and that of the second is taken across the base-pair matches. However, the computation of such function is demanding because the corresponding probabilistic model requires a large number of states to express the complex homology and structure information as described earlier. Therefore, we have adopted an alternative heuristic score function (Equations 2–4) that is formally similar to the formula given earlier but can be computed more easily.
To provide a mathematical definition of our algorithm, we consider a consensus structure annotation 
for each pairwise alignment 
of length L, which consists of sequences x and y of lengths Lx and Ly, respectively.
|
|
is the set of alignment columns without gap characters, and 
is the set of pairs of match columns. We consider only those cases where all the base pairs are formed between the match columns. We also ignore pseudo-knotted structures. We assign a score eL to each loop column 
and a score eS to each column pair 
.
|
| (2) |
|
| (3) |
denotes an element of the base pair substitution matrix. 
L and S are constant coefficients.
For each alignment 
and its consensus structure candidate 
, our heuristic alignment score 
is defined as the sum of the loop match scores eL and the base pair match scores eS.
|
| (4) |
is obtained by taking the maximum (
) of the score among all the alignments and structures.
To compute the maximum of 
, we have adopted the following variant of the Sankoff algorithm.
|
| (5) |
After the DP computation, the maximum of the score is obtained by 
. The computation of Equation (5) requires 
time and 
memory.
Note that the alignment result is defined in terms of the score function 
that depends only on the alignment 
and the structure 
, and is independent of the grammar, or transition rules, of the parsing algorithm (Equation 5). We may use an arbitrary grammar to compute the alignment result provided that the grammar can parse all the alignments and structures and that it does not modify the score system. The latter condition implies that the model cannot have any transition scores and that the left and right emission scores have to be identical. The independence from a particular grammar also indicates that there are no problems with regard to the ambiguity of the grammar. For an ambiguous grammar, two or more parse trees may correspond to the same alignment and structure. Since the score is solely dependent on the alignment and structure, the choice of the parse trees depends upon the detailed order of computations. This indicates that the obtained parse tree has little relevance. However, the alignment and its associated structure are unique, and they are sufficient for our purpose.
In contrast, the computations of the match and pair probabilities are affected by the redundant enumeration of the alignment and structure. However, both the forward–backward algorithm of the model of Figure 1 and the McCaskill algorithm enumerate all the alignments and structures without redundancies. Thus, the whole algorithm is devoid of any redundancy problems.
2.2 Reduction of the DP region
Since the loop match score eL and the base pair match score eS are proportional to the match probability 
, we consider restricting the Lx x Ly DP region to a smaller region that includes all the positions with 
, where 
is a prespecified threshold value.
For each alignment 
, let 
denote the set of match positions in the alignment 
that satisfy 
.
|
|
For a given initial alignment path and a threshold value 
, we then define the restricted DP region as the smallest region in the DP matrix that satisfies the following conditions:
- The region is simply connected, i.e. the region has no holes.
- The region includes the initial alignment path.
- For each alignment path
with 
in the full DP region, there exists an alignment 
' in the restricted DP region that satisfies 
.
We have described the algorithm for computing the restricted DP region in the Supplementary Material. The third condition implies that if all the match probabilities 
that are not greater than are set to zero, then there always exists an alignment in the restricted DP region that has the same score as the optimal score 
in the full DP region. It implies that for a sufficiently low threshold value (we use = 0.0001 throughout the article), restriction of the DP region rarely causes missing of the optimal alignment.
If two sequences are highly similar, the match probabilities concentrate along a specific diagonal in the DP matrix and the reduction of the DP region is quite significant. As shown in the later section, the elapsed time and memory are drastically reduced for similar sequences.
Previous studies (Dowell and Eddy, 2006; Holmes, 2005) have also proposed the algorithms that restrict the DP region using PHMM. In particular, our reduction method is a special case of a more general method proposed by Holmes (2005). However, our method is different from the above-mentioned algorithms in two aspects. First, our method removes only those positions with very small match probability from the DP region rather than selecting only the positions with large match probability as in their methods. This is because the positions with even the slightest match probability might have a large contribution to the DP score of the Sankoff algorithm as the match probability and the score function of the Sankoff algorithm are expected to be only roughly correlated. Second, in our algorithm, the score system of the Sankoff algorithm is more closely related to that of PHMM that is used to reduce the DP region. As defined in the Equations (2) and (3), both the loop match score eL and the pair match score eS are proportional to the match probability 
. This ensures that the total DP score has no contribution from the positions with zero match probability 
. On the other hand, the score systems of the Sankoff algorithm and the PHMM used to restrict the DP region are not directly related in the earlier algorithms. Hence, it is possible that the total alignment score has a large contribution from the positions with diminishing match probabilities. This implies that their restriction methods are at a higher risk of omitting the positions that significantly contribute to the total alignment score.
For these reasons, our restriction method is expected to have less possibility of missing the optimal alignment when compared with those of the previously mentioned algorithms.
2.3 Approximation methods
Most of the alignment softwares based on the Sankoff algorithm provide optional parameters to approximate the DP and to strike a balance between the computational cost and the alignment accuracy (Gorodkin et al., 1997; Havgaard et al., 2005; Hofacker et al., 2004; Holmes, 2005; Mathews and Turner, 2002). Murlet provides two original approximations that constrain the DP region: the strip and skip approximations.
For a given initial alignment path, the strip approximation constrains the DP region to a strip region of fixed width 
around the alignment path. If the strip width is equal to one, then the resulting alignment after the DP computation is the same as the initial alignment, as in the QRNA software (Rivas and Eddy, 2001). If a diagonal path is specified as the initial alignment path, then the strip approximation corresponds to the band alignment that calculates only the region 
for row i and column j in the DP matrix.
The band approximation has been used in the previous version of Dynalign (Mathews and Turner, 2002). The limitation of the band approximation is that the band width cannot be smaller than the difference 
of two sequences. The approximation methods that are adopted by Foldalign and PMMulti also have a similar limitation. The recent version of Dynalign (Uzilov et al., 2006) has adopted an alternative definition of the band region 
that removes this limitation. The strip approximation is more general as compared to these approximations because the initial path can be arbitrarily far from the main diagonal of the DP matrix, and the strip width can be set to one irrespective of the difference of the sequence lengths.
If the restriction of the DP region by match probabilities is not applied, the strip approximation decreases the computational costs by 
times with respect to time and 
times with respect to memory.
The skip approximation constrains the points that are computed during the bifurcation transitions (the last line of Equation 5) to a restricted set of positions in the DP region.
|
| (6) |
That is, the bifurcation calculation is performed only when the end positions 
and 
are in the skip set 
, and the only case considered is the one where the mid position 
is in the skip set 
. The skip set 
is the set of grid positions in the DP region 
that is defined as follows.
|
|
is the set of integers, (i) is a point on the initial alignment path at row i, and 
is a given parameter. 
corresponds to the full bifurcation calculation in the DP region, and in the limit 
, the algorithm can only parse non-bifurcating stem structures similar to the earlier version of Foldalign (Gorodkin et al., 1997). The bifurcation part of computation, which requires 
time and 
memory, decreases by 
times with respect to time and 
times with respect to memory with the skip approximation.
If the skip size 
is three or more, the bifurcation part is not a dominant factor of computation for aligning sequences shorter than 500 bases. In such cases, the total memory consumption is dominated by the 
memory that stores the traceback pointers, for which Murlet requires only one byte per DP recursion. The total time consumption is dominated by the 
calculations of the first seven lines of Equation (5). The order of memory consumption is only 
for these calculations.
The skip approximation is considered because the occurrence frequency of bifurcations in the parse tree is small as compared to the lengths of the RNA sequences despite the fact that the bifurcation calculation is the most compute-intensive part of the Sankoff algorithm. However, the skip approximation may miss a few base pairs if two neighboring stems are close to each other and no skip points are placed between them.
For a given strip width and skip size , the DP region of the Sankoff algorithm is determined as follows (see Fig. 2): First, the initial alignment path is determined (Fig. 2a) by the following DP algorithm, which is an application of the MEA principle to the PHMM.
|
|
|
We refer to the alignment obtained by this computation as the PHMM-MEA alignment. Next, the DP region is constrained to the strip region around the initial alignment path (Fig. 2b). The DP region is further constrained by removing the side regions with low match probabilities
(Fig. 2c). Finally, the skip set 
is determined within the DP region using the initial alignment path (Fig. 2d).
It is tedious to determine the appropriate strip width and skip size for each sequence pair being aligned. Murlet estimates the allocated memory and the computational time for each pairwise alignment and automatically determines the strip width and skip size so that the DP region is maximal under the given memory and time limits specified by the user.
The computation time t is estimated by the following formula.
|
| (7) |
is the size of the 
memory that is required to store traceback information of the Sankoff algorithm, 
is the 
memory that is required to store the scores of the child states of the bifurcation transitions, a and b are fitting parameters, and 
is the estimated number of bifurcation calculations (see Equation 6).
Figure 3 shows a scatter plot of the estimated time (x-axis) and the real time (y-axis). We used the pairwise alignments derived from the dataset of Table 2. We varied the strip width from 
to 
and skip size from 1 to 5 and measured the elapsed time for the computation of the pairwise alignments. As observed in the figure, the computation time can be estimated with a reasonable accuracy.
|
2.4 Probabilistic consistency transformations
For three or more sequences in the same sequence family, Do et al. introduced the probabilistic consistency transformation (PCT) of match probability matrices (Do et al., 2005), which is defined by the formula,
|
|
are the match probabilities after the transformation. This computation requires 
time for N sequences of length L. By this transformation, the match probability 
is increased if there are positions in other sequences that are likely to match with both i and j, and it is decreased if there are no such positions. Thus, the transformation introduces the family specific homology information into the match probabilities.
Here, we propose the PCT of the base pairing probability matrices defined by the formula,
|
|
The computation requires 
time. The corresponding loop probabilities 
are computed by applying Equation (1) to 
. Then, 
assumes a value between 0 and 1.
|
| (8) |
This justifies the consideration of the transformed matrices 
as the pair probability matrices. The proof of the formula 8 is presented in the Supplementary Material. As in the case of match probabilities, the transformation introduces the family specific structure information into the base-pairing probabilities. We show in the later section that the PCT of the match probabilities considerably improves the alignment accuracy.
The PCTs of 
and 
are performed for the sparse matrix representations of the probability matrices to reduce the computation time.
2.5 Multiple alignment procedure
We now describe the multiple alignment procedure. First, the base pairing probability matrices and the match probability matrices are computed for each sequence and each pair of sequences, respectively. Next, PCT is performed for the match probabilities; subsequently, PCT of the base-pairing probabilities using the transformed match probabilities. The similarity between a pair of sequences is defined by the score of the Sankoff algorithm along the PHMM-MEA alignment path. Using this similarity measure, a guide tree is constructed by using the unweighted pair group method (UPGMA) clustering algorithm. The progressive alignment is then performed using the guide tree. To align the two groups of aligned sequences, the base-pairing probabilities are averaged across all the sequences of each group. Further, the match probabilities are averaged across all the pairs of sequences between the two groups. The base-pair substitution score 
in Equation (3) is computed as the sum of the corresponding values for all the pairs of sequences between the groups. We set the proportionality constants L and S (Equations 2 and 3) as dependent on the number of sequences N1 and N2 in the two groups as follows:
|
|
As shown in the Supplementary Material, all the examined multiple alignment programs that make the structure prediction are inferior to Pfold with regard to the accuracy of the predicted structures. It suggests that, at present, it is practical to distinguish between the issue of multiple alignment and that of consensus structure prediction and to use the specialized programs to resolve the latter. Therefore, Murlet does not predict the consensus structure and returns only the aligned sequences.
2.6 The dataset
We collected the test dataset from the Rfam7.0 database (Griffiths-Jones et al., 2003). We used only the hand-curated seed alignments with the consensus structures published in literatures. For each sequence family, we generated up to 1000 random combinations of 10 sequences. We then removed the alignments with mean pairwise sequence identity higher than 95%. Because we are considering the global multiple alignment problem, we removed the alignments that contained more than 30% of the total alignment characters as gap characters. We also removed the alignments that contained < 5% of the total alignment characters as gap characters because the algorithms that merely penalize or forbid the gap insertions show high accuracies for such alignments. We found it difficult to collect completely exclusive alignment set for several sequence families. Therefore, we removed only those alignments sharing more than 30% of sequences with another alignment. Inspecting the number of families and the number of sub-alignments available for each family, we chose the dataset shown in Table 2.
The dataset consists of 85 multiple alignments of 10 sequences. There are 17 sequence families, and there are five alignments for each family. The dataset is reasonably diverse; its mean length varies from 54 bases to 291 bases, and the mean pairwise sequence identities varies from 40 to 94%.
We also used the multiple alignments of BRAlibaseII benchmark dataset for the evaluation (Gardner et al., 2005). The dataset consists of 481 multiple alignments of 5 sequences that are composed of tRNA, Intron_gpII, 5S_rRNA, U5 families in the Rfam5.0 database, and the signal recognition particle RNA family (SRP) in the SRPDB database (Larsen and Zwieb, 1993). As shown in Table 1, approximately half of the alignments have more than 70% sequence identities and few alignments have sequence identities < 50%. Since their dataset does not contain consensus structure annotations to the alignments, we have extracted the consensus structures from the original databases. Since the secondary structures are annotated to all the sequences in SRPDB, we have defined the base pairs that are supported by four or more sequences in the alignment as the consensus base pairs.
|
2.7 Accuracy measures
The accuracy of the alignments is measured by the standard sum-of-pairs score (SPS) (Carillo and Lipman, 1988). To measure the efficiency of the structural alignment, the consensus structures are predicted from the alignment results using the Pfold program (Knudsen and Hein, 2003). The Matthews correlation coefficients (MCC) are then calculated for the predictions (Matthews, 1975). MCC is defined by the formula
|
|
We did not use the consensus structures predicted by Stemloc, PMMulti and RNAforester since the accuracies of their predictions are lower than those of Pfold (see the Supplementary Material).
Since we used the external program Pfold for the computation of MCC, the upper limit of the MCC values is bound by the efficiency of the Pfold program. Furthermore, the results may be skewed by the compatibility of the programs with the Pfold software. To compensate for these inconveniences in our MCC measurement, we also measured the efficiency of structural alignment using the novel indicators sum-of-stem-pairs score (SSS), sum-of-quadruples score (SQS) and pair-column score (PCS) that quantify how well the true stems are aligned to each other. These indicators do not depend on the structure predictions to the alignment results and only use the reference alignments with annotated structure and the subject alignments. They are regarded as analogous to SPS and the column score (or TC score) (Carillo and Lipman, 1988; Thompson et al., 1994), which are frequently used for the evaluation of sequence alignments.
The SQS is defined as the fraction of the count of the pairs of base pairs that are correctly aligned as observed in the reference alignment. The counts are computed for all the pairs of sequences. The base-pairing positions of each sequence are derived from the annotated consensus structure in the obvious manner. The SSS is defined similarly; however, the criterion of a count is less stringent and allows the match of base pairs at different alignment columns in the reference alignment. In other words, it counts one if a base pair is aligned to another base pair irrespective of their alignment columns in the reference alignment. SSS measures how well each stem is aligned to another stem in the multiple alignment solely on the basis of the structural annotation and its values are of practical importance for the consensus structure prediction. The PCS is the fraction of base-pairing columns that are correctly reproduced in the subject alignment. PCS is more strongly dependent on the number of aligned sequences as compared to SQS and SSS and indicates the reliability of the alignment at the level of whole columns.
SQS and PCS take values between 0 and 1, and they are equal to 1 if the subject alignment is identical to the reference alignment. SSS is also a non-negative number, and it is equal to 1 if the alignment is identical to the reference alignment. Additionally, it is 
1 if all the stem regions in the reference alignment do not contain gap characters. However, it might be >1 when two or more sequences have gap characters in the stem regions of the reference alignment. The mathematical definitions and examples of computations for these measures are presented in the Supplementary Material.
|
3 IMPLEMENTATION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSION5 CONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
Murlet was implemented using the C++ language. For the computation of the match probabilities, we used the ProbCons software (version 1.10) (Do et al., 2005). For the computation of the base-pairing probabilities, we used the RNAAlifold program of the Vienna RNA package (version 1.5) (Hofacker et al., 2002; Hofacker, 2003). The base pair substitution matrix was extracted from the Stemloc software in the DART package (Holmes, 2005). Experiments were performed on a cluster of Linux machines equipped with dual AMD Opteron 850 2.4 GHz processors and 6 GB RAM. Due to the formidable time and memory consumption of Stemloc and PMMulti for longer sequence families, we limited the time and the maximal resident physical memory of the process to 500 min and 3.5 GB, respectively. We terminated the computation if the process exceeded the time or memory limit. A stand-alone program of Pfold was obtained for the consensus structure prediction to the alignment results (courtesy of Dr B. Knudsen).
|
4 DISCUSSION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSION5 CONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
4.1 Comparison of programs
Table 2 shows a comparison of the accuracy of the alignment for various alignment algorithms. The first three columns indicate the Rfam family name, mean sequence length and mean pairwise percent identity. The remaining columns show the SPS and MCC values for various algorithms: ClustalW (Thompson et al., 1994) is based on the ordinary DP algorithm of sequence alignment that does not account for the secondary structure. ProbCons (Do et al., 2005) is based on the PHMM-MEA algorithm. Murlet, Stemloc (Holmes, 2005) and PMMulti (Hofacker et al., 2004) are based on the Sankoff algorithm.
|
In Reference (Reeder and Giegerich, 2005), a multiple structural alignment method was proposed as an alternative to the Sankoff algorithm. First, this method predicts the secondary structures that have the same topology or consensus shape for all the unaligned sequences; subsequently, it performs the progressive alignment for the sequences with structure annotation. The secondary structures are predicted by the RNAcast program and the alignments are computed by the RNAforester program, (Hochsmann et al., 2004). For the sake of brevity, we have indicated this method as ‘RNAcast’ in the following tables, though the efficiency depends on both the RNAforester program as well as the RNAcast program.
For Murlet, we set the time and memory limits for each pairwise alignment to 10 min and 2 GB, respectively. The other softwares were used with the default option. If some of the five alignments in the family did not return within the limits of 3.5 GB and 500 min, the fraction of the alignments returned is indicated within parentheses in Table 2.
The last five rows indicate the average values of SPS and MCC for each program. ‘Average (all)’ indicates the average values taken over all the families. ‘Average (Stemloc)’, ‘Average (PMMulti)’, ‘Average (RNAcast)’ and ‘Average (common)’ represent the average values across the partial alignment set for which Stemloc, PMMulti, RNAcast and all the programs returned results, respectively. The ratios of the number of alignments to the whole dataset are indicated in parentheses.
Table 2 shows that among the softwares examined, the performance of Murlet is the best in terms of both the alignment accuracy SPS and the accuracy of the structure prediction MCC. Although the SPS values of ProbCons and the MCC values of Stemloc are relatively close to those of Murlet, the MCC values of ProbCons and the SPS values of Stemloc are much lower than the corresponding values of Murlet. The table also shows that the accuracies of ClustalW, PMMulti and RNAcast are lower than those of the other programs. Within the time and memory limit, Stemloc and PMMulti could not align most of the RNA sequences that were longer than 150 bases. In almost all the cases, the failures of Stemloc and PMMulti are caused by excessive memory requirements.
For the present dataset, RNAcast frequently failed to identify any consensus structures from the sequences. We changed the optional parameter ‘–c’ from 10 (default) to 50, which corresponds to the inclusion of the suboptimal structures that have free energy up to 50% higher than the minimal free energy but the number of correctly returned data remained unchanged.
Table 3 shows the SPS and MCC values for the BRAlibaseII multiple alignment dataset. Although the SPS and MCC values are relatively high and the differences of scores among the programs are smaller than the dataset of Table 2, Murlet still shows the highest accuracies with regard to both the SPS and MCC values.
|
Table 4 shows a comparison of the SSS, SQS and PCS for different softwares. The test sets are the same as those in the last five rows of Table 2. The superiority of Murlet when compared with the other programs is more obvious with respect to these measures. Moreover, Murlet is the only Sankoff-based program that performs better than the PHMM-based ProbCons software in all the accuracy measures. The table indicates that Murlet is the best among the examined programs for the structural alignment of RNA sequences.
|
4.2 Reduction of time and memory
Figure 4 shows the memory and time consumption of the programs. Each data point corresponds to a sequence family shown in Table 2. The x-axis represents the mean sequence length of the sequence family, and the y-axes represent the maximal resident physical memory in MB (left) and the elapsed time in minutes (right). The memory and time consumptions of ClustalW, ProbCons and RNAcast are very small when compared with those of the Sankoff-based programs, and several points for these programs coincide in the figure. The memory consumption of Stemloc and PMMulti drastically increases for sequences that are longer than 100 bases, and these programs cannot align sequences above 200 nts within the limits. In contrast, Murlet can align 10 sequences of the SRP_euk_arch family of mean length 291, within a realistic memory (570 MB) and time (32 min).
|
Figure 5 shows the dependence of the reduction of time and memory requirements on the sequence identities. We used 188 multiple alignments of four sequences collected from the Hammerhead_3 ribozyme family in the Rfam database. We compared the estimated time and the allocated memory between the full DP region and those of the region reduced by the match probabilities. For all 188 alignments, the two cases returned exactly the same alignment results. The mean SPS and MCC values were 0.87 and 0.85, respectively. The ratios of time and memory were binned for each 5% segment of the sequence identity, and the mean value for each bin was plotted. The figure shows the general trend that the time and memory usage decreases with the sequence identity. In particular, for sequence identities larger than 60%, the time and memory requirements are several hundred times smaller than those in the full DP case.
|
4.3 Effects of probabilistic consistency transformations
Figure 6 shows the density plots of the match probability distribution. The probabilities in the left figure are computed using the forward–backward algorithm of PHMM. The sequences are taken from the tRNA family shown in Table 2. The figure on the right represents the probabilities after PCT. Although the dense regions are broadened by the transformation, they are still concentrated around the main diagonal of the DP matrix.
|
Figure 7 shows an example of the true secondary structure of tRNA (left) and the corresponding base pairing probability matrices (right). The base pairing probability matrix as computed by the McCaskill algorithm is shown in the lower-left part of the figure on the right and that obtained after the transformation is shown in the upper-right part of the matrix. As indicated by the arrow in the figure, the McCaskill algorithm fails to identify one of the four stems of tRNA. PCT corrects this failure by adding small probabilities to this region.
|
|
Table 5 shows the effects of the PCTs on the alignment accuracies. For all the measures, the accuracies are the highest when the transformation is performed on both the match and pair probabilities. Further, the PCT of the pair probabilities are more significant than that of the match probabilities, and the latter is only effective when the former is also performed. This indicates that McCaskill algorithm often predicts incorrect base pairs and this results in considerable degradation of the alignment quality.
It is known that alignment errors that occur in the earlier pairwise alignments during the progressive alignment method have a considerable impact on the final alignment result. Therefore, it is important to investigate whether the PCTs improve the alignment quality at the level of pairwise alignment.
Table 6 shows the improvement of the pairwise alignment accuracy with the use of PCTs. We have used the test set that consists of 85 pairs of sequences that are randomly selected from each of the multiple alignments of Table 2. The PCTs have been applied by using the other eight sequences that belong to the same multiple alignment. In order to show the levels of accuracy by comparison, we measured the accuracies of pairwise alignments for several pairwise alignment programs as well as the multiple alignment programs.
|
Foldalign was used with the option that restricts the maximal difference of the segment lengths that are compared to each other to 50 bases. Dynalign was used with the band width of 20 and the gap penalty 0.4 kcal/mol. Murlet was used with the same option as indicated in the multiple case. The other programs were used with their default option. As in the case of multiple alignment, we terminated the program if the computation time or memory exceeded the limits of 500 min and 3.5 GB, respectively. Only Murlet, ProbCons, and ClustalW returned all the data within these limits. The MCC values were calculated for both the Pfold predictions and the original consensus structure predictions (if available) and the better score is listed in Table 6. Only Dynalign demonstrated the better structure predictions than Pfold (original: 0.61, Pfold: 0.60). The first column of Table 6 shows the programs that are compared with Murlet. The second column shows the fraction of the number of the data returned correctly. The third column shows the mean SPS and MCC values for the programs of the first column. The total computation time in minutes is also shown in parentheses. The fourth column shows the pairwise alignment accuracy and the total computation time of Murlet for the dataset that is used to compute the values of the third column. The fifth column is similar except that PCT is applied to the match and base-pair probabilities of each pairwise dataset by using the other eight sequences that belong to the same multiple alignment of Table 2. Since the datasets are different for each row, the comparisons are meaningful only within the row.
Table 6 indicates that Consan is currently the best pairwise alignment program because only Consan shows better scores with regard to SPS and MCC when compared with those of Murlet without PCT. Although the computation by Murlet is very fast if the PCTs are not applied, the alignment accuracies are only modest due to the inaccurate estimation of the match and base-pair probabilities. In the presence of multiple sequences, the inference of probability matrices by Murlet is greatly enhanced by the PCTs, which makes the accuracy of pairwise alignment of Murlet comparable to the best pairwise alignment programs while keeping the computation time 60 times smaller than that of Consan.
Thus, the PCTs efficiently improve the quality of pairwise alignment by using the information of the other sequences, which results in the enhancement of the final multiple alignment.
|
5 CONCLUSION |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSION5 CONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
We have developed an efficient method to align multiple sequences of structural RNAs. First, the method computes the base-pairing probabilities and match probabilities. A simple Sankoff algorithm is then applied to obtain the final alignment by using these probabilities.
We have developed several novel ideas and methods in this article:
- The heuristic score system (Equation 3) that is inspired by the MEA algorithm and is represented by the product of the match probabilities and base-pair probabilities.
- Efficient reduction of the DP region using the match probabilities.
- The strip approximation that restricts the DP region to a strip region around the initial alignment path.
- The skip approximation that limits the compute-intensive bifurcation calculations.
- Dynamical determination of the strip width and the skip size based on the estimation of the memory and time consumption. This method has greatly improved the utility of the program.
- The PCT for the base-pairing probabilities, which has been proved essential for ensuring high alignment quality.
- Novel accuracy measures SSS, SQS and PCS for the structural RNA alignment that measure the fraction of the base pairs aligned in the same manner as observed in the reference alignment.
We have shown that our method has the highest accuracy among the examined programs with regard to both the alignment quality and the structure prediction from our alignment. We have also shown that our program can align relatively long (300 bases) RNA sequences within a realistic memory and time.
We have only optimized the proportionality constants of the loop match score and the stem match score; however, we have not optimized the pair substitution matrix and the parameters of the models that are used to calculate the match and pair probabilities. It will be interesting to investigate the application of machine-learning methods in order to optimize these parameters in an integrated manner.
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACT1 INTRODUCTION2 SYSTEMS AND METHODS3 IMPLEMENTATION4 DISCUSSION5 CONCLUSIONACKNOWLEDGEMENTSREFERENCES |
|---|
This work was partially supported by the ‘Functional RNA Project’ funded by the New Energy and Industrial Technology Development Organization (NEDO) of Japan. The authors thank their colleagues who worked on the project for their discussions and comments. H.K. thanks Michiaki Hamada for inspiring discussions.
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
Associate Editor: Martin Bishop
Received on January 24, 2007; revised on March 19, 2007; accepted on April 10, 2007
|
REFERENCES |
|---|