|Home | About | Journals | Submit | Contact Us | Français|
Noncoding RNAs have emerged as important key players in the cell. Understanding their surprisingly diverse range of functions is challenging for experimental and computational biology. Here, we review computational methods to analyze noncoding RNAs. The topics covered include basic and advanced techniques to predict RNA structures, annotation of noncoding RNAs in genomic data, mining RNA-seq data for novel transcripts and prediction of transcript structures, computational aspects of microRNAs, and database resources.
Noncoding RNAs (ncRNAs) are transcripts that are not translated to proteins but act as functional RNAs. Several well-known ncRNAs such as transfer RNAs or ribosomal RNAs can be found throughout the tree of life. They fulfill central functions in the cell and thus have been studied for a long time.
However, over the past years a few key discoveries have shown that ncRNAs have a much richer functional spectrum than anticipated1. The discovery of microRNAs for example changed our view of how genes are regulated2,3. Another surprising observation revealed by high-throughput methods is that in human 90% of the genome is transcribed at some time in some tissue4. Although the full extent and functional consequences of this pervasive transcription remains highly controversial5,6, the vast amount of transcripts produced suggests that many important ncRNA functions are yet to be discovered.
In particular, long noncoding RNAs (lncRNAs) – transcripts that can be several kilobases in length, spliced and processed like mRNAs but lack obvious coding potential – seem to be a rich source of novel functions7. All these ncRNAs have been suggested to form a hidden layer of regulation that is necessary to establish the complexity of eukaryotic genomes8.
Prokaryotic genomes also contain many surprises. Riboswitches9, small regulatory RNAs10, or completely unknown structured RNAs11 suggest that ncRNAs also form an important functional layer in bacteria.
Understanding the function of ncRNAs – in particular in the age of high-throughput experiments – is clearly not possible without computational approaches. Algorithms to annotate, organize and functionally characterize ncRNAs are of increasing relevance. In this paper, we give a broad overview of programs and resources to analyze many different aspects of ncRNAs (Fig. 1).
For many RNAs there is a close connection between structure and function. Having a good model of the structure of an RNA is thus critical and often the first clue towards elucidating its function. Determining the complete three dimensional structure (“tertiary structure”) of an RNA is a tedious and time-consuming undertaking. Computational methods – either completely de novo or assisted by experimental data – are therefore routinely used to predict structure models. A strong focus lies on the prediction of the secondary structure, i.e. the pattern of intramolecular base-pairs (A·U, G·C and G·U) typically formed in RNAs.
In RNAs, the secondary structural elements are responsible for most of the overall folding energy and can be seen as a coarse-grained approximation of the tertiary structure. This important biophysical property in combination with the fact that secondary structure can easily be formalized as a simple graph (Fig. 2), led to secondary structure being widely studied early on. One of the first attempts to approach the RNA folding problem (i.e. predicting the secondary structure from the primary sequence) was by Nussinov and Jacobson12. They proposed an algorithm to find the secondary structure with the maximum number of base pairs. It is one of the classical examples of dynamic programming algorithms in computational biology and all modern variants of folding algorithm essentially use the same principle (Box 1).
The Dynamic Programming (DP) paradigm is used for many algorithms related to RNA folding. DP breaks down a problem in smaller sub-problems to find the overall solution efficiently. Nussinov’s algorithm is a classical example of a DP algorithm. Let’s assume we want to find the minimum free energy E i,j between the positions i and j of a sequence and already know the solution for a sequence from i + 1 to j, i.e. a sequence that is one base shorter. The new base i can either be unpaired or forms a base-pair with some position k:
The base-pair k divides the problem into two smaller sub-problems, namely finding the solution for Ei+1,k−1 and Ek+1,j. We thus can find the solution using a recursive algorithm:
βi,k is the energy contribution for the base-pair i, k in this simplified energy model.
In practice, however, finding the structure with the maximum number of pairs does not give accurate results. Ideally, we want to find the structure of minimum folding energy. Since most of the folding energy in RNAs is contributed by stacking interactions between neighbouring base pairs, counting single base pairs is not sufficient. Therefore, current folding algorithms use a “nearest neighbour” or also called “loop-based” energy model. A structure is uniquely decomposed into substructural elements (stacked bases, hairpin-loops, bulges, interior-loops, and multi-way-junctions, Fig. 2A). The structural elements are assigned energies which add up to the total folding energy of the structure (Fig. 2B). The energy values are established empirically and typically come from systematic melting experiments on small synthetic RNAs. An up-to-date set of energy parameters is maintained by Douglas Turner’s lab13,14.
A dynamic programming algorithm to find the minimum free energy (MFE) for this more complex energy model was proposed by Zuker and Stiegler15 and forms the basis for modern prediction programs. The most common implementations used are UNAFold16, RNAfold of the Vienna RNA package17 and RNAstructure18. The accuracy of MFE predictions depend on the type of RNA. Although some RNAs can be predicted with high accuracy, in general one has to expect that roughly a third of the predicted base pairs are wrong and one third of true base pairs are missed. It is thus important to keep in mind that even the best currently available prediction methods only give a rough model of the structure. Tertiary interactions, protein context and other inherent limitations of the energy model are all sources of potential errors.
At room temperature, RNAs usually exist in an ensemble of different structures and the MFE structure is not necessarily the biologically relevant structure. There are different algorithms to predict suboptimal structures close to the MFE structure19,20. McCaskill’s algorithm21 allows one to calculate the partition function over all possible structures and subsequently the probability of a particular base pair in the thermodynamic ensemble. Considering the pair-probability matrix of all possible base pairs gives a more comprehensive view of the structural properties of an RNA than just the MFE prediction. It is also possible to obtain individual structure predictions from the pair-probability matrix, either by sampling22 or by finding a structure that maximizes the expected accuracy considering a weighting factor between sensitivity and specificity23,24.
An alternative to the thermodynamic approach of RNA folding is a probabilistic approach based on machine-learning principles. Instead of using energy parameters, folding parameters can be estimated from a training set of known structures and are used to predict structures of unknown sequences. There are several probabilistic frameworks to accomplish parameter estimation and prediction. Stochastic context-free grammars (SCFGs) are a generalization of Hidden Markov models that are widely used in bioinformatics (Box 2). SCFGs allow one to consider nested dependencies, a prerequisite to model RNA structure. They have been used successfully for homology search problems and consensus structure prediction (see below). Although SCFGs can be used for single sequence structure prediction25 they are not widely used for this problem. CONTRAFOLD23, an alternative machine learning approach based on conditional random fields, however, could establish itself as a serious alternative to thermodynamic methods. There are also hybrid approaches that try to enhance the thermodynamic parameters by training on known structures26,27.
Context free grammars (CFG) are concepts from formal language theory. In most simple terms, a CFG is a set production rules V → w where V represent a so-called nonterminal symbol that produces a string of terminal or non-terminal symbols w. An example of a simple RNA grammar would be S → aSâaSSaSS. The grammar has one type of non-terminal symbol S and one type of terminal symbols a A, C, G, T representing the bases. The grammar consists of production rules for unpaired and paired bases (aâ represent two complementary bases). This simple rules allow to produce all possible RNA secondary structures. A stochastic context free grammar (SCFG) extends CFGs by assigning probabilities to all production rules. In the case of our RNA grammar, a full parametrized SCFG would thus describe the probability distribution over all structures and sequences.
Structure probing experiments typically use enzymatic or chemical agents that specifically target paired or unpaired regions28. Most implementations of thermodynamic folding like UNAfold or RNAfold allow for the incorporation of this type of information by restricting the folding to structures consistent with the experimental constraints. As an alternative, experimental information can also be incorporated as “pseudo-energies” into the folding algorithm enforcing regions to be preferentially paired or unpaired reflecting the experimental evidence29–31. Recently, high-throughput sequencing techniques were used to scale up structure probing experiments massively32–34. Dealing with inherently noisy data of this type turned out to be challenging and is still an active field of research.
Another way to improve secondary structure prediction is to consider homologous sequences from related species. If two or more sequences share a common structure but have diverged on the sequence level, typical base substitution patterns that maintain the common structure can be observed (Fig. 3A). A consistent mutation changes one base (e.g. A·U ↔ G·U) while compensatory mutations change both bases in the base pair (e.g. A·U ↔ G·C or C·G ↔ G·C). Clearly, these patterns provide useful information to infer a secondary structure.
The simplest way to exploit this signal is to calculate a mutual information score to find columns that show highly correlated mutation patterns. This method led to surprisingly accurate structures for rRNAs as early as 30 years ago35. Since we rarely have such large datasets of RNAs with extremely conserved structures it is natural to combine co-variance analysis with classical folding algorithms. RNAalifold36 extends Zuker’s folding algorithm to multiple sequence alignments by averaging the energy contribution over the sequences and adding co-variance information in the form of “pseudo-energies”. A probabilistic alternative is Pfold37, which uses a simple stochastic context free grammar combined with an evolutionary model of sequence evolution to infer a consensus structure. PetFold38 extends Pfold by incorporating pair probabilities from thermodynamic folding and thus unifies evolutionary and thermodynamic information. A more recent program, TurboFold39 also uses thermodynamic folding and iteratively refines the energy parameters by incorporating pair probabilities from homologous sequences.
Even unaligned RNAs can provide more information about their common structure than a single sequence. Low sequence homology below of 60% sequence identity40 prohibits the sequence alignment-based approach of the previous section (Fig. 3B), since correct alignment requires information about the structure. Since structure predictions for single sequences are unreliable, folding the sequences followed by structure-based alignment can also fail.
Therefore, the most successful strategies fold and align the RNAs simultaneously. The first such algorithm41, by Sankoff, simultaneously optimizes the alignment’s edit distance and the free energies of both RNA structures applying a loop-based energy model. However, only recent advances made this strategy applicable to practical RNA analysis.
The first complete pairwise Sankoff-implementation Dynalign42 implements a loop-based energy model, but employs a simple banding technique for increasing efficiency. A further pairwise Sankoff-like tool is Foldalign43. Computing a loop-based energy model during the alignment, these algorithms are accurate but computationally expensive; in practice, they compensate this by strong sequence-based heuristic restrictions.
Several less expensive Sankoff-like algorithms are based on simplifications introduced by PMcomp44. PMcomp replaces the loop-based energy model by assigning “pseudo-energies” to single base pairs. This reduces the computational overhead significantly. By computing the pseudo-energies of base pairs from their probabilities in the structure ensembles of the single RNAs, accurate information from the loop-based energy model is fed back into the light-weight algorithm.
Pseudoknotted structures follow the same rules as other secondary structures, but allow non-tree like configurations, e.g. due to an additional level of nested base pairs or pairing between different hairpin loops (kissing hairpin) (Fig. 4). Pseudoknots are prevalent in many ncRNAs; still, most algorithms ignore them for technical reasons: pseudoknot-folding is computationally expensive and accurate empirical energy models are missing.
The run-time of pseudoknot-free structure prediction grows only with the cube of the sequence length. Unfortunately, when general pseudoknots are allowed, the computation time grows much faster, namely exponentially with the sequence length50. Consequently, finding exact solutions is intractable for all but very short RNAs. Note that the “principle of optimality”, which allows dynamic programming (Box 1) in the pseudoknot-free case, is not applicable in the case of general pseudoknots. By this principle, every optimal structure can be composed from optimal structures of its subsequences.
In practice, often heuristic methods are applicable. Among the numerous approaches are ILM53, HotKnots54, KnotSeeker55, and IPknot56. ILM53 applies the classic principle of “hierarchic folding”; it constructs a pseudoknotted structure by iteratively predicting a most likely stem, using pseudoknot-free prediction, which is then added to the structure. HotKnots54 refines the construction of the pseudoknotted structure from likely more general secondary structure elements. TurboKnot57 even predicts conserved pseudoknots from a set of homologous RNAs. Applying a topological classification of RNA structures58, TT2NE59 guarantees to find the best RNA structure regardless of pseudoknot complexity; however, this limits the length of treatable sequences.
Other algorithms restrict the types of pseudoknots, such that dynamic programming can be applied50–52,60–63. These algorithms differ in their computational complexity and the complexity of considered structures. Fig. 4 shows pseudoknots of different complexity.
Rivas and Eddy51 proposed the most general such algorithm. It predicts three-knots (Fig. 4C), but cannot predict more complex pseudoknots such as the one shown in Fig. 4D. Its large time and space requirements prohibit its application to large scale data analysis. The most efficient such algorithm52 has only the same space requirements as pseudoknot-free prediction; its run-time grows with the fourth power of the sequence length, adding only a linear factor over pseudoknot-free folding. However, it predicts only canonical pseudoknots (Fig. 4E), which are formed by two canonical stems: such stems are (i) composed of canonical base pairs (ii) “perfect”, i.e. they do not contain interior loops or bulges, and (iii) maximally extended, i.e. they cannot be extended by canonical base pairs. Further such algorithms are tailored for specific interesting pseudoknot-classes (Fig. 4A and Fig. 4B63). Möhl et al.64 recently managed to speed up such algorithms non-heuristically.
A further challenge of pseudoknot prediction is to find an accurate energy model. The established loop-based energy models for RNA are tailored for pseudoknot-free structures; to date, there are no empirical energy parameters for pseudoknotted structure elements.
Consequently, some algorithms consider only the simplistic case of base-pair maximization65,53. Although some authors argue that important entropy contributions in pseudoknots cannot be covered by a loop-based energy model66, most approaches extend the loop-based energy model for pseudoknot-loops51,67.
While secondary structure is strongly stabilizing the three-dimensional structure, the tertiary structure depends on stabilizing non-canonical base pairs and van-der-Waals interactions. Furthermore, pseudoknots impact the tertiary structure. Therefore, deriving the tertiary structure in a hierarchic way from predicted secondary structure is not straightforward.
There are two main ways to model tertiary RNA structure. One, template-based modeling, employs homology to other RNAs with known structures. The other, de novo prediction, computes structures from physical and knowledge based rules. For example, the MC-fold/MC-sym pipeline68 (Fig. 5) assembles fragments of experimentally determined three-dimensional structures. Based on such structures, the approach builds a library of frequent small secondary structure loop-motifs, called Nucleic Cyclic Motifs (NCMs), together with their three-dimensional configurations. Given a RNA sequence, MC-Fold constructs probable secondary structures by merging NCMs; in this process it assigns likely NCMs to subsequences. MC-sym assigns concrete 3D-structures to the NCMs to generate a consistent three-dimensional structure. The approach has been tested on 13 RNAs of an average size of 30 nucleotides. Running several hours per prediction, the known 3D-structures have been reproduced within 2.3Å at average68. Similar to MC-Fold, the recent RNAwolf69 predicts extended secondary structures considering non-canonical base pairs; the same work reports a more efficient, dynamic programming-based reimplementation of MC-Fold, which improves parameter estimation. A more detailed review and a systematic performance comparison of RNA tertiary structure prediction programs is provided by Laing et al.70.
Many ncRNAs interact with other RNAs by specific base-pairing; most prominently, microRNAs target the untranslated regions of mRNAs. Predicting RNA/RNA interactions can thus elucidate RNA interaction partners and potential functions.
Most generally, one aims to predict the secondary structure of the interaction complex of two RNAs consisting of intramolecular and intermolecular base-pairs (Fig. 6). Alkan et al.71 formalized the problem and showed that – similar to pseudoknot prediction – it cannot be solved efficiently. Therefore, several simplifications and heuristics have been proposed.
Most approaches restrict the possible structures of the interaction complex to enable efficient algorithms using dynamic programming. Fig. 6 shows interaction complexes from several restriction classes. The simplest approaches ignore intramolecular base-pairs and predict only the best set of interacting base-pairs (Fig. 6A); examples are RNAhybrid72 and RNAduplex73.
A more general approach optimizes intra- and intermolecular base-pairs simultaneously in a restricted structure space. “Co-folding” of RNAs, for example implemented by RNAcofold73, concatenates the two RNA sequences and predicts a pseudoknot-free structure for the concatenation. Co-folding leads to a very efficient algorithm but strongly restricts the space of possible structures, such that only external bases can interact (Fig. 6B).
The dynamic-programming algorithms71,74 that predict more general structures (Fig. 6C), forbidding only pseudoknots, crossing interaction, and zig-zags (Fig. 6D), are computationally as expensive as the most complex efficient pseudoknot prediction algorithm51; they are therefore rarely used in practice, albeit their efficiency has been improved recently75.
Several fast methods71,76,77 assume that interactions form in two steps: First, the RNA unfolds partially, which requires certain energy to open the intramolecular base-pairs. Second, the unfolded, now accessible, RNA hybridizes with its partner forming energetically favorable intermolecular base-pairs. RNAup76 computes the energies to unfold each subsequence in the single RNAs and combines the unfolding energies with the hybridization energies to approximate the energy of the interaction complex. IntaRNA77 optimizes this approach and extends it to screen large data sets for potential interaction targets.
Common structure prediction methods assume that the functional RNA structure can be identified solely based on the thermodynamic equilibrium without considering the dynamics of the folding process. Although the true impact of kinetics on functional RNA structures is still unknown, for example RNA-switches80 show the importance of understanding the RNA folding process.
Several groups have studied the folding process of RNAs (reviewed in66,83). The RNA folding process is commonly modeled using energy landscapes84. Such landscapes assign energies to single structures, or states, and define neighborship between states. RNAlocopt85 enables studying the Boltzmann ensemble of local optima in an RNA energy landscape. The folding process iteratively moves from one state to a neighbor; the move probability depends on the energy difference. Studying folding by simulation83 is expensive since it requires averaging over many trajectories. Because the exact model of folding as a Markov process can be solved only for small systems, many methods coarse-grain the energy landscape to enable the analysis of the process. For example “barrier-trees” represent the energy landscape as a tree of local minima connected by their saddle points82. BarMap81 generalizes coarse-graining to non-stationary scenarios like temperature changes or co-transcriptional folding.86 predicts RNA folding pathways based on motion planning techniques from robotics. Kinefold87 simulates single folding paths over seconds to minutes for sequences up to 400 bases.
Another major challenge in understanding the function of ncRNAs is to find and annotate them in complete genomes. We distinguish homology search, i.e. trying to identify new members of already known classes of ncRNAs, and de novo prediction with the aim to discover novel ncRNAs.
Although a general de novo ncRNA finder remains elusive, some progress has been made in the identification of structural RNAs, i.e. ncRNAs that rely on a defined secondary structure for their function.
As a first attempt, one could use normal folding algorithms such RNAfold and hope to find structural RNAs to be thermodynamically more stable than the genomic background. However, although on average structural RNAs are more stable than expected this approach is generally not significant enough to reliably distinguish true structural RNAs from the rest of the genome88,89. Comparative approaches that make use of evolutionary signatures in alignments of related sequences can improve the signal considerably.
The first program that used pairwise alignments to find structured RNAs was QRNA90. Based on stochastic context free grammars it could successfully identify novel ncRNAs in bacteria91,92. MSARI was the first algorithm applying the idea of finding conserved RNA structures to multiple sequence alignments93.
To screen larger genomes higher accuracy was necessary. RNAz94 analyzes multiple sequence alignments and combines evidence from structural conservation and thermodynamic stability. EvoFold95 searches for conserved secondary structures in multiple alignments using a phylogenetic stochastic context-free grammar. Both programs were used to map potential RNA secondary structures in the human95,96 and many other genomes (e.g.97).
Another approach that was used to detect conserved RNA secondary structures in bacteria98 is implemented in CMFinder99. CMFinder builds a covariance model from a set of unaligned sequences by iterative optimization.
A solution is to include structure information in the search. Several motif description languages have been developed that allow one to manually specify sequence and structure properties and subsequently use these patterns to search databases or genomic data. Examples of such descriptor based search algorithm are RNAMOT102 and RNAmotif103.
Another class of programs automatically create a description of a structural RNA from a structure annotated alignment. The most commonly used program of this class is INFERNAL that uses covariance models, a full probabilistic description of an RNA family based on stochastic context-free grammars104. The Rfam database (see below) is based on INFERNAL and provides a curated collection of such covariance models.
In addition to these generic homology search tools, there are several specialized programs for finding ncRNAs of a particular family such as tRNAs105, rRNAs106, snoRNAs107–109, tmRNAs110, signal recognition particle RNAs111.
A complication during ncRNA annotation is the fact that many transcripts appear to be noncoding but in fact have the potential to code for a protein112. For example, short open reading frames can be easily missed and biological ambiguities of transcripts that act both on the level of the RNA and protein can make the annotation difficult113,114. A good overview of different methods to assess the coding potential of RNAs is given by Frith et al. The benchmark study115 found comparative analysis to be one of the most promising approaches. Purifying selection on the protein sequence turned out to be a reliable indicator of coding potential and several programs were developed to exploit this feature90,116,117.
The advent of high-throughput RNA sequencing has provided a robust platform for the development and expansion of several transcriptome-level analyses. RNA-seq is the highly parallelized process of sequencing individual cDNA fragments created from a population of RNA molecules. Here we discuss three challenges that must be addressed to mine RNA-seq data for noncoding transcripts: (i) read mapping to a reference genome (or transcriptome) (ii) transcriptome reconstruction from mapped reads, and (ii) quantification of transcript levels.
The first stage in short-read sequencing data analysis is the alignment of the sequenced reads to a reference genome. The algorithmic details of short read mapping is beyond the scope of this review. Here, it is important to note that transcript reconstruction, requires so-called “spliced aligners”. Spliced aligners such as TopHat118, GSNAP119, and SpliceMap120 identify and map short reads that span exon-exon junctions. From these spliced alignments novel splicing events and subsequently new transcript models can be identified.
A key advantage of RNA-seq over traditional forms of RNA expression analysis is the fact that little to no a priori information on the presence of an RNA sequence is required and, in principle, all required information can be learned directly from the data. This advantage, however, is dependent on the ability to re-construct a transcriptome fragmented into millions of short reads. Common approaches to solving this jigsaw puzzle-like problem focus on one of two different strategies: (i) reference-guided assembly, or (ii) de novo assembly (Fig. 8).
With a reference-guided assembly, reads are initially aligned using a spliced aligner to a reference genome sequence. The requirement for gapped alignments allows for discovery of putative splice junctions at the locations in which a read maps to the reference with a gap across an appropriately sized genomic interval. The two most popular reference-guided transcriptome assembly tools, Scripture121 and Cufflinks122, both treat these gaps as candidate splice junctions, and use this information to construct a graph representation of the transcriptome. In the case of Scripture, the graph represents the exonic connectivity potential of the reference genome. Cufflinks creates independent graph models for each independent genomic interval assumed to be a putative “gene”. In eithercase, the various paths through these connectivity graphs represent independent transcript isoforms. Scripture will attempt to identify all possible paths through the graph that can be explained by the mapped reads for a given gene and in this regard is useful for identifying lowly expressed isoforms, but tends to produce more noise in highly spliced structures. In contrast, Cufflinks produces a set of isoforms that represent the most parsimonious paths that can explain the given mapped reads, which may not report some redundant (but true) isoforms, but does not overburden the results with false positives. Additionally, Cufflinks estimates the read coverage across the paths to assist the selection of the most parsimonious isoforms. The result of either of these two approaches is a reconstructed transcriptome, the detail of which is supported by the read sequences, abundance, and mappability to a reliable reference.
In contrast to these reference-guided approaches, Velvet123 and transABySS124 use the short-read sequences directly and attempt to construct contig-like transcripts124. This approach tends to be significantly more computationally intensive, but is essential in species that do not have a reliable reference genome, or in the case when the expected transcriptome can deviate significantly from the reference genome due to rearrangements.
In RNA-Seq, the number of individual sequenced fragments from a given transcript is used as a proxy for its abundance. Determinations of expression level can be coarsely determined at the gene-level125,126 using a pseudo-model that consists of either the most-abundant isoform model, an intersection model quantifying only the regions present in all predicted isoforms, or a union model. The intersection model has been shown to to reduce the ability to accurately determine differential expression and the union model can under-estimate expression for those genes with alternative splicing127,122. More accurately, gene-level estimates can be determined as the sum of isoform-level abundance estimates122,128 involving a likelihood function to model the various effects encountered in the sequencing process129. The result of fitting these models to the data is a maximum likelihood estimate of the isoform-level abundances for each gene. Gene-level abundance estimates are easily determined by summing the expression levels of individual isoforms.
RNA-seq expression values must be normalized to correct for inherent biases in the data. The “Reads Per Kilobase of transcript per Million mapped” (RPKM) has emerged as a standard metric for reporting of estimated abundance levels. This metric has the advantage of correcting for the two main sources of variability in RNA-seq data: the length of the transcript, and the depth of the libraries.
The robust quantification of transcript levels also allows one to study differential expression. Since most gene-level projections of abundance estimation result in a single RPKM value for each gene, it would be reasonable to directly use most of the many differential expression tests that have been developed for microarray analysis over the past few years. There are, however, additional benefits that can be gained from using RNA-seq data such as the ability to derive a distribution of abundance estimates from a given sample or set of samples. Short read mapping to a given genomic interval can be considered a counting problem, and many differential expression analyses initially attempted to fit read counts to either Poisson or Binomial distributions to determine enriched transcripts. These methods, however, fail to incorporate any information about biological variability. Several more recent applications such as Cuffdiff122, DESeq130, and EdgeR131 incorporate variance information from biological replicates in their differential expression models leading to more rigorous statistics.
MicroRNAs (miRNAs/miRs) are short endogenous regulatory non-coding RNAs found in eukaryotic cells, whose primary function is to post-transcriptionally repress genes132. miRNAs inhibit translation and promote mRNA degradation via sequence-specific binding to the 3′ UTR regions and coding sequences133–135. They are produced from hairpin precursors (pri-miRNAs) that are processed by Drosha to form a pre-miR hairpin and then by Dicer to generate one or more 18- to 23-nt mature microRNAs136. Mature microRNAs are then incorporated into RISC where they hybridize with target sites of the mRNA, which are complementary to the microRNA seed (positions 2-8), leading to post-transcriptional repression. Since their discovery, there has been much interest in the computational identification of miRNAs at a genome wide level using some combination of evolutionary conservation, hairpin structure, thermodynamic stability, genomic context, and more recently the presence of the mature miRNA in sequencing data137.
Some of the first approaches such as miRScan and snarloop use conservation and similarity to known microRNAs for the prediction new examples138,139. Other approaches such as miRSeeker incorporate microRNA-specific patterns of conservation, such as stronger conservation in the hairpin stem compared to the loop140. Additionally, patterns of conservation of target sites have been used to identify novel microRNAs141 and to refine annotations of known microRNAs142. Other approaches combine both sequence and structural alignments to find microRNA homologs143,144.
The secondary structure and thermodynamic stability are important features for the prediction of miRNAs, especially when they are not conserved or orthologs do not exist in known species. Because miRNAs need to form stable hairpins for their processing, many studies have used structural features for their prediction. It has been demonstrated that miRNAs are significantly more stable than randomized sequences of the same nucleotide or dinucleotide composition145, and many studies have used programs like RNAz to predictnovel microRNAs based on this characteristic feature146,147. Other studies have developed machine learning approaches that train classifiers on known miRNAs and subsequently identify novel, and in many cases nonconserved, microRNAs148–150.
Other approaches have looked for features in the surrounding genomic context for the prediction of novel miRNAs and for refining other predictions. Some early work helped to filter predictions with a characteristic motif upstream of and patterns of conservation flanking the pre-miR151. Other studies have used the fact that microRNAs tend to reside within polycistronic clusters of more than one miR to identify novel miRNAs3,152. Some approaches also make use of the fact that regions proximal to microRNAs tend to be devoid of other non-miR small RNAs and when they are flanked by other small RNAs such as miRNA offset RNAs (moRs) the separation from mature microRNA sequences is minimal153.
The analysis of the sequencing of size-selected cDNA libraries has proved to be the most reliable method for the identification of novel microRNAs, in most instances coupled with other features such as structure and conservation to enhance predictions, because it provides validation that the mature sequence is expressed. In addition to novel miRNAs, the analysis of high-throughput sequencing data of small RNAs has led to the elucidation of many other classes of small RNAs including endogenous siRNAs154, piRNAs155,156, and moRs157 among others. There are now a few publicly available software tools for the prediction of microRNAs from high-throughput sequencing data such as miRDeep, MIReNA, and miRTRAP158,153,159.
MicroRNA target prediction is another lively area of computational analysis related to microRNAs. Early approaches such as targetScan identify evolutionary conserved seed matches and later approaches such as PicTar have incorporated target site stability160,155. The topic is related to the problem of predicting RNA/RNA interactions, discussed above. For a comprehensive review on target prediction, see Bartel161.
There are many databases related to ncRNAs and we cannot cover all of them here. A more specialized review162 and the yearly database issue of Nucleic Acids Research163 are good resources to get a more detailed overview.
There are many highly specialized databases that collect RNAs of a specific class. Basically for all well known “classical” RNAs like tRNAs, rRNAs, snoRNAs, SRP RNAs, tmRNAs, group I or II introns a database is available162.
All newly identified ncRNA sequences are usually deposited in general sequence database such as Genbank. However, typically they are not systematically annotated in these databases and consequently a few other databases have emerged that systematically collect ncRNAs (NONCODE164, RNAdb165 fRNAdb166 lncRNAdb167).
Rfam is an important resource for structured RNAs and also includes structured regulatory elements in mRNAs168. It collects hand curated covariance models (see above) that are used to systematically search sequence databases for new members. As of writing this review, Rfam contains 1973 families with a total of 2,756,313 members (Fig. 9). All RNA families in Rfam are manually annotated.
The most extensive database for microRNA sequences, hairpins, and target sites is miRBase (http://www.mirbase.org)169. miRBase has seen rapid growth over the past few years (Fig. 9) and is the official repository and naming authority for newly discovered miRNAs. Other related databases include Tarbase, which is a database of experimentally verified target sites170, and miR2Disease, which is a database that maintains a manually curated set of disease associated microRNA target sites171.
The wide variety of topics covered in this paper reflects the increasing complexity of the field. It also clearly demonstrates the interdisciplinary effort that is necessary to address these problems. It is safe to predict that computational problems related to ncRNAs will remain challenging for the coming years. In particular elucidating the functions of lincRNAs will require new approaches. Many methods for structural analysis, for example, were developed for rather short structured ncRNAs and cannot be directly applied to lncRNAs that can be kilobases in length. Prediction of long range intramolecular interactions within lncRNAs or prediction of intermolecular RNA/RNA or RNA/DNA interactions of lncRNAs will require extensions and improvements of established algorithms. Also the problem of predicting protein-RNA interactions will be of high relevance given the increasing number of examples of lincRNAs that act as scaffolds for protein complexes. Also more accurate and efficient analysis of high-throughput data will be a challenge for the field. We have mentioned analysis of RNA-seq data, but next generation sequencing can also used for a variety of other ncRNA related problems such as high throughput RNA secondary structure probing or mapping RNA/protein interactions. Also new approaches to organize ncRNA data will be important and there is need for new centralized databases and specialized resources172.
This work was supported by the Austrian Fonds zur Förderung der Wissenschaftlichen Forschung [Schrödinger Fellowship J2966-B12 to StW] and by the German Research Foundation (grant WI 3628/1-1 to SW) and NIH award 1RC1CA147187 to DAH.