ConStruct: Improved construction of RNA consensus structures

Prediction of RNA structure as well as searches for homologues in large genomic sequence databases play a prominent role in the era of non-coding RNAs (ncRNAs). Structure prediction may provide insight into RNA function, and pattern-based database searches [1, 2] may reveal new homologues, without the need for time-consuming experiments. Prerequisite for these predictions and searches as well as for inference of phylogeny [3–5] is the existence of an alignment of RNA homologues correct in terms of both sequence and structure. Sequence alignment tools like CLUSTALW [6] often fail to align ncRNA sequences correctly, especially when sequence homology drops below 60 % [7]. One reason is that ncRNA sequences evolve by compensatory base pair changes and ncRNA homologues are more conserved in structure than in sequence. For example, structural elements like thermodynamically extrastable tetraloops (UNCG, GNRA) share no sequence similarity and therefore cannot be correctly aligned by pure sequence alignment programs. Even structure alignment programs (e. g. DYNALIGN [8], FOLDALIGN [9], PMMULTI [10] or STEMLOC [11]) do not necessarily produce high-quality alignments under all conditions [7]. Moreover, these approaches are computationally extremely demanding, not only because they are based on simplified versions of the Sankoff algorithm [12]. Thus, automatically generated alignments often need to be corrected or refined by hand, which is a complex and tedious task. To ease this task a few sophisticated RNA alignment editors exist, e. g. 4SALE [13], SARSE [14] or S2S [15]. One of these tools is CONSTRUCT (construct ion of RNA con sensus struct ures; [16]), which is not only an RNA alignment editor but also allows for a variety of consensus structure predictions.

Here we present the completely revised and largely extended version of this tool and demonstrate some of its new features. CONSTRUCT allows for generation of RNA alignments, which are correct in terms of sequence and structure, by combining thermodynamic RNA structure prediction, several measures for covariation, and any alignment method. By applying this combination, typical shortcomings inherent to the single methods are eliminated; that is, the need of covariation for many, sufficiently divergent sequences is reduced, and the quality of thermodynamic predictions is enhanced. In contrast to tools, which predict a consensus structure automatically from a fixed alignment (e. g. RNAALIFOLD [17] or ILM [18]), CONSTRUCT allows for an interactive modification and optimization of the alignment. The user is able to modify the alignment similar to other RNA alignment editors [19, 13, 14]; the consequences of any alignment modifications are, however, immediately visible in a dotplot showing the probability of all base pairs of all RNA structures of the alignment; i. e., the user is guided during the alignment correction. In addition the user can account for sequence and structure constraints during the correction process. Afterwards optimal and suboptimal consensus structures and tertiary interaction can be predicted using a variety of built-in methods and displayed in several ways.

ConStruct's approach to consensus structure prediction

In the following we will describe the basic approach of CONSTRUCT (see Fig. 1).

1. 1.

   First, an initial sequence alignment needs to be created by means of an alignment program of the user's choice (e. g. by a pure sequence alignment program like MAFFT [20], by a sequence+structure alignment program like STRAL [21] or STEMLOC [11], or by a pure structure alignment program like PROFILE-DYNALIGN [22]).

    
2. 2.

   "Thermodynamic base pairing probability matrices" of all sequences in the alignment are automatically generated by means of a front-end program (CS_FOLD) to RNAfOLD [23]. An alternative to these thermodynamic approaches is creation of dotplots [24] with a minimum length of helices and thermodynamic weighting of their base-pair composition [25].

    
3. 3.

   Gaps from the initial alignment (step 1) are inserted into the matrices (step 2), resulting in identically sized matrices that are superimposed, thus building a consensus matrix. Ideally, homologous base pairs should now possess identical positions. The base pair matrices for each single sequence as well as the consensus matrix are displayed in a graphical user interface (GUI; see green and blue dots in upper triangle of the consensus dotplots shown in Fig. 1 and 2). The probability of a thermodynamic consensus base pair (see red dots in the consensus dotplots) is calculated such that noise from individual base pairs not representing the consensus is reduced and over-representation of sequence families is avoided (for details see [16]). For an overview of colors used in CONSTRUCT see Table S4 in Additional file 1.

    

1. 4.
   The new version of CONSTRUCT now allows to compute either the mutual information content (MI; [26]) or the RNAALIFOLD covariation score [17] to measure the amount of covarying positions (joined nucleotide substitutions, compensatory base pair changes). The results are displayed in the GUI (lower triangle of consensus dotplots in Fig. 1 and Fig. 2), can be filtered and normalized (see below), and are later on used in conjunction with the thermodynamic base pair matrices as the basis for consensus structure prediction. The MI at two aligned nucleotide positions i and j is defined as:

   ${\text{MI}}_{ij}={\displaystyle \sum _{\text{X},\text{Y}}{f}_{ij}(\text{XY})}{\log }_b\frac{f_{ij}(\text{XY})}{f_i(\text{X})\cdot f_j(\text{Y})}$
    

where f _i (X) and f _j (Y) are the frequencies of the nucleotide types X ∈ {A,U,G,C}and Y ∈ {A,U,G,C} at aligned positions i and j, and f _ij (XY) is the joint frequency of finding X at i and Y at j. In addition, he user may apply a normalization method [27], which enhances separation of truly correlated positions from background correlations. That is done by dividing the MI by the joint entropy

$h_{ij}={\displaystyle \sum _{\text{X},\text{Y}}{f}_{ij}(\text{XY})}{\log }_b{f}_{ij}(\text{XY}),$

the upper bound of the MI. For statistical analysis of the MI, maximum likelihood or unbiased probability estimation [28] in nits (b = e) [26] or bits (b = 2) [29] are available.

In comparison to the MI, the covariation score implemented in RNAALIFOLD measures compensations in Watson-Crick and wobble base-pairs [17] only, which is advantageous during search for helices. The meaningfulness of this score can be further improved by taking stacking into account (as shown in [30]), which is also a built-in option of CONSTRUCT. For a comprehensive description of the RNAALIFOLD covariation measure we refer to [17] and [30].

1. 5.

   The alignment of the sequences is displayed in a separate window (see alignment editor in Fig. 1). Position of base pairs from the dot plots is coupled with the position of the corresponding nucleotides in the alignment; i. e., pointing with the mouse to a consensus base pair highlights the corresponding base pairs in the alignment with a color from white to red according to the individual base pairing probabilities (see also Fig. 2); pointing to a base-paired nucleotide in the alignment changes the color of the corresponding base pair in the dot plot. A selected region of a single sequence or multiple sequences, with a gap at either side, may be moved with the mouse towards the gap, and the dot plots are updated correspondingly. Helices not superimposed are easily detectable. Thus, the user is guided during the alignment correction. These functions of the GUI are extremely helpful while correcting positions of structural elements, which were misaligned during the initial sequence alignment (step 1).

    
2. 6.

   Consensus structure prediction is now based upon the weighted and filtered summation of the thermodynamic consensus dotplot (step 2) and the covariation dotplot (step 4), whereas the previous version of CONSTRUCT used the unfiltered thermodynamic consensus dotplot alone.

    

The probability p _c of a thermodynamic consensus base pair at positions i and j is given by

$p_c(i,j)={\left(\frac{{\displaystyle {\sum }_{s=1}^N{w}_s\cdot p_s{(i,j)}^{1/3}}}{{\displaystyle {\sum }_{s=1}^N{w}_s}}\right)}^3$

where N is the total number of sequences and 0 ≤ w _s ≤ 1 is the user-defined weight of sequence s. This weighting can be used to avoid over-representation of a closely related sequence family in comparison to other sequences. The exponentiation helps to reduce low pairing probabilities from individual sequences.

The linear combination of the thermodynamic and the covariation pairing probabilities

$\begin{array}{ccc}{P}_c(i,j)& =& w_{\text{TD}}\cdot \{\begin{array}{ll}{p}_c(i,j)\hfill & \text{if }{p}_c(i,j)>t_{\text{TD}}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\\ +& w_{\text{CV}}\cdot \{\begin{array}{ll}{\text{CV}}_{i,j}\hfill & {\text{if CV}}_{i,j}>t_{\text{CV}}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\end{array}$

allows for thresholds t and a relative weighting (w_TD + w_CV = 1) of thermodynamics and covariation. The thresholds serve to further reduce the statistical noise and to suppress false positive base pairs and can be adjusted by the user. According to our experience, use of only thermodynamics with w_TD = 1.0 and t_TD = 0.03 already results in sensitive and specific predictions of secondary structures. The MI will only give additional information when the alignment contains many and divergent sequences. For tertiary structure predictions or detection of non-canonical base pairs [31], however, the MI must be reasonably high, since thermodynamic data alone are not sufficient.

Prediction of secondary structure is performed by dynamic programming [32] maximizing the weighted combination of the thermodynamic and covariation pairing probability. The new version of CONSTRUCT also allows to predict suboptimal structures [33, 34]. More importantly, routines for predicting tertiary interactions (pseudoknots, triple pairs) by maximum weighted matching (MWM) [35] are now built into CONSTRUCT. Examples for both prediction types are given in the Results section.

The predicted consensus structures can be viewed directly or be stored in several formats (RNAVIZ [36], CONNECT [37] or RNAML [38]). Another newly added feature is the support for structure logos [39], which can be directly requested from within CONSTRUCT. Three different graphical representations are supported:

• The first representation is basically an alignment of the sequences, where the background of the nucleotides is colored according to the nucleotide's structural features (for examples see bottom of "Structure Prediction" in Figs. 1 and Fig. 3C). An additional text output describes the structural alignment in numerical form (numbers of base pairs, consensus base pair changes, mismatches, consensus base pairing probability, MI per helical position, and statistical significance of MI values by χ² test).

• The consensus structure–annotated by an individual sequence or the consensus sequence with or without alignment gaps–can be displayed in different, SQUIGGLES-like ways (for examples see structures in Fig. 2 left). Overlapping of helical regions may be avoided by user interaction; each helix is selectable by the mouse and may be rotated around the upstream loop. Base pair connections can be color annotated according to their probability.

• The consensus structure–annotated by an individual sequence or the consensus sequence–can be viewed as a circular graph (circles plot) with nucleotides as edges and base pairs connected by probability-annotated arcs (for examples see "Structure Prediction" in Fig. 1 and Fig. 3B). Such a plot allows for representation of tertiary interactions; i. e., crossing arcs denote pseudoknots. If chemical or enzymatic mapping data are available, the accessibility of nucleotides can be marked with small triangles.

The two time consuming steps 1 and 2 are executed only once, whereas the remaining steps are computed on the fly. Handling of less than 100 sequences of length below 500 nucleotides is done fluently on a standard desktop PC.

The tool CONSTRUCT combines thermodynamic and statistical methods to predict the consensus structure of a set of homologous RNA molecules. In this respect it is similar to, for example, RNAALIFOLD [17] or ILM [18]. Yet, these programs use fixed alignments as input for structure prediction, whereas CONSTRUCT also allows to correct potentially incorrect alignments beforehand. CONSTRUCT's graphical user interface (GUI) guides the user in optimizing/correcting the alignment with respect to a consensus structure by displaying all base pair probabilities corresponding to the alignment in a consensus dotplot.

Most functions used in CONSTRUCT to extract a consensus structure from a given alignment are well known from the literature [32, 33, 26, 35, 17, 30]. Given a reasonably good structural alignment CONSTRUCT is able to extract the correct consensus structure even without user intervention. Usually almost optimal results can be obtained by using thermodynamic pairing probabilities alone. If an alignment contains 5–10 sequences with an average pairwise sequence identity below ~70 %, then sensitivity as well as specificity [40] for secondary structure predictions are above 80 % (see Fig. S5 and Table S3 in Additional file 1; see also [30]). If additional information either from MI, normalized MI or RNAALIFOLD covariation is used, the mean accuracy of secondary structure predictions is not increased (see Additional file 1). A bonus is, however, the validation of predicted pairings by the statistical information. Nonetheless, covariation scores increase prediction quality for alignments with many sequences, especially when predicting unusual base pairs, and the MI is essential when predicting tertiary interactions (see Example 2 below).

RNA alignment and consensus structure prediction is still a circular problem: A consensus structure needs to be known to create a high-quality alignment, but a high-quality alignment is prerequisite for consensus structure prediction. Thus, despite recent advances on this field [52–55] the need for RNA alignment editors which allow manual refinement based on structural properties is still there. These editors are still widely used and become increasingly sophisticated (see e. g. [13] and [14]). Automatically created RNA alignments and corresponding consensus structure prediction can be optimized in most cases (see this manuscript and [14]).

Our tool, CONSTRUCT, guides the user in correcting structurally misaligned regions. Once the initial alignment is refined, CONSTRUCT is able to predict secondary as well as tertiary consensus structures with high sensitivity and specificity. CONSTRUCT has already been described as an effective and "most elegant" [56] tool for structure alignment generation and RNA structure prediction. One of its strength is the "elaborate GUI" [56] that allows for easy identification and correction of structurally misaligned regions, guides the user in correcting an initial RNA sequence alignment, and allows for setting proper weight and threshold parameters for consensus structure prediction. Structurally misaligned regions are readily identifiable in the thermodynamic consensus dotplot and can be corrected by means of the built-in alignment editor. The example shown in Fig. 2 is typical in the sense that with a pure sequence alignment the "correct" consensus structure is already detectable in the CONSTRUCT dotplot and necessary corrections of the initial alignment are quite obvious.

The gold standard approach for RNA consensus-structure prediction–Comparative Sequence Analysis using covariation and MI [28, 57]–needs many sequences that have to be nearly perfectly aligned, which in turn is almost impossible for most sequence sets. Even given the perfect and large alignment, predictions only based on MI often suffer from non-informative columns (due to either too high sequence conservation or too many gaps) in the alignment. Purely thermodynamic based prediction methods are usually fairly reliable and allow for structure prediction from a few (in the extreme case one) sequences, and gain specificity and sensitivity when more sequences are added. Yet (standard) thermodynamic approaches alone cannot detect tertiary interaction or non-canonical base pairs. By using the (pair-entropy normalized) MI, which makes explicitly no use of base pairing rules, or the RNAALIFOLD covariation function (including stacking), which acknowledges consistent base pair mutations, CONSTRUCT is also able to predict non-canonical base pairs and tertiary interactions (for example see Fig. 3C). The prediction of tertiary interactions or at least certain types of pseudoknots could in principle be enhanced by including data into CONSTRUCT from structure-prediction programs other than RNAfold (for a review of alternatives see [58]).

In case of additional structural knowledge, for example from chemical or enzymatic mapping [63, 64], the initial structure prediction by RNAFOLD can accordingly be constrained (see RNAfold manual and [65]) and thus incorporated into CONSTRUCT. If even information on the three-dimensional structure of one of the sequences from the set is available from X-ray or NMR analysis, the use of S2S [15] in addition to CONSTRUCT is advantageous.

CONSTRUCT version 3 is based on a previously published version [16] and has been rewritten and largely extended. The underlying interpreter (Tcl/Tk 8.4) was extended to speed up the application; the installation process uses the GNU autotools. We successfully tested CONSTRUCT under several Linux distributions (Ubuntu, Debian, SuSE, Red Hat Fedora) as well as under Mac OS X (using Fink). Input of sequences is format independent due to use of the SEQIO package [66]. Graphics output is produced in PostScript. Various output formats for structures and further processing are supported (e. g. RNAML, connect etc.). The CONSTRUCT package (source and Debian package) including a manual can be downloaded at http://www.biophys.uni-duesseldorf.de/construct3/.