Efficient computation of RNA folding dynamics

Assessing the Quality of Cotranscriptional Folding Simulations

Structural changes in RNAs are an important contributor to controlling gene expression not only at the posttranscriptional stage but also during transcription. A subclass of riboswitches and RNA thermometers located in the 5' region of the primary transcript regulates the downstream functional unit - usually an ORF - through premature termination of transcription. Not only such elements occur naturally, but they are also attractive devices in synthetic biology. The possibility to design such riboswitches or RNA thermometers is thus of considerable practical interest. Since these functional RNA elements act already during transcription, it is important to model and understand the dynamics of folding and, in particular, the formation of intermediate structures concurrently with transcription. Cotranscriptional folding simulations are therefore an important step to verify the functionality of design constructs before conducting expensive and labor-intensive wet lab experiments. For RNAs, full-fledged molecular dynamics simulations are far beyond practical reach because of both the size of the molecules and the timescales of interest. Even at the simplified level of secondary structures, further approximations are necessary. The BarMap approach is based on representing the secondary structure landscape for each individual transcription step by a coarse-grained representation that only retains a small set of low-energy local minima and the energy barriers between them. The folding dynamics between two transcriptional elongation steps is modeled as a Markov process on this representation. Maps between pairs of consecutive coarse-grained landscapes make it possible to follow the folding process as it changes in response to transcription elongation. In its original implementation, the BarMap software provides a general framework to investigate RNA folding dynamics on temporally changing landscapes. It is, however, difficult to use in particular for specific scenarios such as cotranscriptional folding. To overcome this limitation, we developed the user-friendly BarMap-QA pipeline described in detail in this contribution. It is illustrated here by an elaborate example that emphasizes the careful monitoring of several quality measures. Using an iterative workflow, a reliable and complete kinetics simulation of a synthetic, transcription-regulating riboswitch is obtained using minimal computational resources. All programs and scripts used in this contribution are free software and available for download as a source distribution for Linux® or as a platform-independent Docker® image including support for Apple macOS® and Microsoft Windows®.

Classified Dynamic Programming in RNA Structure Analysis

Analysis of the folding space of RNA generally suffers from its exponential size. With classified Dynamic Programming algorithms, it is possible to alleviate this burden and to analyse the folding space of RNA in great depth. Key to classified DP is that the search space is partitioned into classes based on an on-the-fly computed feature. A class-wise evaluation is then used to compute class-wide properties, such as the lowest free energy structure for each class, or aggregate properties, such as the class' probability. In this paper we describe the well-known shape and hishape abstraction of RNA structures, their power to help better understand RNA function and related methods that are based on these abstractions.

DrTransformer: heuristic cotranscriptional RNA folding using the nearest neighbor energy model

MOTIVATION

Folding during transcription can have an important influence on the structure and function of RNA molecules, as regions closer to the 5' end can fold into metastable structures before potentially stronger interactions with the 3' end become available. Thermodynamic RNA folding models are not suitable to predict structures that result from cotranscriptional folding, as they can only calculate properties of the equilibrium distribution. Other software packages that simulate the kinetic process of RNA folding during transcription exist, but they are mostly applicable for short sequences.

RESULTS

We present a new algorithm that tracks changes to the RNA secondary structure ensemble during transcription. At every transcription step, new representative local minima are identified, a neighborhood relation is defined and transition rates are estimated for kinetic simulations. After every simulation, a part of the ensemble is removed and the remainder is used to search for new representative structures. The presented algorithm is deterministic (up to numeric instabilities of simulations), fast (in comparison with existing methods), and it is capable of folding RNAs much longer than 200 nucleotides.

AVAILABILITY AND IMPLEMENTATION

This software is open-source and available at https://github.com/ViennaRNA/drtransformer.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Predicting the Structure of a Viroid : Structure, Structure Distribution, Consensus Structure, and Structure Drawing

Viroids are small non-coding RNAs that require a special sequence and structure to be replicated and transported by the host machinery. Many of these features can be predicted and later experimentally verified. Here, we will present workflows to predict viroid structures and draw the predicted structures in a pleasing and descriptive way using recently developed software.

Regulation of the ThiM riboswitch is facilitated by the trapped structure formed during transcription of the wild-type sequence

The ThiM riboswitch from Escherichia coli is a typical mRNA device that modulates downstream gene expression by sensing TPP. The helix-based RNA folding theory is used to investigate its detailed regulatory behaviors in cells. This RNA molecule is transcriptionally trapped in a state with the unstructured SD sequence in the absence of TPP, which induces downstream gene expression. As a key step to turn on gene expression, formation of this trapped state (the genetic ON state) highly depends on the co-transcriptional folding of its wild-type sequence. Instead of stabilities of the genetic ON and OFF states, the transcription rate, pause, and ligand levels are combined to affect the ThiM riboswitch-mediated gene regulation, which is consistent with a kinetic control model.

Review of machine learning methods for RNA secondary structure prediction

Secondary structure plays an important role in determining the function of noncoding RNAs. Hence, identifying RNA secondary structures is of great value to research. Computational prediction is a mainstream approach for predicting RNA secondary structure. Unfortunately, even though new methods have been proposed over the past 40 years, the performance of computational prediction methods has stagnated in the last decade. Recently, with the increasing availability of RNA structure data, new methods based on machine learning (ML) technologies, especially deep learning, have alleviated the issue. In this review, we provide a comprehensive overview of RNA secondary structure prediction methods based on ML technologies and a tabularized summary of the most important methods in this field. The current pending challenges in the field of RNA secondary structure prediction and future trends are also discussed.

Evolutionary conservation of RNA sequence and structure

An RNA structure prediction from a single‐sequence RNA folding program is not evidence for an RNA whose structure is important for function. Random sequences have plausible and complex predicted structures not easily distinguishable from those of structural RNAs. How to tell when an RNA has a conserved structure is a question that requires looking at the evolutionary signature left by the conserved RNA. This question is important not just for long noncoding RNAs which usually lack an identified function, but also for RNA binding protein motifs which can be single stranded RNAs or structures. Here we review recent advances using sequence and structural analysis to determine when RNA structure is conserved or not. Although covariation measures assess structural RNA conservation, one must distinguish covariation due to RNA structure from covariation due to independent phylogenetic substitutions. We review a statistical test to measure false positives expected under the null hypothesis of phylogenetic covariation alone (specificity). We also review a complementary test that measures power, that is, expected covariation derived from sequence variation alone (sensitivity). Power in the absence of covariation signals the absence of a conserved RNA structure. We analyze artifacts that falsely identify conserved RNA structure such as the misuse of programs that do not assess significance, the use of inappropriate statistics confounded by signals other than covariation, or misalignments that induce spurious covariation. Among artifacts that obscure the signal of a conserved RNA structure, we discuss the inclusion of pseudogenes in alignments which increase power but destroy covariation.

This article is categorized under:

RNA Structure and Dynamics > RNA Structure, Dynamics and Chemistry

RNA Evolution and Genomics > Computational Analyses of RNA

RNA Evolution and Genomics > RNA and Ribonucleoprotein Evolution

Simulation of Folding Kinetics for Aligned RNAs

Studying the folding kinetics of an RNA can provide insight into its function and is thus a valuable method for RNA analyses. Computational approaches to the simulation of folding kinetics suffer from the exponentially large folding space that needs to be evaluated. Here, we present a new approach that combines structure abstraction with evolutionary conservation to restrict the analysis to common parts of folding spaces of related RNAs. The resulting algorithm can recapitulate the folding kinetics known for single RNAs and is able to analyse even long RNAs in reasonable time. Our program RNAliHiKinetics is the first algorithm for the simulation of consensus folding kinetics and addresses a long-standing problem in a new and unique way.

RNAxplorer: Harnessing the Power of Guiding Potentials to Sample RNA Landscapes

Motivation

Predicting the folding dynamics of RNAs is a computationally difficult problem, first and foremost due to the combinatorial explosion of alternative structures in the folding space. Abstractions are therefore needed to simplify downstream analyses, and thus make them computationally tractable. This can be achieved by various structure sampling algorithms. However, current sampling methods are still time consuming and frequently fail to represent key elements of the folding space.

Method

We introduce RNAxplorer, a novel adaptive sampling method to efficiently explore the structure space of RNAs. RNAxplorer uses dynamic programming to perform an efficient Boltzmann sampling in the presence of guiding potentials, which are accumulated into pseudo-energy terms and reflect similarity to already well-sampled structures. This way, we effectively steer sampling toward underrepresented or unexplored regions of the structure space.

Results

We developed and applied different measures to benchmark our sampling methods against its competitors. Most of the measures show that RNAxplorer produces more diverse structure samples, yields rare conformations that may be inaccessible to other sampling methods and is better at finding the most relevant kinetic traps in the landscape. Thus, it produces a more representative coarse graining of the landscape, which is well suited to subsequently compute better approximations of RNA folding kinetics.

Availabilityand implementation

https://github.com/ViennaRNA/RNAxplorer/.

Supplementary information

Supplementary data are available at Bioinformatics online.

pourRNA-a time- and memory-efficient approach for the guided exploration of RNA energy landscapes

MOTIVATION

The folding dynamics of ribonucleic acids (RNAs) are typically studied via coarse-grained models of the underlying energy landscape to face the exponential growths of the RNA secondary structure space. Still, studies of exact folding kinetics based on gradient basin abstractions are currently limited to short sequence lengths due to vast memory requirements. In order to compute exact transition rates between gradient basins, state-of-the-art approaches apply global flooding schemes that require to memorize the whole structure space at once. pourRNA tackles this problem via local flooding techniques where memorization is limited to the structure ensembles of individual gradient basins.

RESULTS

Compared to the only available tool for exact gradient basin-based macro-state transition rates (namely barriers), pourRNA computes the same exact transition rates up to 10 times faster and requires two orders of magnitude less memory for sequences that are still computationally accessible for exhaustive enumeration. Parallelized computation as well as additional heuristics further speed up computations while still producing high-quality transition model approximations. The introduced heuristics enable a guided trade-off between model quality and required computational resources. We introduce and evaluate a macroscopic direct path heuristics to efficiently compute refolding energy barrier estimations for the co-transcriptionally trapped RNA sv11 of length 115 nt. Finally, we also show how pourRNA can be used to identify folding funnels and their respective energetically lowest minima.

AVAILABILITY AND IMPLEMENTATION

pourRNA is freely available at https://github.com/ViennaRNA/pourRNA.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Finding the direct optimal RNA barrier energy and improving pathways with an arbitrary energy model

Motivation

RNA folding kinetics plays an important role in the biological functions of RNA molecules. An important goal in the investigation of the kinetic behavior of RNAs is to find the folding pathway with the lowest energy barrier. For this purpose, most of the existing methods use heuristics because the number of possible pathways is huge even if only the shortest (direct) folding pathways are considered.

Results

In this study, we propose a new method using a best-first search strategy to efficiently compute the exact solution of the minimum barrier energy of direct pathways. Using our method, we can find the exact direct pathways within a Hamming distance of 20, whereas the previous methods even miss the exact short pathways. Moreover, our method can be used to improve the pathways found by existing methods for exploring indirect pathways.

Availability and implementation

The source code and datasets created and used in this research are available at https://github.com/eukaryo/czno.

Supplementary information

Supplementary data are available at Bioinformatics online.

Mutations that improve efficiency of a weak-link enzyme are rare compared to adaptive mutations elsewhere in the genome

New enzymes often evolve by gene amplification and divergence. Previous experimental studies have followed the evolutionary trajectory of an amplified gene, but have not considered mutations elsewhere in the genome when fitness is limited by an evolving gene. We have evolved a strain of Escherichia coli in which a secondary promiscuous activity has been recruited to serve an essential function. The gene encoding the ‘weak-link’ enzyme amplified in all eight populations, but mutations improving the newly needed activity occurred in only one. Most adaptive mutations occurred elsewhere in the genome. Some mutations increase expression of the enzyme upstream of the weak-link enzyme, pushing material through the dysfunctional metabolic pathway. Others enhance production of a co-substrate for a downstream enzyme, thereby pulling material through the pathway. Most of these latter mutations are detrimental in wild-type E. coli, and thus would require reversion or compensation once a sufficient new activity has evolved.

Ant colony optimization for predicting RNA folding pathways

RNA folding dynamics plays important roles in various functions of RNAs. To date, coarse-grained modeling has been successfully employed to simulate RNA folding dynamics on the energy landscape composed of secondary structures. In such a modeling, the energy barrier height between metastable structures is a key parameter that crucially affects the simulation results. Although a number of approaches ranging from the exact method to heuristic ones are available to predict the barrier heights, developing an efficient heuristic for this purpose is still an algorithmic challenge. We developed a novel RNA folding pathway prediction method, ACOfoldpath, based on Ant Colony Optimization (ACO). ACO is a widely used powerful combinatorial optimization algorithm inspired from the food-seeking behavior of ants. In ACOfoldpath, to accelerate the folding pathway prediction, we reduce the search space by utilizing originally devised structure generation rules. To evaluate the performance of the proposed method, we benchmarked ACOfoldpath on the known nineteen conformational RNA switches. As a result, ACOfoldpath successfully predicted folding pathways better than or comparable to the previous heuristics. The results of RNA folding dynamics simulations and pseudoknotted pathway predictions are also presented.

Automated sequence-level analysis of kinetics and thermodynamics for domain-level DNA strand-displacement systems

As an engineering material, DNA is well suited for the construction of biochemical circuits and systems, because it is simple enough that its interactions can be rationally designed using Watson-Crick base pairing rules, yet the design space is remarkably rich. When designing DNA systems, this simplicity permits using functional sections of each strand, called domains, without considering particular nucleotide sequences. However, the actual sequences used may have interactions not predicted at the domain-level abstraction, and new rigorous analysis techniques are needed to determine the extent to which the chosen sequences conform to the system's domain-level description. We have developed a computational method for verifying sequence-level systems by identifying discrepancies between the domain-level and sequence-level behaviour. This method takes a DNA system, as specified using the domain-level tool Peppercorn, and analyses data from the stochastic sequence-level simulator Multistrand and sequence-level thermodynamic analysis tool NUPACK to estimate important aspects of the system, such as reaction rate constants and secondary structure formation. These techniques, implemented as the Python package KinDA, will allow researchers to predict the kinetic and thermodynamic behaviour of domain-level systems after sequence assignment, as well as to detect violations of the intended behaviour.

RNA modifications in structure prediction - Status quo and future challenges

Chemical modifications of RNA nucleotides change their identity and characteristics and thus alter genetic and structural information encoded in the genomic DNA. tRNA and rRNA are probably the most heavily modified genes, and often depend on derivatization or isomerization of their nucleobases in order to correctly fold into their functional structures. Recent RNomics studies, however, report transcriptome wide RNA modification and suggest a more general regulation of structuredness of RNAs by this so called epitranscriptome. Modification seems to require specific substrate structures, which in turn are stabilized or destabilized and thus promote or inhibit refolding events of regulatory RNA structures. In this review, we revisit RNA modifications and the related structures from a computational point of view. We discuss known substrate structures, their properties such as sub-motifs as well as consequences of modifications on base pairing patterns and possible refolding events. Given that efficient RNA structure prediction methods for canonical base pairs have been established several decades ago, we review to what extend these methods allow the inclusion of modified nucleotides to model and study epitranscriptomic effects on RNA structures.

Computational approaches for alternative and transient secondary structures of ribonucleic acids

Transient and alternative structures of ribonucleic acids (RNAs) play essential roles in various regulatory processes, such as translation regulation in living cells. Because experimental analyses for RNA structures are difficult and time-consuming, computational approaches based on RNA secondary structures are promising. In this article, we review computational methods for detecting and analyzing transient/alternative secondary structures of RNAs, including static approaches based on probabilistic distributions of RNA secondary structures and dynamic approaches such as kinetic folding and folding pathway predictions.

Efficient computation of co-transcriptional RNA-ligand interaction dynamics

Riboswitches form an abundant class of cis-regulatory RNA elements that mediate gene expression by binding a small metabolite. For synthetic biology applications, they are becoming cheap and accessible systems for selectively triggering transcription or translation of downstream genes. Many riboswitches are kinetically controlled, hence knowledge of their co-transcriptional mechanisms is essential. We present here an efficient implementation for analyzing co-transcriptional RNA-ligand interaction dynamics. This approach allows for the first time to model concentration-dependent metabolite binding/unbinding kinetics. We exemplify this novel approach by means of the recently studied I-A 2^'-deoxyguanosine (2^'dG)-sensing riboswitch from Mesoplasma florum.

In silico design of ligand triggered RNA switches

This contribution sketches a work flow to design an RNA switch that is able to adapt two structural conformations in a ligand-dependent way. A well characterized RNA aptamer, i.e., knowing its K_d and adaptive structural features, is an essential ingredient of the described design process. We exemplify the principles using the well-known theophylline aptamer throughout this work. The aptamer in its ligand-binding competent structure represents one structural conformation of the switch while an alternative fold that disrupts the binding-competent structure forms the other conformation. To keep it simple we do not incorporate any regulatory mechanism to control transcription or translation. We elucidate a commonly used design process by explicitly dissecting and explaining the necessary steps in detail. We developed a novel objective function which specifies the mechanistics of this simple, ligand-triggered riboswitch and describe an extensive in silico analysis pipeline to evaluate important kinetic properties of the designed sequences. This protocol and the developed software can be easily extended or adapted to fit novel design scenarios and thus can serve as a template for future needs.

PATTERNA: transcriptome-wide search for functional RNA elements via structural data signatures

Establishing a link between RNA structure and function remains a great challenge in RNA biology. The emergence of high-throughput structure profiling experiments is revolutionizing our ability to decipher structure, yet principled approaches for extracting information on structural elements directly from these data sets are lacking. We present PATTERNA, an unsupervised pattern recognition algorithm that rapidly mines RNA structure motifs from profiling data. We demonstrate that PATTERNA detects motifs with an accuracy comparable to commonly used thermodynamic models and highlight its utility in automating data-directed structure modeling from large data sets. PATTERNA is versatile and compatible with diverse profiling techniques and experimental conditions.

Computational Methods for Modeling Aptamers and Designing Riboswitches

Riboswitches, which are located within certain noncoding RNA region perform functions as genetic "switches", regulating when and where genes are expressed in response to certain ligands. Understanding the numerous functions of riboswitches requires computation models to predict structures and structural changes of the aptamer domains. Although aptamers often form a complex structure, computational approaches, such as RNAComposer and Rosetta, have already been applied to model the tertiary (three-dimensional (3D)) structure for several aptamers. As structural changes in aptamers must be achieved within the certain time window for effective regulation, kinetics is another key point for understanding aptamer function in riboswitch-mediated gene regulation. The coarse-grained self-organized polymer (SOP) model using Langevin dynamics simulation has been successfully developed to investigate folding kinetics of aptamers, while their co-transcriptional folding kinetics can be modeled by the helix-based computational method and BarMap approach. Based on the known aptamers, the web server Riboswitch Calculator and other theoretical methods provide a new tool to design synthetic riboswitches. This review will represent an overview of these computational methods for modeling structure and kinetics of riboswitch aptamers and for designing riboswitches.

Tractable RNA-ligand interaction kinetics

Background

The binding of small ligands to RNA elements can cause substantial changes in the RNA structure. This constitutes an important, fast-acting mechanism of ligand-controlled transcriptional and translational gene regulation implemented by a wide variety of riboswitches. The associated refolding processes often cannot be explained by thermodynamic effects alone. Instead, they are governed by the kinetics of RNA folding. While the computational analysis of RNA folding can make use of well-established models of the thermodynamics of RNA structures formation, RNA–RNA interaction, and RNA–ligand interaction, kinetic effects pose fundamentally more challenging problems due to the enormous size of the conformation space. The analysis of the combined process of ligand binding and structure formation even for small RNAs is plagued by intractably large state spaces. Moreover, the interaction is concentration-dependent and thus is intrinsically non-linear. This precludes the direct transfer of the strategies previously used for the analysis of RNA folding kinetics.

Results

In our novel, computationally tractable approach to RNA–ligand kinetics, we overcome the two main difficulties by applying a gradient-based coarse graining to RNA–ligand systems and solving the process in a pseudo-first order approximation. The latter is well-justified for the most common case of ligand excess in RNA–ligand systems. We present the approach rigorously and discuss the parametrization of the model based on empirical data. The method supports the kinetic study of RNA–ligand systems, in particular at different ligand concentrations. As an example, we apply our approach to analyze the concentration dependence of the ligand response of the rationally designed, artificial theophylline riboswitch RS3.

Conclusion

This work demonstrates the tractability of the computational analysis of RNA–ligand interaction. Naturally, the model will profit as more accurate measurements of folding and binding parameters become available. Due to this work, computational analysis is available to support tasks like the design of riboswitches; our analysis of RS3 suggests strong co-transcriptional effects for this riboswitch.

The method used in this study is available online, cf. Section “Availability of data and materials”.

In silico approaches to RNA aptamer design

RNA aptamers are ribonucleic acids that bind to specific target molecules. An RNA aptamer for a disease-related protein has great potential for development into a new drug. However, huge time and cost investments are required to develop an RNA aptamer into a pharmaceutical. Recently, SELEX combined with high-throughput sequencers (i.e., HT-SELEX) has been widely used to select candidate RNA aptamers that bind to a target protein with high affinity and specificity. After candidate selection, further optimizations such as shortening and modifying candidate sequences are performed. In these steps, in silico approaches are expected to reduce the time and cost associated with aptamer drug development. In this article, we review existing in silico approaches to RNA aptamer development, including a method for ranking the candidates of RNA aptamers from HT-SELEX data, clustering a huge number of aptamer sequences, and finding motifs amidst a set of significant RNA aptamers. It is expected that further studies in addition to these methods will be utilized for in silico RNA aptamer design, permitting a minimal number of experiments to be performed through the utilization of sophisticated computational methods.

Design of Artificial Riboswitches as Biosensors

RNA aptamers readily recognize small organic molecules, polypeptides, as well as other nucleic acids in a highly specific manner. Many such aptamers have evolved as parts of regulatory systems in nature. Experimental selection techniques such as SELEX have been very successful in finding artificial aptamers for a wide variety of natural and synthetic ligands. Changes in structure and/or stability of aptamers upon ligand binding can propagate through larger RNA constructs and cause specific structural changes at distal positions. In turn, these may affect transcription, translation, splicing, or binding events. The RNA secondary structure model realistically describes both thermodynamic and kinetic aspects of RNA structure formation and refolding at a single, consistent level of modelling. Thus, this framework allows studying the function of natural riboswitches in silico. Moreover, it enables rationally designing artificial switches, combining essentially arbitrary sensors with a broad choice of read-out systems. Eventually, this approach sets the stage for constructing versatile biosensors.

RNA folding kinetics using Monte Carlo and Gillespie algorithms

RNA secondary structure folding kinetics is known to be important for the biological function of certain processes, such as the hok/sok system in E. coli. Although linear algebra provides an exact computational solution of secondary structure folding kinetics with respect to the Turner energy model for tiny ([Formula: see text]20 nt) RNA sequences, the folding kinetics for larger sequences can only be approximated by binning structures into macrostates in a coarse-grained model, or by repeatedly simulating secondary structure folding with either the Monte Carlo algorithm or the Gillespie algorithm. Here we investigate the relation between the Monte Carlo algorithm and the Gillespie algorithm. We prove that asymptotically, the expected time for a K-step trajectory of the Monte Carlo algorithm is equal to [Formula: see text] times that of the Gillespie algorithm, where [Formula: see text] denotes the Boltzmann expected network degree. If the network is regular (i.e. every node has the same degree), then the mean first passage time (MFPT) computed by the Monte Carlo algorithm is equal to MFPT computed by the Gillespie algorithm multiplied by [Formula: see text]; however, this is not true for non-regular networks. In particular, RNA secondary structure folding kinetics, as computed by the Monte Carlo algorithm, is not equal to the folding kinetics, as computed by the Gillespie algorithm, although the mean first passage times are roughly correlated. Simulation software for RNA secondary structure folding according to the Monte Carlo and Gillespie algorithms is publicly available, as is our software to compute the expected degree of the network of secondary structures of a given RNA sequence-see http://bioinformatics.bc.edu/clote/RNAexpNumNbors .

Efficient approximations of RNA kinetics landscape using non-redundant sampling

MOTIVATION

Kinetics is key to understand many phenomena involving RNAs, such as co-transcriptional folding and riboswitches. Exact out-of-equilibrium studies induce extreme computational demands, leading state-of-the-art methods to rely on approximated kinetics landscapes, obtained using sampling strategies that strive to generate the key landmarks of the landscape topology. However, such methods are impeded by a large level of redundancy within sampled sets. Such a redundancy is uninformative, and obfuscates important intermediate states, leading to an incomplete vision of RNA dynamics.

RESULTS

We introduce RNANR, a new set of algorithms for the exploration of RNA kinetics landscapes at the secondary structure level. RNANR considers locally optimal structures, a reduced set of RNA conformations, in order to focus its sampling on basins in the kinetic landscape. Along with an exhaustive enumeration, RNANR implements a novel non-redundant stochastic sampling, and offers a rich array of structural parameters. Our tests on both real and random RNAs reveal that RNANR allows to generate more unique structures in a given time than its competitors, and allows a deeper exploration of kinetics landscapes.

AVAILABILITY AND IMPLEMENTATION

RNANR is freely available at https://project.inria.fr/rnalands/rnanr .

CONTACT

yann.ponty@lix.polytechnique.fr.

NMR Structural Profiling of Transcriptional Intermediates Reveals Riboswitch Regulation by Metastable RNA Conformations

Gene repression induced by the formation of transcriptional terminators represents a prime example for the coupling of RNA synthesis, folding, and regulation. In this context, mapping the changes in available conformational space of transcription intermediates during RNA synthesis is important to understand riboswitch function. A majority of riboswitches, an important class of small metabolite-sensing regulatory RNAs, act as transcriptional regulators, but the dependence of ligand binding and the subsequent allosteric conformational switch on mRNA transcript length has not yet been investigated. We show a strict fine-tuning of binding and sequence-dependent alterations of conformational space by structural analysis of all relevant transcription intermediates at single-nucleotide resolution for the I-A type 2'dG-sensing riboswitch from Mesoplasma florum by NMR spectroscopy. Our results provide a general framework to dissect the coupling of synthesis and folding essential for riboswitch function, revealing the importance of metastable states for RNA-based gene regulation.

The Conserved Dcw Gene Cluster of R. sphaeroides Is Preceded by an Uncommonly Extended 5' Leader Featuring the sRNA UpsM

Cell division and cell wall synthesis mechanisms are similarly conserved among bacteria. Consequently some bacterial species have comparable sets of genes organized in the dcw (division andcellwall) gene cluster. Dcw genes, their regulation and their relative order within the cluster are outstandingly conserved among rod shaped and gram negative bacteria to ensure an efficient coordination of growth and division. A well studied representative is the dcw gene cluster of E. coli. The first promoter of the gene cluster (mraZ1p) gives rise to polycistronic transcripts containing a 38 nt long 5’ UTR followed by the first gene mraZ. Despite reported conservation we present evidence for a much longer 5’ UTR in the gram negative and rod shaped bacterium Rhodobacter sphaeroides and in the family of Rhodobacteraceae. This extended 268 nt long 5’ UTR comprises a Rho independent terminator, which in case of termination gives rise to a non-coding RNA (UpsM). This sRNA is conditionally cleaved by RNase E under stress conditions in an Hfq- and very likely target mRNA-dependent manner, implying its function in trans. These results raise the question for the regulatory function of this extended 5’ UTR. It might represent the rarely described case of a trans acting sRNA derived from a riboswitch with exclusive presence in the family of Rhodobacteraceae.

Novel Approach to Analyzing MFE of Noncoding RNA Sequences

Genomic studies have become noncoding RNA (ncRNA) centric after the study of different genomes provided enormous information on ncRNA over the past decades. The function of ncRNA is decided by its secondary structure, and across organisms, the secondary structure is more conserved than the sequence itself. In this study, the optimal secondary structure or the minimum free energy (MFE) structure of ncRNA was found based on the thermodynamic nearest neighbor model. MFE of over 2600 ncRNA sequences was analyzed in view of its signal properties. Mathematical models linking MFE to the signal properties were found for each of the four classes of ncRNA analyzed. MFE values computed with the proposed models were in concordance with those obtained with the standard web servers. A total of 95% of the sequences analyzed had deviation of MFE values within ±15% relative to those obtained from standard web servers.

Predicting RNA secondary structures from sequence and probing data

RNA secondary structures have proven essential for understanding the regulatory functions performed by RNA such as microRNAs, bacterial small RNAs, or riboswitches. This success is in part due to the availability of efficient computational methods for predicting RNA secondary structures. Recent advances focus on dealing with the inherent uncertainty of prediction by considering the ensemble of possible structures rather than the single most stable one. Moreover, the advent of high-throughput structural probing has spurred the development of computational methods that incorporate such experimental data as auxiliary information.

Computational Design of a Circular RNA with Prionlike Behavior

RNA molecules engineered to fold into predefined conformations have enabled the design of a multitude of functional RNA devices in the field of synthetic biology and nanotechnology. More complex designs require efficient computational methods, which need to consider not only equilibrium thermodynamics but also the kinetics of structure formation. Here we present a novel type of RNA design that mimics the behavior of prions, that is, sequences capable of interaction-triggered autocatalytic replication of conformations. Our design was computed with the ViennaRNA package and is based on circular RNA that embeds domains amenable to intermolecular kissing interactions.

Pseudoknots in RNA folding landscapes

MOTIVATION

The function of an RNA molecule is not only linked to its native structure, which is usually taken to be the ground state of its folding landscape, but also in many cases crucially depends on the details of the folding pathways such as stable folding intermediates or the timing of the folding process itself. To model and understand these processes, it is necessary to go beyond ground state structures. The study of rugged RNA folding landscapes holds the key to answer these questions. Efficient coarse-graining methods are required to reduce the intractably vast energy landscapes into condensed representations such as barrier trees or basin hopping graphs : BHG) that convey an approximate but comprehensive picture of the folding kinetics. So far, exact and heuristic coarse-graining methods have been mostly restricted to the pseudoknot-free secondary structures. Pseudoknots, which are common motifs and have been repeatedly hypothesized to play an important role in guiding folding trajectories, were usually excluded.

RESULTS

We generalize the BHG framework to include pseudoknotted RNA structures and systematically study the differences in predicted folding behavior depending on whether pseudoknotted structures are allowed to occur as folding intermediates or not. We observe that RNAs with pseudoknotted ground state structures tend to have more pseudoknotted folding intermediates than RNAs with pseudoknot-free ground state structures. The occurrence and influence of pseudoknotted intermediates on the folding pathway, however, appear to depend very strongly on the individual RNAs so that no general rule can be inferred.

AVAILABILITY AND IMPLEMENTATION

The algorithms described here are implemented in C++ as standalone programs. Its source code and Supplemental material can be freely downloaded from http://www.tbi.univie.ac.at/bhg.html.

CONTACT

qin@bioinf.uni-leipzig.de

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Network Properties of the Ensemble of RNA Structures

We describe the first dynamic programming algorithm that computes the expected degree for the network, or graph G = (V, E) of all secondary structures of a given RNA sequence a = a₁, …, a_n. Here, the nodes V correspond to all secondary structures of a, while an edge exists between nodes s, t if the secondary structure t can be obtained from s by adding, removing or shifting a base pair. Since secondary structure kinetics programs implement the Gillespie algorithm, which simulates a random walk on the network of secondary structures, the expected network degree may provide a better understanding of kinetics of RNA folding when allowing defect diffusion, helix zippering, and related conformation transformations. We determine the correlation between expected network degree, contact order, conformational entropy, and expected number of native contacts for a benchmarking dataset of RNAs. Source code is available at http://bioinformatics.bc.edu/clotelab/RNAexpNumNbors.

An implementation of the Gillespie algorithm for RNA kinetics with logarithmic time update

In this paper I outline a fast method called KFOLD for implementing the Gillepie algorithm to stochastically sample the folding kinetics of an RNA molecule at single base-pair resolution. In the same fashion as the KINFOLD algorithm, which also uses the Gillespie algorithm to predict folding kinetics, KFOLD stochastically chooses a new RNA secondary structure state that is accessible from the current state by a single base-pair addition/deletion following the Gillespie procedure. However, unlike KINFOLD, the KFOLD algorithm utilizes the fact that many of the base-pair addition/deletion reactions and their corresponding rates do not change between each step in the algorithm. This allows KFOLD to achieve a substantial speed-up in the time required to compute a prediction of the folding pathway and, for a fixed number of base-pair moves, performs logarithmically with sequence size. This increase in speed opens up the possibility of studying the kinetics of much longer RNA sequences at single base-pair resolution while also allowing for the RNA folding statistics of smaller RNA sequences to be computed much more quickly.

Fast, approximate kinetics of RNA folding

In this article, we introduce the software suite Hermes, which provides fast, novel algorithms for RNA secondary structure kinetics. Using the fast Fourier transform to efficiently compute the Boltzmann probability that a secondary structure S of a given RNA sequence has base pair distance x (resp. y) from reference structure A (resp. B), Hermes computes the exact kinetics of folding from A to B in this coarse-grained model. In particular, Hermes computes the mean first passage time from the transition probability matrix by using matrix inversion, and also computes the equilibrium time from the rate matrix by using spectral decomposition. Due to the model granularity and the speed of Hermes, it is capable of determining secondary structure refolding kinetics for large RNA sequences, beyond the range of other methods. Comparative benchmarking of Hermes with other methods indicates that Hermes provides refolding kinetics of accuracy suitable for use in the computational design of RNA, an important area of synthetic biology. Source code and documentation for Hermes are available.

Thermodynamic and kinetic folding of riboswitches

Riboswitches are structured RNA regulatory elements located in the 5'-UTRs of mRNAs. Ligand-binding induces a structural rearrangement in these RNA elements, effecting events in downstream located coding sequences. Since they do not require proteins for their functions, they are ideally suited for computational analysis using the toolbox of RNA structure prediction methods. By their very definition riboswitch function depends on structural change. Methods that consider only the thermodynamic equilibrium of an RNA are therefore of limited use. Instead, one needs to employ computationally more expensive methods that consider the energy landscape and the folding dynamics on that landscape. Moreover, for the important class of kinetic riboswitches, the mechanism of riboswitch function can only be understood in the context of co-transcriptional folding. We present a computational approach to simulate the dynamic behavior of riboswitches during co-transcriptional folding in the presence and absence of a ligand. Our investigations show that the abstraction level of RNA secondary structure in combination with a dynamic folding landscape approach is expressive enough to understand how riboswitches perform their function. We apply our approach to a experimentally validated theophylline-binding riboswitch.

RNA folding: structure prediction, folding kinetics and ion electrostatics

Beyond the "traditional" functions such as gene storage, transport and protein synthesis, recent discoveries reveal that RNAs have important "new" biological functions including the RNA silence and gene regulation of riboswitch. Such functions of noncoding RNAs are strongly coupled to the RNA structures and proper structure change, which naturally leads to the RNA folding problem including structure prediction and folding kinetics. Due to the polyanionic nature of RNAs, RNA folding structure, stability and kinetics are strongly coupled to the ion condition of solution. The main focus of this chapter is to review the recent progress in the three major aspects in RNA folding problem: structure prediction, folding kinetics and ion electrostatics. This chapter will introduce both the recent experimental and theoretical progress, while emphasize the theoretical modelling on the three aspects in RNA folding.

The ViennaRNA web services

The ViennaRNA package is a widely used collection of programs for thermodynamic RNA secondary structure prediction. Over the years, many additional tools have been developed building on the core programs of the package to also address issues related to noncoding RNA detection, RNA folding kinetics, or efficient sequence design considering RNA-RNA hybridizations. The ViennaRNA web services provide easy and user-friendly web access to these tools. This chapter describes how to use this online platform to perform tasks such as prediction of minimum free energy structures, prediction of RNA-RNA hybrids, or noncoding RNA detection. The ViennaRNA web services can be used free of charge and can be accessed via http://rna.tbi.univie.ac.at.

Memory-efficient RNA energy landscape exploration

MOTIVATION

Energy landscapes provide a valuable means for studying the folding dynamics of short RNA molecules in detail by modeling all possible structures and their transitions. Higher abstraction levels based on a macro-state decomposition of the landscape enable the study of larger systems; however, they are still restricted by huge memory requirements of exact approaches.

RESULTS

We present a highly parallelizable local enumeration scheme that enables the computation of exact macro-state transition models with highly reduced memory requirements. The approach is evaluated on RNA secondary structure landscapes using a gradient basin definition for macro-states. Furthermore, we demonstrate the need for exact transition models by comparing two barrier-based approaches, and perform a detailed investigation of gradient basins in RNA energy landscapes.

AVAILABILITY AND IMPLEMENTATION

Source code is part of the C++ Energy Landscape Library available at http://www.bioinf.uni-freiburg.de/Software/.

Basin Hopping Graph: a computational framework to characterize RNA folding landscapes

MOTIVATION

RNA folding is a complicated kinetic process. The minimum free energy structure provides only a static view of the most stable conformational state of the system. It is insufficient to give detailed insights into the dynamic behavior of RNAs. A sufficiently sophisticated analysis of the folding free energy landscape, however, can provide the relevant information.

RESULTS

We introduce the Basin Hopping Graph (BHG) as a novel coarse-grained model of folding landscapes. Each vertex of the BHG is a local minimum, which represents the corresponding basin in the landscape. Its edges connect basins when the direct transitions between them are 'energetically favorable'. Edge weights endcode the corresponding saddle heights and thus measure the difficulties of these favorable transitions. BHGs can be approximated accurately and efficiently for RNA molecules well beyond the length range accessible to enumerative algorithms.

AVAILABILITY AND IMPLEMENTATION

The algorithms described here are implemented in C++ as standalone programs. Its source code and supplemental material can be freely downloaded from http://www.tbi.univie.ac.at/bhg.html.

Phylogeny and evolution of RNA structure

Darwin's conviction that all living beings on Earth are related and the graph of relatedness is tree-shaped has been essentially confirmed by phylogenetic reconstruction first from morphology and later from data obtained by molecular sequencing. Limitations of the phylogenetic tree concept were recognized as more and more sequence information became available. The other path-breaking idea of Darwin, natural selection of fitter variants in populations, is cast into simple mathematical form and extended to mutation-selection dynamics. In this form the theory is directly applicable to RNA evolution in vitro and to virus evolution. Phylogeny and population dynamics of RNA provide complementary insights into evolution and the interplay between the two concepts will be pursued throughout this chapter. The two strategies for understanding evolution are ultimately related through the central paradigm of structural biology: sequence ⇒ structure ⇒ function. We elaborate on the state of the art in modeling both phylogeny and evolution of RNA driven by reproduction and mutation. Thereby the focus will be laid on models for phylogenetic sequence evolution as well as evolution and design of RNA structures with selected examples and notes on simulation methods. In the perspectives an attempt is made to combine molecular structure, population dynamics, and phylogeny in modeling evolution.

Analysing RNA-kinetics based on folding space abstraction

BACKGROUND

RNA molecules, especially non-coding RNAs, play vital roles in the cell and their biological functions are mostly determined by structural properties. Often, these properties are related to dynamic changes in the structure, as in the case of riboswitches, and thus the analysis of RNA folding kinetics is crucial for their study. Exact approaches to kinetic folding are computationally expensive and, thus, limited to short sequences. In a previous study, we introduced a position-specific abstraction based on helices which we termed helix index shapes (hishapes) and a hishape-based algorithm for near-optimal folding pathway computation, called HiPath. The combination of these approaches provides an abstract view of the folding space that offers information about the global features.

RESULTS

In this paper we present HiKinetics, an algorithm that can predict RNA folding kinetics for sequences up to several hundred nucleotides long. This algorithm is based on RNAHeliCes, which decomposes the folding space into abstract classes, namely hishapes, and an improved version of HiPath, namely HiPath2, which estimates plausible folding pathways that connect these classes. Furthermore, we analyse the relationship of hishapes to locally optimal structures, the results of which strengthen the use of the hishape abstraction for studying folding kinetics. Finally, we show the application of HiKinetics to the folding kinetics of two well-studied RNAs.

CONCLUSIONS

HiKinetics can calculate kinetic folding based on a novel hishape decomposition. HiKinetics, together with HiPath2 and RNAHeliCes, is available for download at http://www.cyanolab.de/software/RNAHeliCes.htm.

RNA folding pathways and kinetics using 2D energy landscapes

RNA folding pathways play an important role in various biological processes, such as (i) the hok/sok (host-killing/suppression of killing) system in E. coli to check for sufficient plasmid copy number, (ii) the conformational switch in spliced leader (SL) RNA from Leptomonas collosoma, which controls trans splicing of a portion of the '5 exon, and (iii) riboswitches--portions of the 5' untranslated region of messenger RNA that regulate genes by allostery. Since RNA folding pathways are determined by the energy landscape, we describe a novel algorithm, FFTbor2D, which computes the 2D projection of the energy landscape for a given RNA sequence. Given two metastable secondary structures A, B for a given RNA sequence, FFTbor2D computes the Boltzmann probability p(x, y) = Z(x,y)/Z that a secondary structure has base pair distance x from A and distance y from B. Using polynomial interpolationwith the fast Fourier transform,we compute p(x, y) in O(n(5)) time and O(n(2)) space, which is an improvement over an earlier method, which runs in O(n(7)) time and O(n(4)) space. FFTbor2D has potential applications in synthetic biology, where one might wish to design bistable switches having target metastable structures A, B with favorable pathway kinetics. By inverting the transition probability matrix determined from FFTbor2D output, we show that L. collosoma spliced leader RNA has larger mean first passage time from A to B on the 2D energy landscape, than 97.145% of 20,000 sequences, each having metastable structures A, B. Source code and binaries are freely available for download at http://bioinformatics.bc.edu/clotelab/FFTbor2D. The program FFTbor2D is implemented in C++, with optional OpenMP parallelization primitives.

Predicting RNA structure: advances and limitations

RNA secondary structures can be predicted using efficient algorithms. A widely used software package implementing a large number of computational methods is the ViennaRNA Package. This chapter describes how to use programs from the ViennaRNA Package to perform common tasks such as prediction of minimum free-energy structures, suboptimal structures, or base pairing probabilities, and generating secondary structure plots with reliability annotation. Moreover, we present recent methods to assess the folding kinetics of an RNA via 2D projections of the energy landscape, identification of local minima and energy barriers, or simulation of RNA folding as a Markov process.

Energy-directed RNA structure prediction

In this chapter we present the classic dynamic programming algorithms for RNA structure prediction by energy minimization, as well as variations of this approach that allow to compute suboptimal foldings, or even the partition function over all possible secondary structures. The latter are essential in order to deal with the inaccuracy of minimum free energy (MFE) structure prediction, and can be used, for example, to derive reliability measures that assign a confidence value to all or part of a predicted structure. In addition, we discuss recently proposed alternatives to the MFE criterion such as the use of maximum expected accuracy (MEA) or centroid structures. The dynamic programming algorithms implicitly assume that the RNA molecule is in thermodynamic equilibrium. However, especially for long RNAs, this need not be the case. In the last section we therefore discuss approaches for predicting RNA folding kinetics and co-transcriptional folding.

Accessibility of microRNA binding sites in metastable RNA secondary structures in the presence of SNPs

MOTIVATION

We study microRNA (miRNA) bindings to metastable RNA secondary structures close to minimum free energy conformations in the context of single nucleotide polymorphisms (SNPs) and messenger RNA (mRNA) concentration levels, i.e. whether features of miRNA bindings to metastable conformations could provide additional information supporting the differences in expression levels of the two sequences defined by a SNP. In our study, the instances [mRNA/3'UTR; SNP; miRNA] were selected based on strong expression level analyses, SNP locations within binding regions and the computationally feasible identification of metastable conformations.

RESULTS

We identified 14 basic cases [mRNA; SNP; miRNA] of 3' UTR-lengths ranging from 124 up to 1078 nt reported in recent literature, and we analyzed the number, structure and miRNA binding to metastable conformations within an energy offset above mfe conformations. For each of the 14 instances, the miRNA binding characteristics are determined by the corresponding STarMir output. Among the different parameters we introduced and analyzed, we found that three of them, related to the average depth and average opening energy of metastable conformations, may provide supporting information for a stronger separation between miRNA bindings to the two alleles defined by a given SNP.

AVAILABILITY AND IMPLEMENTATION

At http://kks.inf.kcl.ac.uk/MSbind.html the MSbind tool is available for calculating features of metastable conformations determined by putative miRNA binding sites.

Computational analysis of noncoding RNAs

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.

Efficient procedures for the numerical simulation of mid-size RNA kinetics

Motivation

Methods for simulating the kinetic folding of RNAs by numerically solving the chemical master equation have been developed since the late 90's, notably the programs Kinfold and Treekin with Barriers that are available in the Vienna RNA package. Our goal is to formulate extensions to the algorithms used, starting from the Gillespie algorithm, that will allow numerical simulations of mid-size (~ 60–150 nt) RNA kinetics in some practical cases where numerous distributions of folding times are desired. These extensions can contribute to analyses and predictions of RNA folding in biologically significant problems.

Results

By describing in a particular way the reduction of numerical simulations of RNA folding kinetics into the Gillespie stochastic simulation algorithm for chemical reactions, it is possible to formulate extensions to the basic algorithm that will exploit memoization and parallelism for efficient computations. These can be used to advance forward from the small examples demonstrated to larger examples of biological interest.

Software

The implementation that is described and used for the Gillespie algorithm is freely available by contacting the authors, noting that the efficient procedures suggested may also be applicable along with Vienna's Kinfold.

Approximating the set of local minima in partial RNA folding landscapes

MOTIVATION

We study a stochastic method for approximating the set of local minima in partial RNA folding landscapes associated with a bounded-distance neighbourhood of folding conformations. The conformations are limited to RNA secondary structures without pseudoknots. The method aims at exploring partial energy landscapes pL induced by folding simulations and their underlying neighbourhood relations. It combines an approximation of the number of local optima devised by Garnier and Kallel (2002) with a run-time estimation for identifying sets of local optima established by Reeves and Eremeev (2004).

RESULTS

The method is tested on nine sequences of length between 50 nt and 400 nt, which allows us to compare the results with data generated by RNAsubopt and subsequent barrier tree calculations. On the nine sequences, the method captures on average 92% of local minima with settings designed for a target of 95%. The run-time of the heuristic can be estimated by O(n(2)Dνlnν), where n is the sequence length, ν is the number of local minima in the partial landscape pL under consideration and D is the maximum number of steepest descent steps in attraction basins associated with pL.

Computing the partition function for kinetically trapped RNA secondary structures

An RNA secondary structure is locally optimal if there is no lower energy structure that can be obtained by the addition or removal of a single base pair, where energy is defined according to the widely accepted Turner nearest neighbor model. Locally optimal structures form kinetic traps, since any evolution away from a locally optimal structure must involve energetically unfavorable folding steps. Here, we present a novel, efficient algorithm to compute the partition function over all locally optimal secondary structures of a given RNA sequence. Our software, RNAlocopt runs in time and space. Additionally, RNAlocopt samples a user-specified number of structures from the Boltzmann subensemble of all locally optimal structures. We apply RNAlocopt to show that (1) the number of locally optimal structures is far fewer than the total number of structures – indeed, the number of locally optimal structures approximately equal to the square root of the number of all structures, (2) the structural diversity of this subensemble may be either similar to or quite different from the structural diversity of the entire Boltzmann ensemble, a situation that depends on the type of input RNA, (3) the (modified) maximum expected accuracy structure, computed by taking into account base pairing frequencies of locally optimal structures, is a more accurate prediction of the native structure than other current thermodynamics-based methods. The software RNAlocopt constitutes a technical breakthrough in our study of the folding landscape for RNA secondary structures. For the first time, locally optimal structures (kinetic traps in the Turner energy model) can be rapidly generated for long RNA sequences, previously impossible with methods that involved exhaustive enumeration. Use of locally optimal structure leads to state-of-the-art secondary structure prediction, as benchmarked against methods involving the computation of minimum free energy and of maximum expected accuracy. Web server and source code available at http://bioinformatics.bc.edu/clotelab/RNAlocopt/.

Efficient exploration of discrete energy landscapes

Many physical and chemical processes, such as folding of biopolymers, are best described as dynamics on large combinatorial energy landscapes. A concise approximate description of the dynamics is obtained by partitioning the microstates of the landscape into macrostates. Since most landscapes of interest are not tractable analytically, the probabilities of transitions between macrostates need to be extracted numerically from the microscopic ones, typically by full enumeration of the state space or approximations using the Arrhenius law. Here, we propose to approximate transition probabilities by a Markov chain Monte Carlo method. For landscapes of the number partitioning problem and an RNA switch molecule, we show that the method allows for accurate probability estimates with significantly reduced computational cost.