De novo prediction of structured RNAs from genomic sequences

Introduction

Non-coding RNAs (ncRNAs) are functional transcripts that do not encode proteins. A handful of examples, such as transfer and ribosomal RNAs, have been well known since the dawn of molecular biology, and probably have existed since the dawn of life. These few examples have critical and deeply central roles, but, ironically, elucidating the full spectrum of ncRNA activity has received relatively little attention. Within the last 10-15 years, however, a number of striking discoveries including RNA interference, microRNAs and riboswitches have demonstrated that RNAs have unexpectedly diverse, sophisticated and important roles in all living organisms, sparking renewed interest in the “modern RNA world” [1]. To hint at the scope of the issue, only 1.2% of the human genome encodes protein [2], but recent data suggest that 90% of the genome is transcribed on one or both strands, to some extent, at some time, in some tissue [3]. The significance of this observation remains unclear and controversial, but a growing body of evidence points to the presence of many functionally important non-coding transcripts [4]. Even if the bulk of the non-coding transcription is “noise,” this still provides a vast substrate upon which natural selection may have acted to generate a multitude of biologically important non-coding RNAs. Thus, even striking ncRNA-related discoveries, such as microRNAs, may only be the “tip of the iceberg.”

Interest in conducting computational screens for functional RNAs has risen with their increased prominence [5]. However, in contrast to protein coding genes, whose regular codon structure provides strong signals aiding their computational recognition within nucleotide sequences, the signals for ncRNAs are subtler. For example, many ncRNAs are believed to be processed out of longer primary transcripts—including introns of protein-coding genes—hence lack features such as proximal promoters in addition to codon structure [1]. There is, however, one general characteristic shared by many (but not all) known ncRNAs: they fold into complex shapes that are critical to function and thus conserved. While prediction of RNA 3D structure is less well developed than protein structure prediction, the prediction of RNA secondary structure—the set of intramolecular, largely Watson-Crick, base pairs defining the fundamental units from which the tertiary structure is formed—is reasonably well understood and computationally tractable. Furthermore, secondary structure also tends to be conserved, while primary sequence often evolves more rapidly (see Figure 1) and is of limited use except when searching for close homologs. These facts make secondary structure the key signal to be exploited in ncRNA prediction. Although secondary structure (henceforth, simply structure) prediction is tractable, it is not trivial (e.g., involving interactions between nucleotides at variable and sometimes large distances in the primary sequence), which makes the problem challenging intellectually and computationally [5].

Most RNA will fold into some secondary structure, but given the importance of secondary structure in functional ncRNAs, it is reasonable to ask whether they have more stable structures than random genomic sequences. A negative answer to this question was provided about a decade ago. That is, the folding energy of RNA genes is not easily distinguished from that of appropriately chosen background sequences, such as randomly shuffled versions of the RNA gene itself [6, 7]. Therefore, a simple screening approach that seeks unusually stable genomic segments, will generally not work. MicroRNAs, whose precursors form unusually stable stem-loops, are a notable exception [8, 9]

Although searching for possible structured elements in one nucleotide sequence is generally ineffective, as discussed above, searching in multiple (orthologous) sequences can be highly effective, since evolutionary conservation highlights functionally important regions of all kinds. Importantly, such searches leverage the rapidly increasing body of comparative genomic sequence data. A key issue, however, is that the evolutionary signature of an RNA gene is quite different from that of a protein-coding gene. In particular, as noted earlier, the nucleotide sequence of an ncRNA might evolve relatively rapidly, making identification and alignment of orthologous sequences difficult or impossible, especially when using tools that focus only on sequence. However, patterns of compensating base changes, for example, an A–U base pair in a human RNA sequence corresponding to a C–G pair in mouse, can provide evidence supporting both the existence of a conserved RNA structure (without requiring conserved sequence; see Figure 1) and insight into the structure itself. Indeed, short of X-ray crystallography, this type of comparative analysis, done carefully by human experts [10], has been the gold standard for RNA secondary structure prediction for more than forty years. In a nutshell, this highlights the key challenges in computational prediction of ncRNAs: to find orthologous regions, expose the common structure therein, and do so both rapidly and accurately.

The comparative approach is important for an additional reason. Over the next few years, we expect that emerging technologies such as high throughput sequencing of RNA (RNAseq [11]) will reveal the transcriptomes of many organisms with unprecedented depth and precision. Yet, given the extensive breadth of genomic transcription now observed at least in mammals, evidence of transcription can no longer be taken as proof of functional importance. Any observed transcript might just be “noise” or incompletely degraded detritus arising from the expression of some nearby, functionally important RNA. Furthermore, experimental protocols will remain limited with regard to the diversity of species, cell types, states, growth and stress conditions that are probed. Consequently, lack of measured expression of a given genomic segment is not proof of lack of function. In contrast, evolutionary conservation strongly suggests functional importance, whether or not expression has already been experimentally verified. This does not deny the value of experimental evidence, of course, nor deny the existence of functionally important species-specific ncRNAs, but merely argues that detection of evolutionarily conserved ncRNAs by comparative genomics is a powerful tool, and belongs in any effort to understand living systems.

Computational search for conserved RNA structure does have certain important limitations. We highlight two of them here, since they color much of what follows. First, these searches are computationally expensive, principally due to the nature of the underlying RNA folding algorithms [12, 13] that need to be applied to these multiple sequences. Even single-sequence folding algorithms have run times that grow as the cube of the sequence length. Applied naively, however fast screening of a sequence of one kilobase is, it will be a thousand times slower for 10 kilobases and a million times slower for 100 kilobases. Hence, all successful programs in this arena are carefully engineered to control runtime, which entails some, hopefully modest, loss of accuracy on long genomic sequences. For example, one simple, widely used strategy is the “sliding window” approach, wherein the genomic sequence is cut into multiple, overlapping, fixed-length segments (“windows”) that are processed separately. This obviously limits the cubic runtime penalty to the length of the window, but, unfortunately, also limits the lengths of discoverable structures and risks arbitrarily truncating them. Even using substantially more sophisticated techniques, genome-scale ncRNA analyses often consume tens to hundreds of computer-years. These high computational costs are one reason why ncRNA gene finding is still in its infancy.

The second significant limitation of the general searches for conserved RNA secondary structures we will describe is more conceptual. It is natural to want to think of each element discovered as a “non-coding RNA gene”, but the truth is more complex since the approaches described here might generate only a partial picture for each ncRNA. For example, technical limitations related to “window” boundaries or splicing could result in partial or fragmentary predictions. More intrinsically, some ncRNAs lack conserved secondary structures, or may have only patches of conserved structure embedded in longer, largely unstructured transcripts. Additionally, conserved, functionally important RNA structures, such as the selenocysteine insertion sequence (SECIS element), are known to exist in messenger RNAs, usually in their untranslated regions. Therefore, identification of RNA structures in genomic data should trigger post-processing steps and follow up experiments to more precisely characterize transcript boundaries (and function!). For reasons of simplicity, however, in the following we will refer to individual conserved structures as “ncRNA genes.”

This review focuses explicitly on computational prediction of non-coding RNA elements by comparative genomics, i.e., the discovery of conserved structured elements in multiple genomic sequences. Other methods for de novo prediction have succeeded in some contexts (e.g., exploiting organism-specific differences in mono- or dinucleotide frequencies of ncRNAs versus background [14-16]), but the comparative approach appears to be the most broadly applicable. We will say little about the related ncRNA homology search problem—finding new instances of a particular RNA family given one or more examples—but this equally important task comes with its own set of issues [17], especially the difficulty of finding homologs outside the phylogenetic range of known examples.

From RNA folding to gene finding

Even though RNA structure cannot reliably be detected by merely folding single sequences, the principles obtained from folding single sequences are fundamental and often constitute an implicit part of more elaborate methods. For example, to date, no large-scale RNA structure screens have accounted for so-called pseudoknots, since the underlying RNA folding algorithms do not do so. Without pseudoknots, RNA secondary structure can be represented by nested parentheses. For example, the hairpin structure of the sequence AAAAUUCGGCAAUUUU with base pairs exactly between A's and U's can be written as ((((((….)))))). Although most large RNA structures do contain pseudoknots, they make up a relatively small portion of the overall structure. Algorithms that account for these more general cases are dramatically more complex (and slower) for seemingly modest gains in accuracy.

Folding of sequences is typically performed through energy minimization, whereby different structural elements are associated with experimentally measured energy parameters. Well-known programs for this task are mfold [18] and RNAfold [19]. Rather than minimizing the free energy, folding can also be carried out in a probabilistic framework, wherein the structure with the highest probability is sought. This is achieved through a so-called Stochastic Context-Free Grammar (SCFG), see e.g., [20, 21], typically using probabilities extracted from curated rRNA and tRNA alignments [22, 23]. The probabilistic and energetic frameworks are intimately connected, since one expects energetically destabilizing substitutions to be evolutionarily disfavored. Indeed, recent work shows that statistics derived from good structural alignments can be used to improve upon experimentally measured energy parameters [24]. Both energy-based and SCFG based methods are now part of the core inventory of ncRNA gene prediction methods. These principles for folding single sequences are variously employed on multiple sequences as well, often together with schemes that simultaneously search for compensating base changes.

In silico screening for RNA structures

Predicting RNA structure in genomic sequence is of course closely related to the existing methods for RNA structure prediction. As indicated above, multiple sequences are needed to reliably predict RNA structures. Several strategies exist for RNA structure prediction based on multiple sequences [25], which can be loosely categorized as “align-first,” “fold-first” and “joint.” As the name suggests, align-first strategies start by aligning all sequences using standard multiple sequence alignment tools, followed by inference of their presumed common structure. That inference typically proceeds by some combination of folding energy prediction and detection of compensatory base changes, e.g., as quantified by high mutual information [20, 22] between pairs of columns in the alignment. These methods work best when sequence conservation is sufficiently low that compensatory changes are not rare, but high enough that they are generally bracketed by well-conserved patches, which constrain the changes to be correctly aligned (and hence visible as a high mutual information score between paired columns). Dual to align-first strategies, fold-first strategies fold individual sequences separately, and then align the structures to each other. These strategies can be expected to excel where sequence conservation is low, but structure is clear-cut, stable and well-conserved. Joint strategies, such as approaches based on Sankoff's algorithm [26], simultaneously align and fold their input sequences. In a sense, the align-first and fold-first strategies are aimed at opposite evolutionary extremes. The joint approach, with an appropriate joint model for sequence and structure evolution, subsumes the other two, as shown in [27]. In practice, however, both high computational costs and difficulties in appropriately specifying such an evolutionary model have hindered the joint methods, leaving viable niches for all three approaches. Current approaches for de novo genomic screening mostly utilize methods based on the fold-first and joint strategies, as exemplified in Figure 2, in which prediction is based on finding structures within relatively short sequences locally, rather than globally. Meanwhile, fold-first methods have mainly been applied in homology searches [28, 29].

Besides ignoring pseudoknots, another crucial limitation of current genomic screening methods is their inability to cope with structural variation. That is, they assume that the individual RNA sequences share a conserved structure with little variation. In particular, structural inserts, resulting in large length differences, are poorly handled by current programs. Such structural inserts are, however, frequently observed in evolutionarily old ncRNAs, such as RNaseP [30]. Even the state-of-the-art SCFG-based homology search program Infernal deals only indirectly with structural inserts through its ability to match a local substructure [31]. Recent programs for comparison of RNA structure, such as LocARNA [32], point out some ways to address the problem. Additionally, the presence of two nearby structural alignments might suggest an intervening structural insertion/deletion in some of the sequences, a situation that potentially could be detected in a post-processing step. However, none of these approaches fully overcomes this limitation of the current de novo screening methods.

In the following we provide an overview of the main methodologies applied for de novo RNA structure discovery. These are also listed in Table 1.

Conclusions and perspectives

To date, annotation in the genome databases relates almost exclusively to protein coding genes. In contrast, the computational screens for non-coding RNAs described here suggest that a large number of functional ncRNAs still remain to be found. This is consistent with the observation that most of the non-coding mammalian genome is transcribed. Computational de novo discovery of non-coding RNAs within genomic sequences is still in its infancy. Nevertheless, these screens provide a valuable starting point for subsequent studies of specific genomic regions and their functional characterization.

In silico screens are useful as the information generated can be correlated to other studies, or to obtain a more general picture of the putative RNA structure content. Currently, screens essentially only address RNA secondary structure, with the numerous limitations discussed above, leaving much room for improvement. Another major issue is their high computational cost (i.e., their time and memory requirements). Their high false discovery rates also make them inappropriate for systematic annotation and necessitate experimental follow-up. There is clearly a demand for faster, better computational approaches. Nonetheless, present methods have already demonstrated their value and have opened an entirely new branch within comparative genomics. For instance, simple extrapolation from the ENCODE studies suggests prediction of ∼2 million RNA structures within the human genome, albeit potentially with many structures per transcript, high FDR and other caveats discussed above. Computational methods should, however, be developed beyond searching for local RNA structures, in part because some ncRNAs might either have no structure (e.g., piRNAs [67]), or only partial structures (e.g., Evf-1 [68, 69]; see also [70]). Additionally, ncRNAs are known to vary considerably in size, ranging from ∼20 nucleotides to ∼100,000 nucleotides, contributing significantly to the challenge at hand. Fusion transcripts also pose a serious computational challenge that the community has not yet addressed. An example of such a transcript (the cDNA KIAA0510) involving segments from different chromosomes was reported a few years ago [71], a story echoed by the ENCODE project [3]. Experimental approaches continue to lead to novel discoveries. The lincRNAs [72] are one recent example. Identified by a combination of ChIP-Seq and custom tiling array technologies, they include 1600 large, multiexon RNAs, evidently under purifying selection and expressed in multiple mouse cell types. Information obtained from such studies can in turn be incorporated into computational methods, which then might lead to further discoveries of novel ncRNAs. Methods to predict mRNA-like ncRNAs [73] promise to contribute to this. Another challenge is to unravel the significance of conserved structures in UTRs of protein coding genes, which constitute a sizeable fraction of predictions. For example, CMfinder candidate r-6354 ([64], see http://genome.ku.dk/resources/cmf_encode/pages/candidate.php?id=r-6354) found in a 5′ UTR in several mammals, might have a role in regulating its own mRNA, but is also a predicted microRNA precursor. In this situation, comparing to known elements in UTR databases and/or attempting de novo motif finding in UTRs are plausible approaches [74, 75].

Putative ncRNAs are available from various sources, including the UCSC browser [60] and (with some degree of overlap) RNAdb [76]. In our online supplement at http://genome.ku.dk/resources/rsgs, we provide links to files defining the genomic locations of many of the ncRNAs candidates described here. These files may be directly passed to the UCSC genome browser, so that the candidates may be viewed in context with other genomic features such as gene predictions and sequence-based multispecies conservation information.

Novel strategies are also emerging for post-clustering of predicted RNA structures. In particular, LocARNA [32], a Sankoff-based approach for structural alignment, grouped structurally related ncRNAs in an RNAz screen of Ciona intestinalis and C. savignyi, identifying known RNA families, such as tRNAs, microRNAs, and spliceosomal RNAs, and furthermore suggested additional groups consisting of novel, so far unannotated RNAs. Moreover, clustering based on FoldalignM detects novel relationships among microRNA family members [77]. Another clustering approach using RNAforester suggests some of the structural characteristics that distinguish microRNA precursors from the enormous number of other transcribed RNA stem-loops in mammalian genomes [78].

The increased availability of genome sequences has exposed sufficient nucleotide variation to allow the detection of patterns of RNA structure in sequence based-alignments, but has also revealed that sequence-based alignments alone might not be sufficient for detection of ncRNAs. The computational screens performed so far have pushed the limits of comparative genomics and pointed (perhaps unsurprisingly) towards structural alignments as an essential corequisite for detection of ncRNAs. As we have seen, this is not only the case for those orthologous genomic regions whose sequences are poorly conserved across related organisms, but might also be true for small portions of otherwise well conserved regions. This in turn poses new challenges in addressing comparative genomics from an RNA perspective. Furthermore, despite remarkable progress in both computer hardware and algorithms, even faster tools are highly desirable. With this in mind, tools for integrating new genomic sequence with existing de novo ncRNA predictions would be valuable, both for annotating the new genome and for strengthening or dismissing previously marginal predictions. As these predictions only need to be carried out in limited regions, they should be fast.

Even though further advances of methods for predicting RNA structures within genomic sequences are expected, the field is also developing tools to predict RNA interactions, with many utilizing the same principles as for RNA folding. Being able to predict RNA interactions will be an essential component in assigning a putative functional context for predicted RNA structures in the future.