A set of nearest neighbor parameters for predicting the enthalpy chang
http://www.100md.com
《核酸研究医学期刊》
Department of Biochemistry & Biophysics and Center for Pediatric Biomedical Research, University of Rochester Medical Center 601 Elmwood Avenue, Box 712, Rochester, NY 14642, USA 1 Department of Chemistry, University of Rochester Box 0216, Rochester, NY 14627, USA
*To whom correspondence should be addressed. Tel: 1 585 275 1734; Fax: 1 585 275 6007; Email: david_mathews@urmc.rochester.edu
ABSTRACT
A complete set of nearest neighbor parameters to predict the enthalpy change of RNA secondary structure formation was derived. These parameters can be used with available free energy nearest neighbor parameters to extend the secondary structure prediction of RNA sequences to temperatures other than 37°C. The parameters were tested by predicting the secondary structures of sequences with known secondary structure that are from organisms with known optimal growth temperatures. Compared with the previous set of enthalpy nearest neighbor parameters, the sensitivity of base pair prediction improved from 65.2 to 68.9% at optimal growth temperatures ranging from 10 to 60°C. Base pair probabilities were predicted with a partition function and the positive predictive value of structure prediction is 90.4% when considering the base pairs in the lowest free energy structure with pairing probability of 0.99 or above. Moreover, a strong correlation is found between the predicted melting temperatures of RNA sequences and the optimal growth temperatures of the host organism. This indicates that organisms that live at higher temperatures have evolved RNA sequences with higher melting temperatures.
INTRODUCTION
RNA is more than a simple single-stranded sequence carrying genetic information as in the Central Dogma of Biology. For example, it can form tertiary structures that, such as proteins, can be catalytic. Natural and engineered RNA molecules are widely used as functional tools in enzymatic catalysis and genetic control (1–5). One current problem is how to predict the structures of functional RNA sequences.
Secondary structure, the sum of canonical base pairs, is stronger (6–9) and forms faster (10) than tertiary structure. Therefore, secondary structure can largely be determined without knowledge of tertiary structure. Comparative sequence analysis is a standard technique for determining the secondary structure of homologous RNA sequences (11–13). When only a few or even a single sequence is available, the secondary structure at 37°C can be predicted by free energy minimization algorithms (14–17) using a set of empirical free energy parameters, determined from optical melting experiments (17–21). Each parameter only depends on the sequence identity of nucleotides in the motif and in adjacent base pairs and the total free energy is the sum of nearest neighbor terms. The average sensitivity (the percentage of known base pairs that are correctly predicted) of free energy minimization prediction has been benchmarked as high as 72.8 ± 9.4% for a diverse database of sequences having fewer than 800 nt (17). Furthermore, experimentally determined constraints can improve this accuracy of prediction up to 84% (17,18) for sequences with <6% pseudoknotted (non-nested) base pairs (17). Partition function prediction of base pair probabilities can be used to identify base pairs in the predicted lowest free energy structure that are much more likely than average to be in the known secondary structure (22,23). For example, 91.0% of base pairs in the lowest free energy structure with pairing probability of 0.99 or higher are contained in the known structure, on average (22). The high accuracy of thermodynamic structure prediction (17) demonstrates that many RNA secondary structures can be determined from sequences, without knowledge of any tertiary contacts or protein interactions.
The current set of free energy nearest neighbor parameters for predicting the free energy of RNA secondary structure, however, is limited to application at 37°C. Many organisms, thermophiles and psychrophiles, live at temperatures far from 37°C and many experiments are conducted at other temperatures. The prediction of secondary structure of RNA at arbitrary temperature would expand our knowledge of structure and evolution in the RNA world. Moreover, it would facilitate studying and designing functional RNA molecules at temperatures other than 37°C. The enthalpy nearest neighbor parameters can be used in conjunction with available free energy nearest neighbor parameters for 37°C to determine free energy nearest neighbors at other temperatures. But the most recent enthalpy parameters were derived in 1995 using a simple model (24). At that time, no themes had emerged for the sequence-dependent stability of internal loops. Subsequently, the nearest neighbor model for free energy change at 37°C was significantly improved (17) using experimental results. Therefore, we applied the principles of the current free energy nearest neighbor model (17,18) to determine a complete set of enthalpy nearest neighbor parameters using the available optical melting data.
MATERIALS AND METHODS
Database of experiments
The database of experimental data for derivation of enthalpy parameters is included in Supplementary Data. It includes 130 hairpin loops (25–31), 37 bulge loops (32,33), 337 internal loops (17,18,34–49) (99 of which are 2 x 2 internal loops), 74 multibranch loops (50,51) and 43 coaxial stacking models (52–55).
Derivation and refinement of enthalpy parameters
Canonical base pairs
The enthalpies of Watson–Crick and GU base pairs were derived by Xia et al. (21) and Mathews et al. (18), respectively.
Dangling ends and terminal mismatches
Dangling ends are unpaired nucleotides adjacent to canonical pairs and their enthalpy parameters were compiled previously (24). Dangling ends on terminal GU pairs are treated similar to dangling ends on terminal AU pairs. Terminal mismatches are non-canonical pairs at the end of helixes. The enthalpy parameters of terminal mismatches are taken from another compilation (20), with the exception of mismatches on terminal GU pairs, which were measured recently (30).
If a terminal mismatch has the potential to pair canonically, the values of A–C and C–A mismatches are used for the purine–pyrimidine mismatch and pyrimidine–purine mismatches, respectively. This is important for partition function calculations, where all possible secondary structures are considered.
Hairpin loops
The experimental enthalpies of hairpin loop formation are calculated from published experimental data (25–31) with the following equation:
where is the experimental value for unfolding the hairpin loop with stem, is calculated by the INN-HB parameters (18,21), without an intermolecular initiation term.
The hairpin loop enthalpy parameters are estimated by linear regression using the same model as free energy nearest neighbor parameters (17), except that the GG first mismatch bonus observed for free energy does not apply for enthalpy because the bonus was not statistically significant for enthalpy. The GG stability bonus is therefore entropic in nature, consistent with the observation that GG mismatches are dynamic (56), i.e. they sample more than one single microstate on short timescales.
The enthalpies of hairpin loops are estimated by the following equation:
where n is the number of unpaired nucleotides in the loop. Hairpins with fewer than 3 unpaired nucleotides are not allowed by the model. When n = 3, only the initiation term is considered without any bonus and penalty terms, except a penalty for hairpin loops with three Cs. When n > 3, the special GU closure bonus applies to GU closed hairpins in which a 5' closing G is preceded by two G residues; and (UU or GA first mismatch but not AG) is applied to loops with first mismatches of UU or GA (G on the 5' side and A on 3' side of loop). The oligo-C penalty applies only to loops composed of all C residues and, if n > 3, is calculated with (oligo-C loops, n > 3) = An + B. For hairpinloops composed entirely of 3 C residues, the (oligo-C loops, n = 3) is applied.
The enthalpy parameters are listed in Table 1 and the database of measured loop enthalpies is available as Supplementary Data. In the absence of data, for hairpin loops longer than 9 nt, the initiation enthalpy is approximated with the initiation term for a hairpin of 9 nt. This assumes that additional instability of hairpin loops as the loop lengthens derives from the entropy (57).
Table 1 Hairpin loop enthalpy parametersa
The measured free energies at 37°C of some special hairpin loops of 3, 4 or 6 unpaired nucleotides (30,31,34–36) are either more or less stable by 0.9 kcal/mol than the model predicts. The enthalpies for each of these sequences are listed in a separate lookup table (Table 2), to be consistent with the free energy parameters.
Table 2 Lookup table for unstable triloops and stable tetraloops and hexaloops
Bulge loops
RNA secondary structure is destabilized by bulge loops, which are an interruption of helical structure in one strand only (32,37,38). The initiation terms, for bulge loops of 1–3 nt, are listed in Table 3. They are the average values of experimental data (32,33), calculated using the following equation:
where the enthalpy of the duplex without bulge is the experimental value of the sequence of the duplex without the bulge or as calculated with INN-HB parameters (21) if the experimental values were not available. is the stacking enthalpy of the base pairs in the duplex without the bulge that flank the bulge loop in the duplex with the bulge. Because the difference of initiation enthalpies between 2 and 3 nt bulges is almost zero, it is assumed that the increasing instability for longer bulges (n 4) comes from the entropy of the loop closure (39,57). Thus, the initiation enthalpy for bulges longer than 3 nt is approximated as the 3 nt bulge enthalpy.
Table 3 Bulge loop initiation enthalpy parametersa
Assuming that helical stacking is continuous between the adjacent helices for single bulges, but is interrupted by longer bulges (39,40), the enthalpies of bulge loops are calculated with the following equation:
The calculation of enthalpies for the adjacent helices would include the terminal AU/GU penalty (21) for AU/GU pairs adjacent to the bulge loops that are longer than 1 nt. is the canonical helix stacking enthalpy applied for the two closing base pairs as though the helix was not interrupted by the bulge loop.
Internal loops
Internal loop enthalpies were calculated from experimental data (17,18,34–49) using the following equation:
The range of measured enthalpies differs for internal loops of different size and symmetry; therefore, different enthalpy models are used to predict different loop types. The models are similar to those used to model free energies (17).
1 x 1 Internal loops (single mismatches)
For single non-canonical pairs (1 x 1 internal loops), the loop enthalpies are approximated by the following equation:
where is the enthalpy of initiation for a single non-canonical pair; is the penalty for each AU or GU closing base pair; is a bonus for a GG pair in a 1 x 1 loop; and is a bonus for a 5'RU/3'YU stack in a 1 x 1 loop, where R is a purine and Y is a pyrimidine.
2 x 2 Internal loops (tandem mismatches)
The 2 x 2 internal loops, also called tandem mismatches, interrupt helical RNA with two opposing unpaired nucleotides on each strand. Many of the sequence-symmetric 2 x 2 loops have been studied experimentally (17,18, 34–49) and their enthalpies are assembled in a ‘periodic table’ (Table 4). Symmetric sequences that have not been measured are approximated by averaging the most adjacent columns that have been measured. For asymmetric 2 x 2 loops, the enthalpies are approximated using the following equation:
where (12.5 ± 2.7 kcal/mol) is applied to loops with a GG pair adjacent to an AA or any non-canonical pair with a pyrimidine and p (2.4 ± 3.1 kcal/mol) is applied to loops with an AG or GA pair adjacent to a UC, CU or CC pair or with a UU pair adjacent to an AA pair.
Table 4 The periodic table of tandem mismatch (2 x 2 internal loop) enthalpya
Other internal loops
The enthalpies of other internal loops are approximated using the following equation:
where is the enthalpy of initiation for a loop of n nucleotides; is a penalty for loops with unequal numbers of nucleotides on each side, with n1 and n2 the number of nucleotides on each side; is applied for each sequence-specific first mismatch (Table 5), but it is not applied to loops of the form 1 x (n – 1) with n > 3 (n is the total number of unpaired bases). Special first mismatch bonuses were determined for 2 x 3 and 1 x 2 internal loops with separate linear regressions.
Table 5 Approximations for internal loop enthalpy parameters at 37°C (in kcal/mol)
Moreover, the free energy parameters (Table 6) were updated for internal loops based on recent experimental measurements. The free energy parameters were obtained using the method of Mathews et al. (17). The recent data include the 3 x 3 loops from Chen et al. (41), but excluding the 3 x 3 loops with a middle GA pair. The middle GA pair is shown to enhance stability and this extra stability cannot be predicted by the nearest neighbor parameter set used in this work (41).
Table 6 Updated internal loop free energy parameters at 37°C (in kcal/mol)
Coaxial stacking
Coaxial stacking, which is a favorable interaction of two helices stacked end to end, occurs in multibranch loops and exterior loops. Stability increments for coaxial stacking were measured with a structure composed of a short oligonucleotide bound to a single-stranded end of a stem–loop structure, creating a helical interface (52–55). The enthalpy of coaxial stacking is quantified as follows:
where H°(correction) is the enthalpy for displacing a 3' dangling end on the stem–loop structure if one is present.
When the helixes have no intervening mismatches, the enthalpy bonus is approximated by the nearest neighbor parameter (21) of a base pair in a helix. The excess enthalpy above the helical stacking nearest neighbor from Xia et al. (21), , for each measured interface was calculated. With flush interfaces, i.e. with no intervening mismatch, and no strand extensions beyond the interface, the average excess enthalpy is –1.53 ± 1.45 kcal/mol. For interfaces followed by strand extensions, the excess enthalpy is 1.82 ± 1.13 kcal/mol. As the excess enthalpy changes are not statistically significant, coaxial stacking of helices with no intervening nucleotides is modeled with the enthalpy parameter in a helix.
With one intervening nucleotide from each strand, two helices can stack with an intervening mismatch between them. There are two stack increments: one is the mismatch stack at the end of one helix with continuous backbone, which is equal to the mismatch stacking parameter on a helix, and the other is the mismatch stack with discontinuous backbone, which is modeled as sequence independent. The average enthalpy of sequence independent stacks is –8.46 ± 2.75 kcal/mol. In addition to this, an enthalpy bonus of –0.4 or –0.2 kcal/mol are applied to intervening mismatches composed of nucleotides that could form a Watson–Crick or a GU base pair, respectively. These bonuses are identical to free energy increments that are used and are empirically found to improve structure prediction accuracy.
Multibranch loops
The parameters are determined by linear regression of experimental data for three- and four-way multibranch loops (50,51). In a nearest neighbor model, the bimolecular enthalpy () for the formation of the duplex with a multibranch loop is given by the following equation:
where helix 1 and helix 2 are the intermolecular paired helices with H° predicted from nearest neighbor parameters for Watson–Crick pairs (without including bimolecular initiation so that appears only once). The is a term that accounts for the stacking enthalpy increment of the nucleotides that can stack on the hairpin loop stems to form a modified motif after the two strands have dissociated. This is the most favorable configuration with coaxial stacking of helixes (in the case of four-way multibranch loops) or of the stacking of unpaired nucleotides. is the experimental value which is taken from versus ln(CT/4) plots. The multibranch loop enthalpy initiation term () can be calculated from the above equation. The enthalpy of multibranch loops () is then modeled as the sum of two terms, initiation and stacking:
The stacking term is the favorable enthalpy of coaxial stacking, terminal mismatch and/or dangling end stacking. It is determined from the stacking conformation that gives the lowest free energy, as determined by free energy nearest neighbors (50). The initiation term can be approximated by the following equation:
where a, b and c are parameters determined from linear regression (Table 7) and h is the number of branching helices. is a strain enthalpy that only applies to three-way multibranch loops with fewer than two unpaired nucleotides. The asym term is the average asymmetry that reflects the distribution of unpaired nucleotides, which is defined by the following equation:
Table 7 Enthalpy parameters for multibranch loop initiation
The average asymmetry is limited to 2.0, following the rules suggested by free energy parameters. Asymmetry cannot be applied, however, by dynamic programming algorithms for secondary structure prediction (17,22). Thus, the b term was excluded for secondary structure prediction and the parameters a and c were optimized by finding the parameters that lead to the highest average sensitivity of secondary structure prediction by free energy minimization. The maximum sensitivity of prediction was found with a = 30.0 kcal/mol and c = –2.2 kcal/mol.
Database of RNA secondary structures
The revised enthalpy nearest neighbor model was tested with RNA sequences with known secondary structure from organisms with known optimal growth temperature. The structures were taken from comparative analysis databases (42–49,58,59). Small (16S) subunit rRNA sequences are divided into domains as defined by Jaeger et al. (39). Large (23S) subunit rRNA sequences are divided into domains of fewer than 700 nt each (18). The optimal growth temperatures of different organisms were taken from the Prokaryotic Growth Temperature Database (http://pgtdb.csie.ncu.edu.tw/) and the DSMZ German Collection of Microorganisms and Cell Cultures website (http://www.dsmz.de/). Only the RNA sequences of mesophiles (organisms living at temperatures between 10 and 60°C, but with organisms living at 37°C excluded) were chosen to test the sensitivity and positive predictive value (PPV) of secondary structure prediction. Considering that posttranscriptional modification (60) and high pressure (61) in the thermophiles and hyperthermophiles (organism living above 60°C) would change the thermodynamics of secondary structure formation, sequences from these organisms were excluded. A list of sequences and optimal growth temperatures used are available in Supplementary Data.
Accuracy of secondary structure prediction
The accuracy of structure prediction is determined by the sum of the canonical base pairs correctly predicted. A base pair is considered correctly predicted even if it is shifted by 1 nt on one side. For example a base pair between nucleotides i and j is considered to be correctly predicted if any of these base pairs is predicted: i to j, i to j – 1, i to j + 1, i – 1 to j or i + 1 to j. The predicted base pair between i – 1 and j + 1, however, is not considered to be correct. This scoring scheme reflects the uncertainty of exact base pair matches in comparative sequence analysis and the possibility for dynamics in base pairing. The values of sensitivity and PPV of this scoring scheme are 2–3% higher than when determined with exact base pairing only, where only the i to j base pair is considered to be correct. The prediction accuracies are shown in Supplementary Tables 11 and 12. Each table includes accuracies determined when pairs can be shifted and when pairs must be an exact match.
Availability of parameters
Machine-readable tables of the enthalpy parameters are available on the Mathews lab website (http://rna.urmc.rochester.edu/).
RESULTS
Nearest neighbor model parameters
In the nearest neighbor model of free energy (17,18), the parameters for Watson–Crick base pairs are well determined at 37°C with errors <10%, or 0.1–0.2 kcal/mol (21). For other motifs such as loops and GU base pairs, individual nearest neighbor free energy increments are often determined with an error <0.5 kcal/mol (17,18). In order to extend the current model to predict free energy at temperatures other than 37°C, enthalpy parameters consistent with the current nearest neighbor model are required. The free energy at arbitrary temperature for each parameter is then
(1)
where the enthalpy (H°) and entropy (S°) are assumed to be temperature independent. As described in Materials and Methods, parameters for enthalpy prediction, compatible with the free energy model, were determined using available experimental data from optical melting experiments.
Experimental studies consistently demonstrate that enthalpy and entropy measurements have considerably larger percent error than free energy measurements. Free energy at 37°C is determined with greater precision because of correlation between errors in enthalpy and entropy (21). The larger experimental errors in enthalpy result in larger percent errors for enthalpy nearest neighbor parameters than free energy parameters. The enthalpy of RNA secondary structure is known to be a function of temperature. A linear model for heat capacity change predicts the following:
(2)
(3)
where is a constant heat capacity change and T0 is a chosen reference temperature. It is hypothesized that the heat capacity change arises from the extent of stacking increasing with decreasing temperature. Thus, is negative because single strands are more organized at low rather than high temperature (62–67). The can be estimated by linear fits of enthalpy and entropy changes as a function of melting temperature (50,51,62) or determined by isothermal titration calorimetry at multiple temperatures (68,69). However, the effects of heat capacity change on enthalpy and entropy are antagonistic in terms of free energy change:
(4)
Therefore, for certain T (T = T – T0), can be neglected because the effects are compensated in terms of free energy. To calculate the compensation for a set of RNA duplexes (62), the free energy, G°, was derived directly from Equation 4 assuming that the entropy and enthalpy were independent of temperature. Then the temperature-dependent free energy, , was calculated with the measured non-zero from Equations 2–4. The free energy difference, G° = – G°, increases with the deviation of temperature from T0 (37°C) (Figure 1). The exact G° for each duplex is shown in Table 8 for different temperatures. The experimental error in individual loop free energy nearest neighbor parameters at 37°C is as large as 0.5 kcal/mol (17), which corresponds to roughly a factor of 2 in equilibrium constant. Thus, the small G° for helices suggests that the approximation of = 0 is reasonable for predictions from 10 to 60°C. Therefore, the enthalpy parameters derived here assume = 0 and are most accurate at predicting free energy change close to 37°C.
Figure 1 (A) Free energy difference of RNA duplex CCGGUp. G° (dashed line) was derived from Equation 3, where enthalpy and entropy were averaged from the optical melting curve fits, assuming that they were independent of the temperature. (solid line) was calculated from Equations 1–3, where the heat capacity was accounted. (B) Free energy difference is G° = (62).
Table 8 Free energy differences of RNA duplexes
Dynamic programming algorithm for RNA secondary structure prediction
RNAstructure is a program for RNA secondary structure prediction and analysis. It includes prediction of secondary structure by free energy minimization (17), prediction of base pair probabilities using a partition function (22), the efn2 function for predicting the free energy change of folding given a sequence and secondary structure (18), and the Dynalign algorithm for finding the secondary structure common to two sequences (70). RNAstructure was revised to make predictions at user-defined temperature. Because large internal loops are more likely at high temperature, the previous limitation on internal loop size (fewer than 30 unpaired nucleotides) (17,18,22) was removed by implementing the method of Lyngs? et al. (71). This provides an O(N3) algorithm that can predict internal loops of arbitrary size. Benchmarks for calculation time and memory requirement with and without this revision are shown in Table 9.
Table 9 Calculation time and memory size of dynamic programming for sequences of different length
Sensitivities and PPVs of structure predictions
The enthalpy nearest neighbor parameters were compared with the previous parameters and model for enthalpy and free energy assembled by Serra and Turner (24) by predicting the secondary structures of RNA sequences with known secondary structures. Sensitivities, the percent of known base pairs that are correctly predicted, using both sets of parameters are shown in Figure 2 (detailed numbers are in Supplementary Table 11A) for different types of structural RNA sequences. The known structures of these sequences were taken from comparative analysis databases (42–49,58,59). The average sensitivity is improved from 65.2 to 68.9% using the new parameters assembled here. Sensitivities are improved for most types of the RNA. The exceptions are 5S rRNA and Group II introns.
Figure 2 Improvement of prediction at optimal growth temperatures. The sequences are those from mesophiles (optimal growth temperature from 10 to 60°C) without organisms with optimal growth at 37°C. The lowest free energy secondary structures were predicted at the organims' optimal growth temperatures using two models. The previous model and parameters are those of Serra and Turner (24), which are widely used. The improved prediction uses the model and parameters presented in this work. The small and large subunits of rRNA sequences are divided into domains of <700 nt. The total sensitivity is the average of sensitivities of different types of RNA.
To test the enthalpy parameters, the accuracy of secondary structure prediction at optimal growth temperature was compared to the accuracy of structure prediction at 37°C for organisms that do not grow optimally at 37°C for several types of RNAs (Table 10). The comparison of predictions was shown in different groups divided by optimal growth temperature. The organisms in each group grow optimally in a certain range of temperatures. Compared to the prediction at 37°C, structure prediction at optimal growth temperature performs better for the organism living at temperatures between 22 and 37°C, but is worse at other optimal growth temperatures. This suggests that when enthalpy parameters are assumed to be temperature independent, their utility as a tool for deriving free energy parameters for use in predicting the lowest free energy structure is limited to a narrow temperature range. Small errors in enthalpy change parameters have a larger effect on free energy change parameter determination (Equation 1), the farther the temperature is from 37°C.
Table 10 Prediction sensitivities of the lowest free energy structurea
Figure 3 shows the PPV for base pairs from the lowest free energy structure for base pairs with different pairing probabilities (see detailed numbers in Supplementary Table 12A). They are predicted using a partition function calculation at optimal growth temperature (22). PPV is the percentage of predicted base pairs that are found in the known structure. The average PPV of all pairs in the lowest free energy structures is only 62.0%, which is lower than the sensitivity (68.9%). This suggests that the model over-predicts base pairs and/or that the base pairs may not be annotated completely in the structures from comparative analysis (22). For example, if a base pair is completely conserved, then it is sometimes not annotated by comparative analysis (42–49,58,59). Base pair probabilities for all possible pairs are calculated with a partition function and grouped by different thresholds. The PPV is significantly higher for predicted base pairs in the lowest free energy structure with higher pairing probability. The average PPV is up to 90.4% for those known base pairs having probability of 0.99 or above. It has been demonstrated previously that base pair probabilities predicted at 37°C can be used to find pairs with high PPV (22). The fact that this holds true at other temperatures shows that the enthalpy parameters are robust for base pair probability prediction.
Figure 3 PPV for optimal structure and base pairs with different pairing probabilities. PPV equals the number predicted base pairs in that are in the known structure divided by total number of predicted base pairs. Pairs in the optimal structures are grouped by different thresholds of pairing probabilities. The pairing probabilities were calculated with a partition function calculation (22) at organisms' optimal growth temperatures, using the model and parameters presented in Materials and Methods. The small and large subunits of rRNA sequences are divided into domains of <700 nt. The sequences of different type of RNA are those from mesophiles (living from 10 to 60°C) without organisms living at 37°C.
The fact that the accuracy of secondary structure prediction is sensitive to the accuracy of the nearest neighbor parameters, but the base pair probabilities remain a robust measure of confidence for a wide variety of temperatures is consistent with a previous work. Layton and Bundschuh (72) demonstrated that the predicted lowest free energy structure was often changed in repeated structure predictions after random adjustments of the nearest neighbor parameters within the limits of their error. Base pair probabilities, however, were less perturbed by changes in the parameters (72). With the extrapolation of nearest neighbor parameters to temperatures far from 37°C, the accuracy of the predicted lowest free energy structure is often reduced as compared to structure prediction at 37°C. The ability of the partition function predicted base pair probabilities to determine base pairs predicted with a higher confidence is unchanged with secondary structure prediction at temperatures far from 37°C. This is because the determination of base pair probabilities is not as perturbed by errors in the nearest neighbor parameters.
An example of secondary structure prediction at 37°C and at optimal growth temperature of 30°C is shown in Figure 4 for a tRNA sequence. The base pairs with higher predicted pairing probability (color annotated according to pairing probability in Figure 4B and C) are pairs predicted with greater confidence. For this sequence, secondary structure prediction is more accurate and the fidelity of structure prediction (as judged by the percent of high probability pairs) is improved at optimal growth temperature.
Figure 4 Secondary structure prediction of Saccharomyces cerevisiae tRNA (RM4000) at optimal growth temperature (30°C) (B) and at 37°C (C) with the presented nearest neighbor parameters. Base pairs in the original structure (A) are derived from the comparative analysis database (42–49,58,59). Structures are also color annotated to indicate predicted base pair probabilities (Pbp) for each helix: red, Pbp 0.95; yellow, 0.95 > Pbp 0.7; green, 0.7 > Pbp 0.3; blue, 0.3 > Pbp. The structures were drawn with XRNA (http://rna.ucsc.edu/rnacenter/xrna/xrna.html) and Adobe Illustrator.
Correlation between melting temperature and optimal growth temperature
Melting temperature, Tm, is defined as the temperature at which half of strands are unpaired. Assuming that an RNA melts with a two-state transition, the melting temperature (in Kelvins) of a single-stranded RNA structure can be predicted by Tm = H°/S° (73). For example, the predicted melting temperatures (°C) for all hairpins in the database of optically melted sequences (Supplementary Data) (25–31) are plotted in Figure 5 as a function of experimentally determined Tm. This shows that the parameters adequately reflect the thermal stabilities of RNA sequences with known Tm. Better correlation was found at higher temperatures. This is expected because most hairpins were measured with high melting temperatures in experiments (25–31).
Figure 5 Experimental (Supplementary Data) (25–31) versus predicted (Tm = H°/S° – 273.15) melting temperatures of hairpin stem–loop structures. The line shows the ideal location of points, predicted Tm = measured Tm. The root mean squared deviation (r.m.s.d.) of prediction compared to experiment is 5.86°C. The new enthalpy parameters provide improved Tm prediction compared to the previous compilation of parameters (24), which have an r.m.s.d. of 7.58°C as compared to experiment for this dataset.
Melting temperature reflects the thermal stability of a structure. Therefore RNA structures in organisms living at higher temperature are expected to have higher melting temperatures. Figure 6A shows a plot of predicted melting temperatures of the lowest free energy structure versus organism optimal growth temperature (10–90°C). A strong correlation (linear correlation coefficient of 0.797) is found between the melting temperature and the optimal growth temperature for different types of RNA structures. On the other hand, there appears to be less correlation between nucleotide content and optimal growth temperature (Figure 6B–D) for diverse types of RNA, although uracil content of 16S rRNA of thermophiles and psychrophiles were found recently to correlate inversely with their optimal growth temperatures (74). Evidently, the thermal stability of RNA structure is not simply controlled by base content. Organisms that grow at high temperature have apparently evolved RNA secondary structures with a combination of motifs that provide thermal stability.
Figure 6 Relationships of melting temperatures, nucleotide contents and optimal growth temperatures of different types of RNA in different organisms with optimal growth temperature from 10 to 90°C: (A) Predicted melting temperature; (B) G–C pair content; (C) G content; and (D) U content versus optimal growth temperature. Melting temperatures are predicted for different types of RNA sequences from comparative analysis databases (42–49,58,59) with a two-state transition assumption.
DISCUSSION
The nearest neighbor parameters for enthalpy were derived here using similar rules as for free energy nearest neighbor parameters at 37°C (17). This makes these parameters useful for determining free energy parameters at arbitrary temperature that are compatible with dynamic programming algorithms for secondary structure prediction. Some of the enthalpy parameters have large percent standard errors as compared with the parameters of free energy. This reflects the larger errors in the experimental results of enthalpy than free energy, but it also suggests that enthalpy may be more sequence dependent than free energy. This sequence dependence cannot be determined using the currently available database of optical melting experiments and suggests a need for further optical melting experiments on model RNA systems.
Another source of error comes from the assumption that the enthalpy and entropy are independent of the temperature in both the model and in the analysis of optical melting experiments. When the temperature is too far from 37°C, the sensitivity of prediction is expected to be worse than 68.9% on average because of the approximation of . For example, experiments demonstrate cold denaturation of RNA (68,69), but the nearest neighbor model does not reproduce those results. Further experiments by isothermal titration calorimetry would be needed to provide the data for a model that can include a non-zero heat capacity change.
There are common error sources that should be considered for the prediction of base pairs. Free energy minimization assumes that the secondary structure is at equilibrium. The nearest neighbor model is an incomplete representation of structural free energy. The parameters average some sequence-specific effects and were derived from a limited set of experiments. Some RNA sequences, in particular mRNA, may sample multiple structures at equilibrium. The parameters are derived from experimental data at 1 M NaCl, whereas the salt concentration in different organisms may be very different.
In spite of all these limitations, the nearest neighbor model predicts secondary structures with a 72.8% average sensitivity (17). Recent experimental results on the self-folding of the 16S rRNA 5' domain (75) support the assumption of thermodynamic control of folding pathway. Moreover, the base pair prediction with the partition function can be used to determine pairs predicted with greater confidence (22).
In spite of the fact that the enthalpy parameters have larger percent errors than the free energy parameters for 37°C, the enthalpy parameters are able to predict optical melting temperatures for small model sequences. Predicted melting temperatures for structural RNA sequences correlate well with optimal growth temperature, suggesting that these parameters capture many of the sequence-dependent features of RNA folding enthalpy change.
SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.
ACKNOWLEDGEMENTS
The authors thank Rahul Tyagi and Andrew V. Uzilov for helpful discussions. D.H.M. is an Alfred P. Sloan Research Fellow. This work was supported by National Institutes of Health Grants GM22939 to D.H.T. and GM076485 to D.H.M. Funding to pay the Open Access publication charges for this article was provided by National Institutes of Health.
REFERENCES
Nelson, P., Kiriakidou, M., Sharma, A., Maniataki, E., Mourelatos, Z. (2003) The microRNA world: small is mighty Trends Biochem. Sci, . 28, 534–540 .
Doudna, J. and Cech, T. (2002) The chemical repertoire of natural ribozymes Nature, 418, 222–228 .
Walter, P. and Blobel, G. (1982) Signal recognition particle contains a 7S RNA essential for protein translocation across the endoplasmic reticulum Nature, 299, 691–698 .
Lau, N.C., Lim, L.P., Weinstein, E.G., Bartel, D.P. (2001) An abundant class of tiny RNAs with probable regulatory roles in Caenorhabditis elegans Science, 294, 858–862 .
Lagos-Quintana, M., Rauhut, R., Lendeckel, W., Tuschl, T. (2001) Identification of novel genes coding for small expressed RNAs Science, 294, 853–858 .
Onoa, B. and Tinoco, I., Jr. (2004) RNA folding and unfolding Curr. Opin. Struct. Biol, . 14, 374–379 .
Crothers, D.M., Cole, P.E., Hilbers, C.W., Schulman, R.G. (1974) The molecular mechanism of thermal unfolding of Escherichia coli formylmethionine transfer RNA J. Mol. Biol, . 87, 63–88 .
Mathews, D.H., Banerjee, A.R., Luan, D.D., Eickbush, T.H., Turner, D.H. (1997) Secondary structure model of the RNA recognized by the reverse transcriptase from the R2 retrotransposable element RNA, 3, 1–16 .
Banerjee, A.R., Jaeger, J.A., Turner, D.H. (1993) Thermal unfolding of a group I ribozyme: the low temperature transition is primarily a disruption of tertiary structure Biochemistry, 32, 153–163 .
Woodson, S.A. (2000) Recent insights on RNA folding mechanisms from catalytic RNA Cell Mol. Life Sci, . 57, 796–808 .
James, B.D., Olsen, G.J., Pace, N.R. (1989) Phylogenetic comparative analysis of RNA secondary structure Methods Enzymol, . 180, 227–239 .
Pace, N.R., Thomas, B.C., Woese, C.R. (1999) Probing RNA structure, function, and history by comparative analysis In Gesteland, R.F., Cech, T.R., Atkins, J.F. (Eds.). The RNA World, 2nd edn Cold Spring Harbor Cold Spring Harbor Laboratory Press pp. 113–141 .
Woese, C.R., Gutell, R.R., Gupta, R., Noller, H.F. (1983) Detailed analysis of the higher order structure of 16S-like ribosomal ribonucleic acids Microbiol. Rev, . 47, 621–669 .
Andronescu, M., Aguirre-Hernandez, R., Condon, A., Hoos, H.H. (2003) RNAsoft: A suite of RNA secondary structure prediction and design software tools Nucleic Acids Res, . 31, 3416–3422 .
Hofacker, I.L. (2003) Vienna RNA secondary structure server Nucleic Acids Res, . 31, 3429–3431 .
Zuker, M. (2003) Mfold web server for nucleic acid folding and hybridization prediction Nucleic Acids Res, . 31, 3406–3415 .
Mathews, D.H., Disney, M.D., Childs, J.L., Schroeder, S.J., Zuker, M., Turner, D.H. (2004) Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure Proc. Natl Acad. Sci. USA, 101, 7287–7292 .
Mathews, D.H., Sabina, J., Zuker, M., Turner, D.H. (1999) Expanded sequence dependence of thermodynamic parameters provides improved prediction of RNA Secondary Structure J. Mol. Biol, . 288, 911–940 .
Turner, D.H. (2000) Conformational changes In Bloomfield, V., Crothers, D., Tinoco, I. (Eds.). Nucleic Acids, Sausalito, CA University Science Books pp. 259–334 .
Xia, T., Mathews, D.H., Turner, D.H. (1999) Thermodynamics of RNA secondary structure formation In Soll, D.G., Nishimura, S., Moore, P.B. (Eds.). Prebiotic Chemistry, Molecular Fossils, Nucleosides, and RNA, NY Elsevier pp. 21–47 .
Xia, T., SantaLucia, J., Jr, Burkard, M.E., Kierzek, R., Schroeder, S.J., Jiao, X., Cox, C., Turner, D.H. (1998) Thermodynamic parameters for an expanded nearest-neighbor model for formation of RNA duplexes with Watson–Crick pairs Biochemistry, 37, 14719–14735 .
Mathews, D.H. (2004) Using an RNA secondary structure partition function to determine confidence in base pairs predicted by free energy minimization RNA, 10, 1178–1190 .
McCaskill, J.S. (1990) The equilibrium partition function and base pair probabilities for RNA secondary structure Biopolymers, 29, 1105–1119 .
Serra, M.J. and Turner, D.H. (1995) Predicting thermodynamic properties of RNA Methods Enzymol, . 259, 242–261 .
Giese, M.R., Betschart, K., Dale, T., Riley, C.K., Rowan, C., Sprouse, K.J., Serra, M.J. (1998) Stability of RNA hairpins closed by wobble base pairs Biochemistry, 37, 1094–1100 .
Serra, M.J., Lyttle, M.H., Axenson, T.J., Schadt, C.A., Turner, D.H. (1993) RNA hairpin loop stability depends on closing pair Nucleic Acids Res, . 21, 3845–3849 .
Serra, M.J., Axenson, T.J., Turner, D.H. (1994) A model for the stabilities of RNA hairpins based on a study of the sequence dependence of stability for hairpins of six nucleotides Biochemistry, 33, 14289–14296 .
Serra, M.J., Barnes, T.W., Betschart, K., Gutierrez, M.J., Sprouse, K.J., Riley, C.K., Stewart, L., Temel, R.E. (1997) Improved parameters for the prediction of RNA hairpin stability Biochemistry, 36, 4844–4851 .
Groebe, D.R. and Uhlenbeck, O.C. (1988) Characterization of RNA hairpin loop stability Nucleic Acids Res, . 16, 11725–11735 .
Dale, T., Smith, R., Serra, M.J. (2000) A test of the model to predict unusually stable RNA hairpin loop stability RNA, 6, 608–615 .
Antao, V.P. and Tinoco, I., Jr. (1992) Thermodynamic parameters for loop formation in RNA and DNA hairpin tetraloops Nucleic Acids Res, . 20, 819–824 .
Longfellow, C.E., Kierzek, R., Turner, D.H. (1990) Thermodynamic and spectroscopic study of bulge loops in oligoribonucleotides Biochemistry, 29, 278–285 .
Znosko, B.M., Silvestri, S.B., Volkman, H., Boswell, B., Serra, M.J. (2002) Thermodynamic parameters for an expanded nearest-neighbor model for the formation of RNA duplexes with single nucleotide bulges Biochemistry, 41, 10406–10417 .
Shu, Z. and Bevilacqua, P.C. (1999) Isolation and characterization of thermodynamically stable and unstable RNA hairpins from a triloop combinatorial library Biochemistry, 38, 15369–15379 .
Proctor, D.J., Schaak, J.E., Bevilacqua, J.M., Falzone, C.J., Bevilacqua, P.C. (2002) Isolation and characterization of stable tetraloops with the motif YNMG that participates in tertiary interactions Biochemistry, 41, 12062–12075 .
Laing, L.G. and Hall, K.B. (1996) A model of the iron responsive element RNA hairpin loop structure determined from NMR and thermodynamic data Biochemistry, 35, 13586–13596 .
Groebe, D.R. and Uhlenbeck, O.C. (1989) Thermal stability of RNA hairpins containing a four-membered loop and a bulge nucleotide Biochemistry, 28, 742–747 .
Fink, T.R. and Crothers, D.M. (1972) Free energy of imperfect nucleic acid helices,I. The bulge defect J. Mol. Biol, . 66, 1–12 .
Jaeger, J.A., Turner, D.H., Zuker, M. (1989) Improved predictions of secondary structures for RNA Proc. Natl Acad. Sci. USA, 86, 7706–7710 .
Weeks, K.M. and Crothers, D.M. (1993) Major groove accessibility of RNA Science, 261, 1574–1577 .
Chen, G., Znosko, B.M., Jiao, X., Turner, D.H. (2004) Factors affecting thermodynamic stabilities of RNA 3 x 3 internal loops Biochemistry, 43, 12865–12876 .
Gutell, R.R. (1994) Collection of small subunit (16S- and 16S-like) ribosomal RNA structures Nucleic Acids Res, . 22, 3502–3507 .
Gutell, R.R., Gray, M.W., Schnare, M.N. (1993) A compilation of large subunit (23S- and 23S-like) ribosomal RNA structures Nucleic Acids Res, . 21, 3055–3074 .
Schnare, M.N., Damberger, S.H., Gray, M.W., Gutell, R.R. (1996) Comprehensive comparison of structural characteristics in eukaryotic cytoplasmic large subunit (23S-like) ribosomal RNA J. Mol. Biol, . 256, 701–719 .
Szymanski, M., Specht, T., Barciszewska, M.Z., Barciszewski, J., Erdmann, V.A. (1998) 5S rRNA data bank Nucleic Acids Res, . 26, 156–159 .
Sprinzl, M., Horn, C., Brown, M., Ioudovitch, A., Steinberg, S. (1998) Compilation of tRNA sequences and sequences of tRNA genes Nucleic Acids Res, . 26, 148–153 .
Larsen, N., Samuelsson, T., Zwieb, C. (1998) The signal recognition particle database (SRPDB) Nucleic Acids Res, . 26, 177–178 .
Brown, J.W. (1998) The ribonuclease P database Nucleic Acids Res, . 26, 351–352 .
Damberger, S.H. and Gutell, R.R. (1994) A comparative database of group I intron structures Nucleic Acids Res, . 22, 3508–3510 .
Mathews, D.H. and Turner, D.H. (2002) Experimentally derived nearest neighbor parameters for the stability of RNA three- and four-way multibranch loops Biochemistry, 41, 869–880 .
Diamond, J.M., Turner, D.H., Mathews, D.H. (2001) Thermodynamics of three-way multibranch loops in RNA Biochemistry, 40, 6971–6981 .
Walter, A.E., Wu, M., Turner, D.H. (1994) The stability and structure of tandem GA mismatches in RNA depend on closing base pairs Biochemistry, 33, 11349–11354 .
Walter, A.E., Turner, D.H., Kim, J., Lyttle, M.H., Miller, P., Mathews, D.H., Zuker, M. (1994) Coaxial stacking of helixes enhances binding of oligoribonucleotides and improves predictions of RNA folding Proc. Natl Acad. Sci. USA, 91, 9218–9222 .
Kim, J., Walter, A.E., Turner, D.H. (1996) Thermodynamics of coaxially stacked helices with GA and CC mismatches Biochemistry, 35, 13753–13761 .
SantaLucia, J., Jr, Kierzek, R., Turner, D.H. (1991) Functional group substitutions as probes of hydrogen bonding between GA mismatches in RNA internal loops J. Am. Chem. Soc, . 113, 4313–4322 .
Burkard, M.E. and Turner, D.H. (2000) NMR structures of r(GCAGGCGUGC)2 and determinants of stability for single guanosine–guanosine base pairs Biochemistry, 39, 11748–11762 .
Jacobson, H. and Stockmayer, W.H. (1950) Intramolecular reaction in polycondensations. I. The theory of linear systems J. Chem. Phys, . 18, 1600–1606 .
Michel, F., Umesono, K., Ozeki, H. (1989) Comparative and functional anatomy of group II catalytic introns—a review Gene, 82, 5–30 .
Waring, R.B. and Davies, R.W. (1984) Assessment of a model for intron RNA secondary structure relevant to RNA self-splicing—a review Gene, 28, 277–291 .
Kowalak, J.A., Dalluge, J.J., McCloskey, J.A., Stetter, K.O. (1994) The role of posttranscriptional modification in stabilization of transfer RNA from hyperthermophiles Biochemistry, 33, 7869–7876 .
Dubins, D.N., Lee, A., Macgregor, R.B., Jr, Chalikian, T.V. (2001) On the stability of double stranded nucleic acids J. Am. Chem. Soc, . 123, 9254–9259 .
Petersheim, M. and Turner, D.H. (1983) Base-stacking and base-pairing contributions to helix stability: thermodynamics of double-helix formation with CCGG, CCGGp, CCGGAp, ACCGGp, CCGGUp, and ACCGGUp Biochemistry, 22, 256–268 .
Holbrook, J.A., Capp, M.W., Saecker, R.M., Record, M.T., Jr. (1999) Enthalpy and heat capacity changes for formation of an oligomeric DNA duplex: interpretation in terms of coupled processes of formation and association of single-stranded helices Biochemistry, 38, 8409–8422 .
Suurkuusk, J., Alvarez, J., Freire, E., Biltonen, R. (1977) Calorimetric determination of the heat capacity changes associated with the conformational transitions of polyriboadenylic acid and polyribouridylic acid Biopolymers, 16, 2641–2652 .
P?rschke, D., Uhlenbeck, O.C., Martin, F.H. (1973) Thermodynamics and kinetics of the helix-coil transition of oligomers containing GC base pairs Biopolymers, 12, 1313–1335 .
Appleby, D.W. and Kallenbach, N.R. (1973) Theory of oligonucleotide stabilization. I. The effect of single-strand stacking Biopolymers, 12, 2093–2120 .
Freier, S.M., Hill, K.D., Dewey, T.G., Marky, L.A., Breslauer, K.J., Turner, D.H. (1981) Solvent effects on the kinetics and thermodynamics of stacking in poly(cytidylic acid) Biochemistry, 20, 1419–1426 .
Takach, J.C., Mikulecky, P.J., Feig, A.L. (2004) Salt-dependent heat capacity changes for RNA duplex formation J. Am. Chem. Soc, . 126, 6530–6531 .
Mikulecky, P.J., Takach, J.C., Feig, A.L. (2004) Entropy-driven folding of an RNA helical junction: an isothermal titration calorimetric analysis of the hammerhead ribozyme Biochemistry, 43, 5870–5881 .
Mathews, D.H. (2005) Predicting a set of minimal free energy RNA secondary structures common to two sequences Bioinformatics, 21, 2246–2253 .
Lyngs?, R., Zuker, M., Pederson, C. (1999) Fast evaluation of internal loops in RNA secondary structure prediction Bioinformatics, 15, 440–445 .
Layton, D.M. and Bundschuh, R. (2005) A statistical analysis of RNA folding algorithms through thermodynamic parameter perturbation Nucleic Acids Res, . 33, 519–524 .
Borer, P.N., Dengler, B., Tinoco, I., Jr, Uhlenbeck, O.C. (1974) Stability of ribonucleic acid double-stranded helices J. Mol. Biol, . 86, 843–853 .
Khachane, A.N., Timmis, K.N., dos Santos, V.A. (2005) Uracil content of 16S rRNA of thermophilic and psychrophilic prokaryotes correlates inversely with their optimal growth temperatures Nucleic Acids Res, . 33, 4016–4022 .
Adilakshmi, T., Ramaswamy, P., Woodson, S.A. (2005) Protein-independent folding pathway of the 16S rRNA 5' domain J. Mol. Biol, . 351, 508–519 .
Wu, M., McDowell, J.A., Turner, D.H. (1995) A periodic table of symmetric tandem mismatches in RNA Biochemistry, 34, 3204–3211 .
SantaLucia, J., Jr, Kierzek, R., Turner, D.H. (1991) Stabilities of consecutive A·C, C·C, G·G, U·C, and U·U mismatches in RNA internal loops: evidence for stable hydrogen-bonded U·U and C·C+ pairs Biochemistry, 30, 8242–8251 .
Znosko, B.M., Burkard, M.E., Krugh, T.R., Turner, D.H. (2002) Molecular recognition in purine-rich internal loops: thermodynamic, structural, and dynamic consequences of purine for adenine substitutions in 5' (rGGCAAGCCU)2 Biochemistry, 41, 14978–14987 .
Schroeder, S.J. and Turner, D.H. (2001) Thermodynamic stabilities of internal loops with GU closing pairs in RNA Biochemistry, 40, 11509–11517 .
Xia, T., McDowell, J.A., Turner, D.H. (1997) Thermodynamics of nonsymmetric tandem mismatches adjacent to G·C base pairs in RNA Biochemistry, 36, 12486–12487 .(Zhi John Lu, Douglas H. Turner1 and Davi)
*To whom correspondence should be addressed. Tel: 1 585 275 1734; Fax: 1 585 275 6007; Email: david_mathews@urmc.rochester.edu
ABSTRACT
A complete set of nearest neighbor parameters to predict the enthalpy change of RNA secondary structure formation was derived. These parameters can be used with available free energy nearest neighbor parameters to extend the secondary structure prediction of RNA sequences to temperatures other than 37°C. The parameters were tested by predicting the secondary structures of sequences with known secondary structure that are from organisms with known optimal growth temperatures. Compared with the previous set of enthalpy nearest neighbor parameters, the sensitivity of base pair prediction improved from 65.2 to 68.9% at optimal growth temperatures ranging from 10 to 60°C. Base pair probabilities were predicted with a partition function and the positive predictive value of structure prediction is 90.4% when considering the base pairs in the lowest free energy structure with pairing probability of 0.99 or above. Moreover, a strong correlation is found between the predicted melting temperatures of RNA sequences and the optimal growth temperatures of the host organism. This indicates that organisms that live at higher temperatures have evolved RNA sequences with higher melting temperatures.
INTRODUCTION
RNA is more than a simple single-stranded sequence carrying genetic information as in the Central Dogma of Biology. For example, it can form tertiary structures that, such as proteins, can be catalytic. Natural and engineered RNA molecules are widely used as functional tools in enzymatic catalysis and genetic control (1–5). One current problem is how to predict the structures of functional RNA sequences.
Secondary structure, the sum of canonical base pairs, is stronger (6–9) and forms faster (10) than tertiary structure. Therefore, secondary structure can largely be determined without knowledge of tertiary structure. Comparative sequence analysis is a standard technique for determining the secondary structure of homologous RNA sequences (11–13). When only a few or even a single sequence is available, the secondary structure at 37°C can be predicted by free energy minimization algorithms (14–17) using a set of empirical free energy parameters, determined from optical melting experiments (17–21). Each parameter only depends on the sequence identity of nucleotides in the motif and in adjacent base pairs and the total free energy is the sum of nearest neighbor terms. The average sensitivity (the percentage of known base pairs that are correctly predicted) of free energy minimization prediction has been benchmarked as high as 72.8 ± 9.4% for a diverse database of sequences having fewer than 800 nt (17). Furthermore, experimentally determined constraints can improve this accuracy of prediction up to 84% (17,18) for sequences with <6% pseudoknotted (non-nested) base pairs (17). Partition function prediction of base pair probabilities can be used to identify base pairs in the predicted lowest free energy structure that are much more likely than average to be in the known secondary structure (22,23). For example, 91.0% of base pairs in the lowest free energy structure with pairing probability of 0.99 or higher are contained in the known structure, on average (22). The high accuracy of thermodynamic structure prediction (17) demonstrates that many RNA secondary structures can be determined from sequences, without knowledge of any tertiary contacts or protein interactions.
The current set of free energy nearest neighbor parameters for predicting the free energy of RNA secondary structure, however, is limited to application at 37°C. Many organisms, thermophiles and psychrophiles, live at temperatures far from 37°C and many experiments are conducted at other temperatures. The prediction of secondary structure of RNA at arbitrary temperature would expand our knowledge of structure and evolution in the RNA world. Moreover, it would facilitate studying and designing functional RNA molecules at temperatures other than 37°C. The enthalpy nearest neighbor parameters can be used in conjunction with available free energy nearest neighbor parameters for 37°C to determine free energy nearest neighbors at other temperatures. But the most recent enthalpy parameters were derived in 1995 using a simple model (24). At that time, no themes had emerged for the sequence-dependent stability of internal loops. Subsequently, the nearest neighbor model for free energy change at 37°C was significantly improved (17) using experimental results. Therefore, we applied the principles of the current free energy nearest neighbor model (17,18) to determine a complete set of enthalpy nearest neighbor parameters using the available optical melting data.
MATERIALS AND METHODS
Database of experiments
The database of experimental data for derivation of enthalpy parameters is included in Supplementary Data. It includes 130 hairpin loops (25–31), 37 bulge loops (32,33), 337 internal loops (17,18,34–49) (99 of which are 2 x 2 internal loops), 74 multibranch loops (50,51) and 43 coaxial stacking models (52–55).
Derivation and refinement of enthalpy parameters
Canonical base pairs
The enthalpies of Watson–Crick and GU base pairs were derived by Xia et al. (21) and Mathews et al. (18), respectively.
Dangling ends and terminal mismatches
Dangling ends are unpaired nucleotides adjacent to canonical pairs and their enthalpy parameters were compiled previously (24). Dangling ends on terminal GU pairs are treated similar to dangling ends on terminal AU pairs. Terminal mismatches are non-canonical pairs at the end of helixes. The enthalpy parameters of terminal mismatches are taken from another compilation (20), with the exception of mismatches on terminal GU pairs, which were measured recently (30).
If a terminal mismatch has the potential to pair canonically, the values of A–C and C–A mismatches are used for the purine–pyrimidine mismatch and pyrimidine–purine mismatches, respectively. This is important for partition function calculations, where all possible secondary structures are considered.
Hairpin loops
The experimental enthalpies of hairpin loop formation are calculated from published experimental data (25–31) with the following equation:
where is the experimental value for unfolding the hairpin loop with stem, is calculated by the INN-HB parameters (18,21), without an intermolecular initiation term.
The hairpin loop enthalpy parameters are estimated by linear regression using the same model as free energy nearest neighbor parameters (17), except that the GG first mismatch bonus observed for free energy does not apply for enthalpy because the bonus was not statistically significant for enthalpy. The GG stability bonus is therefore entropic in nature, consistent with the observation that GG mismatches are dynamic (56), i.e. they sample more than one single microstate on short timescales.
The enthalpies of hairpin loops are estimated by the following equation:
where n is the number of unpaired nucleotides in the loop. Hairpins with fewer than 3 unpaired nucleotides are not allowed by the model. When n = 3, only the initiation term is considered without any bonus and penalty terms, except a penalty for hairpin loops with three Cs. When n > 3, the special GU closure bonus applies to GU closed hairpins in which a 5' closing G is preceded by two G residues; and (UU or GA first mismatch but not AG) is applied to loops with first mismatches of UU or GA (G on the 5' side and A on 3' side of loop). The oligo-C penalty applies only to loops composed of all C residues and, if n > 3, is calculated with (oligo-C loops, n > 3) = An + B. For hairpinloops composed entirely of 3 C residues, the (oligo-C loops, n = 3) is applied.
The enthalpy parameters are listed in Table 1 and the database of measured loop enthalpies is available as Supplementary Data. In the absence of data, for hairpin loops longer than 9 nt, the initiation enthalpy is approximated with the initiation term for a hairpin of 9 nt. This assumes that additional instability of hairpin loops as the loop lengthens derives from the entropy (57).
Table 1 Hairpin loop enthalpy parametersa
The measured free energies at 37°C of some special hairpin loops of 3, 4 or 6 unpaired nucleotides (30,31,34–36) are either more or less stable by 0.9 kcal/mol than the model predicts. The enthalpies for each of these sequences are listed in a separate lookup table (Table 2), to be consistent with the free energy parameters.
Table 2 Lookup table for unstable triloops and stable tetraloops and hexaloops
Bulge loops
RNA secondary structure is destabilized by bulge loops, which are an interruption of helical structure in one strand only (32,37,38). The initiation terms, for bulge loops of 1–3 nt, are listed in Table 3. They are the average values of experimental data (32,33), calculated using the following equation:
where the enthalpy of the duplex without bulge is the experimental value of the sequence of the duplex without the bulge or as calculated with INN-HB parameters (21) if the experimental values were not available. is the stacking enthalpy of the base pairs in the duplex without the bulge that flank the bulge loop in the duplex with the bulge. Because the difference of initiation enthalpies between 2 and 3 nt bulges is almost zero, it is assumed that the increasing instability for longer bulges (n 4) comes from the entropy of the loop closure (39,57). Thus, the initiation enthalpy for bulges longer than 3 nt is approximated as the 3 nt bulge enthalpy.
Table 3 Bulge loop initiation enthalpy parametersa
Assuming that helical stacking is continuous between the adjacent helices for single bulges, but is interrupted by longer bulges (39,40), the enthalpies of bulge loops are calculated with the following equation:
The calculation of enthalpies for the adjacent helices would include the terminal AU/GU penalty (21) for AU/GU pairs adjacent to the bulge loops that are longer than 1 nt. is the canonical helix stacking enthalpy applied for the two closing base pairs as though the helix was not interrupted by the bulge loop.
Internal loops
Internal loop enthalpies were calculated from experimental data (17,18,34–49) using the following equation:
The range of measured enthalpies differs for internal loops of different size and symmetry; therefore, different enthalpy models are used to predict different loop types. The models are similar to those used to model free energies (17).
1 x 1 Internal loops (single mismatches)
For single non-canonical pairs (1 x 1 internal loops), the loop enthalpies are approximated by the following equation:
where is the enthalpy of initiation for a single non-canonical pair; is the penalty for each AU or GU closing base pair; is a bonus for a GG pair in a 1 x 1 loop; and is a bonus for a 5'RU/3'YU stack in a 1 x 1 loop, where R is a purine and Y is a pyrimidine.
2 x 2 Internal loops (tandem mismatches)
The 2 x 2 internal loops, also called tandem mismatches, interrupt helical RNA with two opposing unpaired nucleotides on each strand. Many of the sequence-symmetric 2 x 2 loops have been studied experimentally (17,18, 34–49) and their enthalpies are assembled in a ‘periodic table’ (Table 4). Symmetric sequences that have not been measured are approximated by averaging the most adjacent columns that have been measured. For asymmetric 2 x 2 loops, the enthalpies are approximated using the following equation:
where (12.5 ± 2.7 kcal/mol) is applied to loops with a GG pair adjacent to an AA or any non-canonical pair with a pyrimidine and p (2.4 ± 3.1 kcal/mol) is applied to loops with an AG or GA pair adjacent to a UC, CU or CC pair or with a UU pair adjacent to an AA pair.
Table 4 The periodic table of tandem mismatch (2 x 2 internal loop) enthalpya
Other internal loops
The enthalpies of other internal loops are approximated using the following equation:
where is the enthalpy of initiation for a loop of n nucleotides; is a penalty for loops with unequal numbers of nucleotides on each side, with n1 and n2 the number of nucleotides on each side; is applied for each sequence-specific first mismatch (Table 5), but it is not applied to loops of the form 1 x (n – 1) with n > 3 (n is the total number of unpaired bases). Special first mismatch bonuses were determined for 2 x 3 and 1 x 2 internal loops with separate linear regressions.
Table 5 Approximations for internal loop enthalpy parameters at 37°C (in kcal/mol)
Moreover, the free energy parameters (Table 6) were updated for internal loops based on recent experimental measurements. The free energy parameters were obtained using the method of Mathews et al. (17). The recent data include the 3 x 3 loops from Chen et al. (41), but excluding the 3 x 3 loops with a middle GA pair. The middle GA pair is shown to enhance stability and this extra stability cannot be predicted by the nearest neighbor parameter set used in this work (41).
Table 6 Updated internal loop free energy parameters at 37°C (in kcal/mol)
Coaxial stacking
Coaxial stacking, which is a favorable interaction of two helices stacked end to end, occurs in multibranch loops and exterior loops. Stability increments for coaxial stacking were measured with a structure composed of a short oligonucleotide bound to a single-stranded end of a stem–loop structure, creating a helical interface (52–55). The enthalpy of coaxial stacking is quantified as follows:
where H°(correction) is the enthalpy for displacing a 3' dangling end on the stem–loop structure if one is present.
When the helixes have no intervening mismatches, the enthalpy bonus is approximated by the nearest neighbor parameter (21) of a base pair in a helix. The excess enthalpy above the helical stacking nearest neighbor from Xia et al. (21), , for each measured interface was calculated. With flush interfaces, i.e. with no intervening mismatch, and no strand extensions beyond the interface, the average excess enthalpy is –1.53 ± 1.45 kcal/mol. For interfaces followed by strand extensions, the excess enthalpy is 1.82 ± 1.13 kcal/mol. As the excess enthalpy changes are not statistically significant, coaxial stacking of helices with no intervening nucleotides is modeled with the enthalpy parameter in a helix.
With one intervening nucleotide from each strand, two helices can stack with an intervening mismatch between them. There are two stack increments: one is the mismatch stack at the end of one helix with continuous backbone, which is equal to the mismatch stacking parameter on a helix, and the other is the mismatch stack with discontinuous backbone, which is modeled as sequence independent. The average enthalpy of sequence independent stacks is –8.46 ± 2.75 kcal/mol. In addition to this, an enthalpy bonus of –0.4 or –0.2 kcal/mol are applied to intervening mismatches composed of nucleotides that could form a Watson–Crick or a GU base pair, respectively. These bonuses are identical to free energy increments that are used and are empirically found to improve structure prediction accuracy.
Multibranch loops
The parameters are determined by linear regression of experimental data for three- and four-way multibranch loops (50,51). In a nearest neighbor model, the bimolecular enthalpy () for the formation of the duplex with a multibranch loop is given by the following equation:
where helix 1 and helix 2 are the intermolecular paired helices with H° predicted from nearest neighbor parameters for Watson–Crick pairs (without including bimolecular initiation so that appears only once). The is a term that accounts for the stacking enthalpy increment of the nucleotides that can stack on the hairpin loop stems to form a modified motif after the two strands have dissociated. This is the most favorable configuration with coaxial stacking of helixes (in the case of four-way multibranch loops) or of the stacking of unpaired nucleotides. is the experimental value which is taken from versus ln(CT/4) plots. The multibranch loop enthalpy initiation term () can be calculated from the above equation. The enthalpy of multibranch loops () is then modeled as the sum of two terms, initiation and stacking:
The stacking term is the favorable enthalpy of coaxial stacking, terminal mismatch and/or dangling end stacking. It is determined from the stacking conformation that gives the lowest free energy, as determined by free energy nearest neighbors (50). The initiation term can be approximated by the following equation:
where a, b and c are parameters determined from linear regression (Table 7) and h is the number of branching helices. is a strain enthalpy that only applies to three-way multibranch loops with fewer than two unpaired nucleotides. The asym term is the average asymmetry that reflects the distribution of unpaired nucleotides, which is defined by the following equation:
Table 7 Enthalpy parameters for multibranch loop initiation
The average asymmetry is limited to 2.0, following the rules suggested by free energy parameters. Asymmetry cannot be applied, however, by dynamic programming algorithms for secondary structure prediction (17,22). Thus, the b term was excluded for secondary structure prediction and the parameters a and c were optimized by finding the parameters that lead to the highest average sensitivity of secondary structure prediction by free energy minimization. The maximum sensitivity of prediction was found with a = 30.0 kcal/mol and c = –2.2 kcal/mol.
Database of RNA secondary structures
The revised enthalpy nearest neighbor model was tested with RNA sequences with known secondary structure from organisms with known optimal growth temperature. The structures were taken from comparative analysis databases (42–49,58,59). Small (16S) subunit rRNA sequences are divided into domains as defined by Jaeger et al. (39). Large (23S) subunit rRNA sequences are divided into domains of fewer than 700 nt each (18). The optimal growth temperatures of different organisms were taken from the Prokaryotic Growth Temperature Database (http://pgtdb.csie.ncu.edu.tw/) and the DSMZ German Collection of Microorganisms and Cell Cultures website (http://www.dsmz.de/). Only the RNA sequences of mesophiles (organisms living at temperatures between 10 and 60°C, but with organisms living at 37°C excluded) were chosen to test the sensitivity and positive predictive value (PPV) of secondary structure prediction. Considering that posttranscriptional modification (60) and high pressure (61) in the thermophiles and hyperthermophiles (organism living above 60°C) would change the thermodynamics of secondary structure formation, sequences from these organisms were excluded. A list of sequences and optimal growth temperatures used are available in Supplementary Data.
Accuracy of secondary structure prediction
The accuracy of structure prediction is determined by the sum of the canonical base pairs correctly predicted. A base pair is considered correctly predicted even if it is shifted by 1 nt on one side. For example a base pair between nucleotides i and j is considered to be correctly predicted if any of these base pairs is predicted: i to j, i to j – 1, i to j + 1, i – 1 to j or i + 1 to j. The predicted base pair between i – 1 and j + 1, however, is not considered to be correct. This scoring scheme reflects the uncertainty of exact base pair matches in comparative sequence analysis and the possibility for dynamics in base pairing. The values of sensitivity and PPV of this scoring scheme are 2–3% higher than when determined with exact base pairing only, where only the i to j base pair is considered to be correct. The prediction accuracies are shown in Supplementary Tables 11 and 12. Each table includes accuracies determined when pairs can be shifted and when pairs must be an exact match.
Availability of parameters
Machine-readable tables of the enthalpy parameters are available on the Mathews lab website (http://rna.urmc.rochester.edu/).
RESULTS
Nearest neighbor model parameters
In the nearest neighbor model of free energy (17,18), the parameters for Watson–Crick base pairs are well determined at 37°C with errors <10%, or 0.1–0.2 kcal/mol (21). For other motifs such as loops and GU base pairs, individual nearest neighbor free energy increments are often determined with an error <0.5 kcal/mol (17,18). In order to extend the current model to predict free energy at temperatures other than 37°C, enthalpy parameters consistent with the current nearest neighbor model are required. The free energy at arbitrary temperature for each parameter is then
(1)
where the enthalpy (H°) and entropy (S°) are assumed to be temperature independent. As described in Materials and Methods, parameters for enthalpy prediction, compatible with the free energy model, were determined using available experimental data from optical melting experiments.
Experimental studies consistently demonstrate that enthalpy and entropy measurements have considerably larger percent error than free energy measurements. Free energy at 37°C is determined with greater precision because of correlation between errors in enthalpy and entropy (21). The larger experimental errors in enthalpy result in larger percent errors for enthalpy nearest neighbor parameters than free energy parameters. The enthalpy of RNA secondary structure is known to be a function of temperature. A linear model for heat capacity change predicts the following:
(2)
(3)
where is a constant heat capacity change and T0 is a chosen reference temperature. It is hypothesized that the heat capacity change arises from the extent of stacking increasing with decreasing temperature. Thus, is negative because single strands are more organized at low rather than high temperature (62–67). The can be estimated by linear fits of enthalpy and entropy changes as a function of melting temperature (50,51,62) or determined by isothermal titration calorimetry at multiple temperatures (68,69). However, the effects of heat capacity change on enthalpy and entropy are antagonistic in terms of free energy change:
(4)
Therefore, for certain T (T = T – T0), can be neglected because the effects are compensated in terms of free energy. To calculate the compensation for a set of RNA duplexes (62), the free energy, G°, was derived directly from Equation 4 assuming that the entropy and enthalpy were independent of temperature. Then the temperature-dependent free energy, , was calculated with the measured non-zero from Equations 2–4. The free energy difference, G° = – G°, increases with the deviation of temperature from T0 (37°C) (Figure 1). The exact G° for each duplex is shown in Table 8 for different temperatures. The experimental error in individual loop free energy nearest neighbor parameters at 37°C is as large as 0.5 kcal/mol (17), which corresponds to roughly a factor of 2 in equilibrium constant. Thus, the small G° for helices suggests that the approximation of = 0 is reasonable for predictions from 10 to 60°C. Therefore, the enthalpy parameters derived here assume = 0 and are most accurate at predicting free energy change close to 37°C.
Figure 1 (A) Free energy difference of RNA duplex CCGGUp. G° (dashed line) was derived from Equation 3, where enthalpy and entropy were averaged from the optical melting curve fits, assuming that they were independent of the temperature. (solid line) was calculated from Equations 1–3, where the heat capacity was accounted. (B) Free energy difference is G° = (62).
Table 8 Free energy differences of RNA duplexes
Dynamic programming algorithm for RNA secondary structure prediction
RNAstructure is a program for RNA secondary structure prediction and analysis. It includes prediction of secondary structure by free energy minimization (17), prediction of base pair probabilities using a partition function (22), the efn2 function for predicting the free energy change of folding given a sequence and secondary structure (18), and the Dynalign algorithm for finding the secondary structure common to two sequences (70). RNAstructure was revised to make predictions at user-defined temperature. Because large internal loops are more likely at high temperature, the previous limitation on internal loop size (fewer than 30 unpaired nucleotides) (17,18,22) was removed by implementing the method of Lyngs? et al. (71). This provides an O(N3) algorithm that can predict internal loops of arbitrary size. Benchmarks for calculation time and memory requirement with and without this revision are shown in Table 9.
Table 9 Calculation time and memory size of dynamic programming for sequences of different length
Sensitivities and PPVs of structure predictions
The enthalpy nearest neighbor parameters were compared with the previous parameters and model for enthalpy and free energy assembled by Serra and Turner (24) by predicting the secondary structures of RNA sequences with known secondary structures. Sensitivities, the percent of known base pairs that are correctly predicted, using both sets of parameters are shown in Figure 2 (detailed numbers are in Supplementary Table 11A) for different types of structural RNA sequences. The known structures of these sequences were taken from comparative analysis databases (42–49,58,59). The average sensitivity is improved from 65.2 to 68.9% using the new parameters assembled here. Sensitivities are improved for most types of the RNA. The exceptions are 5S rRNA and Group II introns.
Figure 2 Improvement of prediction at optimal growth temperatures. The sequences are those from mesophiles (optimal growth temperature from 10 to 60°C) without organisms with optimal growth at 37°C. The lowest free energy secondary structures were predicted at the organims' optimal growth temperatures using two models. The previous model and parameters are those of Serra and Turner (24), which are widely used. The improved prediction uses the model and parameters presented in this work. The small and large subunits of rRNA sequences are divided into domains of <700 nt. The total sensitivity is the average of sensitivities of different types of RNA.
To test the enthalpy parameters, the accuracy of secondary structure prediction at optimal growth temperature was compared to the accuracy of structure prediction at 37°C for organisms that do not grow optimally at 37°C for several types of RNAs (Table 10). The comparison of predictions was shown in different groups divided by optimal growth temperature. The organisms in each group grow optimally in a certain range of temperatures. Compared to the prediction at 37°C, structure prediction at optimal growth temperature performs better for the organism living at temperatures between 22 and 37°C, but is worse at other optimal growth temperatures. This suggests that when enthalpy parameters are assumed to be temperature independent, their utility as a tool for deriving free energy parameters for use in predicting the lowest free energy structure is limited to a narrow temperature range. Small errors in enthalpy change parameters have a larger effect on free energy change parameter determination (Equation 1), the farther the temperature is from 37°C.
Table 10 Prediction sensitivities of the lowest free energy structurea
Figure 3 shows the PPV for base pairs from the lowest free energy structure for base pairs with different pairing probabilities (see detailed numbers in Supplementary Table 12A). They are predicted using a partition function calculation at optimal growth temperature (22). PPV is the percentage of predicted base pairs that are found in the known structure. The average PPV of all pairs in the lowest free energy structures is only 62.0%, which is lower than the sensitivity (68.9%). This suggests that the model over-predicts base pairs and/or that the base pairs may not be annotated completely in the structures from comparative analysis (22). For example, if a base pair is completely conserved, then it is sometimes not annotated by comparative analysis (42–49,58,59). Base pair probabilities for all possible pairs are calculated with a partition function and grouped by different thresholds. The PPV is significantly higher for predicted base pairs in the lowest free energy structure with higher pairing probability. The average PPV is up to 90.4% for those known base pairs having probability of 0.99 or above. It has been demonstrated previously that base pair probabilities predicted at 37°C can be used to find pairs with high PPV (22). The fact that this holds true at other temperatures shows that the enthalpy parameters are robust for base pair probability prediction.
Figure 3 PPV for optimal structure and base pairs with different pairing probabilities. PPV equals the number predicted base pairs in that are in the known structure divided by total number of predicted base pairs. Pairs in the optimal structures are grouped by different thresholds of pairing probabilities. The pairing probabilities were calculated with a partition function calculation (22) at organisms' optimal growth temperatures, using the model and parameters presented in Materials and Methods. The small and large subunits of rRNA sequences are divided into domains of <700 nt. The sequences of different type of RNA are those from mesophiles (living from 10 to 60°C) without organisms living at 37°C.
The fact that the accuracy of secondary structure prediction is sensitive to the accuracy of the nearest neighbor parameters, but the base pair probabilities remain a robust measure of confidence for a wide variety of temperatures is consistent with a previous work. Layton and Bundschuh (72) demonstrated that the predicted lowest free energy structure was often changed in repeated structure predictions after random adjustments of the nearest neighbor parameters within the limits of their error. Base pair probabilities, however, were less perturbed by changes in the parameters (72). With the extrapolation of nearest neighbor parameters to temperatures far from 37°C, the accuracy of the predicted lowest free energy structure is often reduced as compared to structure prediction at 37°C. The ability of the partition function predicted base pair probabilities to determine base pairs predicted with a higher confidence is unchanged with secondary structure prediction at temperatures far from 37°C. This is because the determination of base pair probabilities is not as perturbed by errors in the nearest neighbor parameters.
An example of secondary structure prediction at 37°C and at optimal growth temperature of 30°C is shown in Figure 4 for a tRNA sequence. The base pairs with higher predicted pairing probability (color annotated according to pairing probability in Figure 4B and C) are pairs predicted with greater confidence. For this sequence, secondary structure prediction is more accurate and the fidelity of structure prediction (as judged by the percent of high probability pairs) is improved at optimal growth temperature.
Figure 4 Secondary structure prediction of Saccharomyces cerevisiae tRNA (RM4000) at optimal growth temperature (30°C) (B) and at 37°C (C) with the presented nearest neighbor parameters. Base pairs in the original structure (A) are derived from the comparative analysis database (42–49,58,59). Structures are also color annotated to indicate predicted base pair probabilities (Pbp) for each helix: red, Pbp 0.95; yellow, 0.95 > Pbp 0.7; green, 0.7 > Pbp 0.3; blue, 0.3 > Pbp. The structures were drawn with XRNA (http://rna.ucsc.edu/rnacenter/xrna/xrna.html) and Adobe Illustrator.
Correlation between melting temperature and optimal growth temperature
Melting temperature, Tm, is defined as the temperature at which half of strands are unpaired. Assuming that an RNA melts with a two-state transition, the melting temperature (in Kelvins) of a single-stranded RNA structure can be predicted by Tm = H°/S° (73). For example, the predicted melting temperatures (°C) for all hairpins in the database of optically melted sequences (Supplementary Data) (25–31) are plotted in Figure 5 as a function of experimentally determined Tm. This shows that the parameters adequately reflect the thermal stabilities of RNA sequences with known Tm. Better correlation was found at higher temperatures. This is expected because most hairpins were measured with high melting temperatures in experiments (25–31).
Figure 5 Experimental (Supplementary Data) (25–31) versus predicted (Tm = H°/S° – 273.15) melting temperatures of hairpin stem–loop structures. The line shows the ideal location of points, predicted Tm = measured Tm. The root mean squared deviation (r.m.s.d.) of prediction compared to experiment is 5.86°C. The new enthalpy parameters provide improved Tm prediction compared to the previous compilation of parameters (24), which have an r.m.s.d. of 7.58°C as compared to experiment for this dataset.
Melting temperature reflects the thermal stability of a structure. Therefore RNA structures in organisms living at higher temperature are expected to have higher melting temperatures. Figure 6A shows a plot of predicted melting temperatures of the lowest free energy structure versus organism optimal growth temperature (10–90°C). A strong correlation (linear correlation coefficient of 0.797) is found between the melting temperature and the optimal growth temperature for different types of RNA structures. On the other hand, there appears to be less correlation between nucleotide content and optimal growth temperature (Figure 6B–D) for diverse types of RNA, although uracil content of 16S rRNA of thermophiles and psychrophiles were found recently to correlate inversely with their optimal growth temperatures (74). Evidently, the thermal stability of RNA structure is not simply controlled by base content. Organisms that grow at high temperature have apparently evolved RNA secondary structures with a combination of motifs that provide thermal stability.
Figure 6 Relationships of melting temperatures, nucleotide contents and optimal growth temperatures of different types of RNA in different organisms with optimal growth temperature from 10 to 90°C: (A) Predicted melting temperature; (B) G–C pair content; (C) G content; and (D) U content versus optimal growth temperature. Melting temperatures are predicted for different types of RNA sequences from comparative analysis databases (42–49,58,59) with a two-state transition assumption.
DISCUSSION
The nearest neighbor parameters for enthalpy were derived here using similar rules as for free energy nearest neighbor parameters at 37°C (17). This makes these parameters useful for determining free energy parameters at arbitrary temperature that are compatible with dynamic programming algorithms for secondary structure prediction. Some of the enthalpy parameters have large percent standard errors as compared with the parameters of free energy. This reflects the larger errors in the experimental results of enthalpy than free energy, but it also suggests that enthalpy may be more sequence dependent than free energy. This sequence dependence cannot be determined using the currently available database of optical melting experiments and suggests a need for further optical melting experiments on model RNA systems.
Another source of error comes from the assumption that the enthalpy and entropy are independent of the temperature in both the model and in the analysis of optical melting experiments. When the temperature is too far from 37°C, the sensitivity of prediction is expected to be worse than 68.9% on average because of the approximation of . For example, experiments demonstrate cold denaturation of RNA (68,69), but the nearest neighbor model does not reproduce those results. Further experiments by isothermal titration calorimetry would be needed to provide the data for a model that can include a non-zero heat capacity change.
There are common error sources that should be considered for the prediction of base pairs. Free energy minimization assumes that the secondary structure is at equilibrium. The nearest neighbor model is an incomplete representation of structural free energy. The parameters average some sequence-specific effects and were derived from a limited set of experiments. Some RNA sequences, in particular mRNA, may sample multiple structures at equilibrium. The parameters are derived from experimental data at 1 M NaCl, whereas the salt concentration in different organisms may be very different.
In spite of all these limitations, the nearest neighbor model predicts secondary structures with a 72.8% average sensitivity (17). Recent experimental results on the self-folding of the 16S rRNA 5' domain (75) support the assumption of thermodynamic control of folding pathway. Moreover, the base pair prediction with the partition function can be used to determine pairs predicted with greater confidence (22).
In spite of the fact that the enthalpy parameters have larger percent errors than the free energy parameters for 37°C, the enthalpy parameters are able to predict optical melting temperatures for small model sequences. Predicted melting temperatures for structural RNA sequences correlate well with optimal growth temperature, suggesting that these parameters capture many of the sequence-dependent features of RNA folding enthalpy change.
SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.
ACKNOWLEDGEMENTS
The authors thank Rahul Tyagi and Andrew V. Uzilov for helpful discussions. D.H.M. is an Alfred P. Sloan Research Fellow. This work was supported by National Institutes of Health Grants GM22939 to D.H.T. and GM076485 to D.H.M. Funding to pay the Open Access publication charges for this article was provided by National Institutes of Health.
REFERENCES
Nelson, P., Kiriakidou, M., Sharma, A., Maniataki, E., Mourelatos, Z. (2003) The microRNA world: small is mighty Trends Biochem. Sci, . 28, 534–540 .
Doudna, J. and Cech, T. (2002) The chemical repertoire of natural ribozymes Nature, 418, 222–228 .
Walter, P. and Blobel, G. (1982) Signal recognition particle contains a 7S RNA essential for protein translocation across the endoplasmic reticulum Nature, 299, 691–698 .
Lau, N.C., Lim, L.P., Weinstein, E.G., Bartel, D.P. (2001) An abundant class of tiny RNAs with probable regulatory roles in Caenorhabditis elegans Science, 294, 858–862 .
Lagos-Quintana, M., Rauhut, R., Lendeckel, W., Tuschl, T. (2001) Identification of novel genes coding for small expressed RNAs Science, 294, 853–858 .
Onoa, B. and Tinoco, I., Jr. (2004) RNA folding and unfolding Curr. Opin. Struct. Biol, . 14, 374–379 .
Crothers, D.M., Cole, P.E., Hilbers, C.W., Schulman, R.G. (1974) The molecular mechanism of thermal unfolding of Escherichia coli formylmethionine transfer RNA J. Mol. Biol, . 87, 63–88 .
Mathews, D.H., Banerjee, A.R., Luan, D.D., Eickbush, T.H., Turner, D.H. (1997) Secondary structure model of the RNA recognized by the reverse transcriptase from the R2 retrotransposable element RNA, 3, 1–16 .
Banerjee, A.R., Jaeger, J.A., Turner, D.H. (1993) Thermal unfolding of a group I ribozyme: the low temperature transition is primarily a disruption of tertiary structure Biochemistry, 32, 153–163 .
Woodson, S.A. (2000) Recent insights on RNA folding mechanisms from catalytic RNA Cell Mol. Life Sci, . 57, 796–808 .
James, B.D., Olsen, G.J., Pace, N.R. (1989) Phylogenetic comparative analysis of RNA secondary structure Methods Enzymol, . 180, 227–239 .
Pace, N.R., Thomas, B.C., Woese, C.R. (1999) Probing RNA structure, function, and history by comparative analysis In Gesteland, R.F., Cech, T.R., Atkins, J.F. (Eds.). The RNA World, 2nd edn Cold Spring Harbor Cold Spring Harbor Laboratory Press pp. 113–141 .
Woese, C.R., Gutell, R.R., Gupta, R., Noller, H.F. (1983) Detailed analysis of the higher order structure of 16S-like ribosomal ribonucleic acids Microbiol. Rev, . 47, 621–669 .
Andronescu, M., Aguirre-Hernandez, R., Condon, A., Hoos, H.H. (2003) RNAsoft: A suite of RNA secondary structure prediction and design software tools Nucleic Acids Res, . 31, 3416–3422 .
Hofacker, I.L. (2003) Vienna RNA secondary structure server Nucleic Acids Res, . 31, 3429–3431 .
Zuker, M. (2003) Mfold web server for nucleic acid folding and hybridization prediction Nucleic Acids Res, . 31, 3406–3415 .
Mathews, D.H., Disney, M.D., Childs, J.L., Schroeder, S.J., Zuker, M., Turner, D.H. (2004) Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure Proc. Natl Acad. Sci. USA, 101, 7287–7292 .
Mathews, D.H., Sabina, J., Zuker, M., Turner, D.H. (1999) Expanded sequence dependence of thermodynamic parameters provides improved prediction of RNA Secondary Structure J. Mol. Biol, . 288, 911–940 .
Turner, D.H. (2000) Conformational changes In Bloomfield, V., Crothers, D., Tinoco, I. (Eds.). Nucleic Acids, Sausalito, CA University Science Books pp. 259–334 .
Xia, T., Mathews, D.H., Turner, D.H. (1999) Thermodynamics of RNA secondary structure formation In Soll, D.G., Nishimura, S., Moore, P.B. (Eds.). Prebiotic Chemistry, Molecular Fossils, Nucleosides, and RNA, NY Elsevier pp. 21–47 .
Xia, T., SantaLucia, J., Jr, Burkard, M.E., Kierzek, R., Schroeder, S.J., Jiao, X., Cox, C., Turner, D.H. (1998) Thermodynamic parameters for an expanded nearest-neighbor model for formation of RNA duplexes with Watson–Crick pairs Biochemistry, 37, 14719–14735 .
Mathews, D.H. (2004) Using an RNA secondary structure partition function to determine confidence in base pairs predicted by free energy minimization RNA, 10, 1178–1190 .
McCaskill, J.S. (1990) The equilibrium partition function and base pair probabilities for RNA secondary structure Biopolymers, 29, 1105–1119 .
Serra, M.J. and Turner, D.H. (1995) Predicting thermodynamic properties of RNA Methods Enzymol, . 259, 242–261 .
Giese, M.R., Betschart, K., Dale, T., Riley, C.K., Rowan, C., Sprouse, K.J., Serra, M.J. (1998) Stability of RNA hairpins closed by wobble base pairs Biochemistry, 37, 1094–1100 .
Serra, M.J., Lyttle, M.H., Axenson, T.J., Schadt, C.A., Turner, D.H. (1993) RNA hairpin loop stability depends on closing pair Nucleic Acids Res, . 21, 3845–3849 .
Serra, M.J., Axenson, T.J., Turner, D.H. (1994) A model for the stabilities of RNA hairpins based on a study of the sequence dependence of stability for hairpins of six nucleotides Biochemistry, 33, 14289–14296 .
Serra, M.J., Barnes, T.W., Betschart, K., Gutierrez, M.J., Sprouse, K.J., Riley, C.K., Stewart, L., Temel, R.E. (1997) Improved parameters for the prediction of RNA hairpin stability Biochemistry, 36, 4844–4851 .
Groebe, D.R. and Uhlenbeck, O.C. (1988) Characterization of RNA hairpin loop stability Nucleic Acids Res, . 16, 11725–11735 .
Dale, T., Smith, R., Serra, M.J. (2000) A test of the model to predict unusually stable RNA hairpin loop stability RNA, 6, 608–615 .
Antao, V.P. and Tinoco, I., Jr. (1992) Thermodynamic parameters for loop formation in RNA and DNA hairpin tetraloops Nucleic Acids Res, . 20, 819–824 .
Longfellow, C.E., Kierzek, R., Turner, D.H. (1990) Thermodynamic and spectroscopic study of bulge loops in oligoribonucleotides Biochemistry, 29, 278–285 .
Znosko, B.M., Silvestri, S.B., Volkman, H., Boswell, B., Serra, M.J. (2002) Thermodynamic parameters for an expanded nearest-neighbor model for the formation of RNA duplexes with single nucleotide bulges Biochemistry, 41, 10406–10417 .
Shu, Z. and Bevilacqua, P.C. (1999) Isolation and characterization of thermodynamically stable and unstable RNA hairpins from a triloop combinatorial library Biochemistry, 38, 15369–15379 .
Proctor, D.J., Schaak, J.E., Bevilacqua, J.M., Falzone, C.J., Bevilacqua, P.C. (2002) Isolation and characterization of stable tetraloops with the motif YNMG that participates in tertiary interactions Biochemistry, 41, 12062–12075 .
Laing, L.G. and Hall, K.B. (1996) A model of the iron responsive element RNA hairpin loop structure determined from NMR and thermodynamic data Biochemistry, 35, 13586–13596 .
Groebe, D.R. and Uhlenbeck, O.C. (1989) Thermal stability of RNA hairpins containing a four-membered loop and a bulge nucleotide Biochemistry, 28, 742–747 .
Fink, T.R. and Crothers, D.M. (1972) Free energy of imperfect nucleic acid helices,I. The bulge defect J. Mol. Biol, . 66, 1–12 .
Jaeger, J.A., Turner, D.H., Zuker, M. (1989) Improved predictions of secondary structures for RNA Proc. Natl Acad. Sci. USA, 86, 7706–7710 .
Weeks, K.M. and Crothers, D.M. (1993) Major groove accessibility of RNA Science, 261, 1574–1577 .
Chen, G., Znosko, B.M., Jiao, X., Turner, D.H. (2004) Factors affecting thermodynamic stabilities of RNA 3 x 3 internal loops Biochemistry, 43, 12865–12876 .
Gutell, R.R. (1994) Collection of small subunit (16S- and 16S-like) ribosomal RNA structures Nucleic Acids Res, . 22, 3502–3507 .
Gutell, R.R., Gray, M.W., Schnare, M.N. (1993) A compilation of large subunit (23S- and 23S-like) ribosomal RNA structures Nucleic Acids Res, . 21, 3055–3074 .
Schnare, M.N., Damberger, S.H., Gray, M.W., Gutell, R.R. (1996) Comprehensive comparison of structural characteristics in eukaryotic cytoplasmic large subunit (23S-like) ribosomal RNA J. Mol. Biol, . 256, 701–719 .
Szymanski, M., Specht, T., Barciszewska, M.Z., Barciszewski, J., Erdmann, V.A. (1998) 5S rRNA data bank Nucleic Acids Res, . 26, 156–159 .
Sprinzl, M., Horn, C., Brown, M., Ioudovitch, A., Steinberg, S. (1998) Compilation of tRNA sequences and sequences of tRNA genes Nucleic Acids Res, . 26, 148–153 .
Larsen, N., Samuelsson, T., Zwieb, C. (1998) The signal recognition particle database (SRPDB) Nucleic Acids Res, . 26, 177–178 .
Brown, J.W. (1998) The ribonuclease P database Nucleic Acids Res, . 26, 351–352 .
Damberger, S.H. and Gutell, R.R. (1994) A comparative database of group I intron structures Nucleic Acids Res, . 22, 3508–3510 .
Mathews, D.H. and Turner, D.H. (2002) Experimentally derived nearest neighbor parameters for the stability of RNA three- and four-way multibranch loops Biochemistry, 41, 869–880 .
Diamond, J.M., Turner, D.H., Mathews, D.H. (2001) Thermodynamics of three-way multibranch loops in RNA Biochemistry, 40, 6971–6981 .
Walter, A.E., Wu, M., Turner, D.H. (1994) The stability and structure of tandem GA mismatches in RNA depend on closing base pairs Biochemistry, 33, 11349–11354 .
Walter, A.E., Turner, D.H., Kim, J., Lyttle, M.H., Miller, P., Mathews, D.H., Zuker, M. (1994) Coaxial stacking of helixes enhances binding of oligoribonucleotides and improves predictions of RNA folding Proc. Natl Acad. Sci. USA, 91, 9218–9222 .
Kim, J., Walter, A.E., Turner, D.H. (1996) Thermodynamics of coaxially stacked helices with GA and CC mismatches Biochemistry, 35, 13753–13761 .
SantaLucia, J., Jr, Kierzek, R., Turner, D.H. (1991) Functional group substitutions as probes of hydrogen bonding between GA mismatches in RNA internal loops J. Am. Chem. Soc, . 113, 4313–4322 .
Burkard, M.E. and Turner, D.H. (2000) NMR structures of r(GCAGGCGUGC)2 and determinants of stability for single guanosine–guanosine base pairs Biochemistry, 39, 11748–11762 .
Jacobson, H. and Stockmayer, W.H. (1950) Intramolecular reaction in polycondensations. I. The theory of linear systems J. Chem. Phys, . 18, 1600–1606 .
Michel, F., Umesono, K., Ozeki, H. (1989) Comparative and functional anatomy of group II catalytic introns—a review Gene, 82, 5–30 .
Waring, R.B. and Davies, R.W. (1984) Assessment of a model for intron RNA secondary structure relevant to RNA self-splicing—a review Gene, 28, 277–291 .
Kowalak, J.A., Dalluge, J.J., McCloskey, J.A., Stetter, K.O. (1994) The role of posttranscriptional modification in stabilization of transfer RNA from hyperthermophiles Biochemistry, 33, 7869–7876 .
Dubins, D.N., Lee, A., Macgregor, R.B., Jr, Chalikian, T.V. (2001) On the stability of double stranded nucleic acids J. Am. Chem. Soc, . 123, 9254–9259 .
Petersheim, M. and Turner, D.H. (1983) Base-stacking and base-pairing contributions to helix stability: thermodynamics of double-helix formation with CCGG, CCGGp, CCGGAp, ACCGGp, CCGGUp, and ACCGGUp Biochemistry, 22, 256–268 .
Holbrook, J.A., Capp, M.W., Saecker, R.M., Record, M.T., Jr. (1999) Enthalpy and heat capacity changes for formation of an oligomeric DNA duplex: interpretation in terms of coupled processes of formation and association of single-stranded helices Biochemistry, 38, 8409–8422 .
Suurkuusk, J., Alvarez, J., Freire, E., Biltonen, R. (1977) Calorimetric determination of the heat capacity changes associated with the conformational transitions of polyriboadenylic acid and polyribouridylic acid Biopolymers, 16, 2641–2652 .
P?rschke, D., Uhlenbeck, O.C., Martin, F.H. (1973) Thermodynamics and kinetics of the helix-coil transition of oligomers containing GC base pairs Biopolymers, 12, 1313–1335 .
Appleby, D.W. and Kallenbach, N.R. (1973) Theory of oligonucleotide stabilization. I. The effect of single-strand stacking Biopolymers, 12, 2093–2120 .
Freier, S.M., Hill, K.D., Dewey, T.G., Marky, L.A., Breslauer, K.J., Turner, D.H. (1981) Solvent effects on the kinetics and thermodynamics of stacking in poly(cytidylic acid) Biochemistry, 20, 1419–1426 .
Takach, J.C., Mikulecky, P.J., Feig, A.L. (2004) Salt-dependent heat capacity changes for RNA duplex formation J. Am. Chem. Soc, . 126, 6530–6531 .
Mikulecky, P.J., Takach, J.C., Feig, A.L. (2004) Entropy-driven folding of an RNA helical junction: an isothermal titration calorimetric analysis of the hammerhead ribozyme Biochemistry, 43, 5870–5881 .
Mathews, D.H. (2005) Predicting a set of minimal free energy RNA secondary structures common to two sequences Bioinformatics, 21, 2246–2253 .
Lyngs?, R., Zuker, M., Pederson, C. (1999) Fast evaluation of internal loops in RNA secondary structure prediction Bioinformatics, 15, 440–445 .
Layton, D.M. and Bundschuh, R. (2005) A statistical analysis of RNA folding algorithms through thermodynamic parameter perturbation Nucleic Acids Res, . 33, 519–524 .
Borer, P.N., Dengler, B., Tinoco, I., Jr, Uhlenbeck, O.C. (1974) Stability of ribonucleic acid double-stranded helices J. Mol. Biol, . 86, 843–853 .
Khachane, A.N., Timmis, K.N., dos Santos, V.A. (2005) Uracil content of 16S rRNA of thermophilic and psychrophilic prokaryotes correlates inversely with their optimal growth temperatures Nucleic Acids Res, . 33, 4016–4022 .
Adilakshmi, T., Ramaswamy, P., Woodson, S.A. (2005) Protein-independent folding pathway of the 16S rRNA 5' domain J. Mol. Biol, . 351, 508–519 .
Wu, M., McDowell, J.A., Turner, D.H. (1995) A periodic table of symmetric tandem mismatches in RNA Biochemistry, 34, 3204–3211 .
SantaLucia, J., Jr, Kierzek, R., Turner, D.H. (1991) Stabilities of consecutive A·C, C·C, G·G, U·C, and U·U mismatches in RNA internal loops: evidence for stable hydrogen-bonded U·U and C·C+ pairs Biochemistry, 30, 8242–8251 .
Znosko, B.M., Burkard, M.E., Krugh, T.R., Turner, D.H. (2002) Molecular recognition in purine-rich internal loops: thermodynamic, structural, and dynamic consequences of purine for adenine substitutions in 5' (rGGCAAGCCU)2 Biochemistry, 41, 14978–14987 .
Schroeder, S.J. and Turner, D.H. (2001) Thermodynamic stabilities of internal loops with GU closing pairs in RNA Biochemistry, 40, 11509–11517 .
Xia, T., McDowell, J.A., Turner, D.H. (1997) Thermodynamics of nonsymmetric tandem mismatches adjacent to G·C base pairs in RNA Biochemistry, 36, 12486–12487 .(Zhi John Lu, Douglas H. Turner1 and Davi)