Kengo SatopersName, Manato AkiyamapersName, Yasubumi SakakibarapersName (2021date). “RNA secondary structure prediction using deep learning with thermodynamic integration..” Nature Communicationstitle·j.
Object of study — deep-learning prediction of RNA secondary structureterm·key with thermodynamic integrationterm·mth (MXfold2term·mth)
Research question — whether integrating Turner free-energy parameters into DNN folding scores reduces overfittingterm·mth
Contribution — MXfold2term·mth: DNN folding scores integrated with nearest-neighborterm·key free energies and trained with thermodynamic regularizationterm·mth
phenomenon · claim-type · cross-domain-transfer · mechanism · claim-relation · replication-status · limitation — grounded in the full text| phenomenon | overfitting of richly parameterized RNA structure predictors to training families — “However, rich parameterization can easily cause overfitting to the training data, thus preventing robust predictions for a wide range of RNA sequences 15 .” | |
| claim-type | method — “Here we show that overfitting can be minimized when RNA folding scores learnt using a deep neural network are integrated together with Turner’s nearest-neighbor free energy parameters.” | |
| method — evidence | — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” | |
| object of study — evidence | — “Here we show that overfitting can be minimized when RNA folding scores learnt using a deep neural network are integrated together with Turner’s nearest-neighbor free energy parameters.” | |
| research question — evidence | — “Although machine learning-based models have achieved high performance in terms of prediction accuracy, overfitting is a common risk for such highly parameterized models.” | |
| contribution — evidence | — “Here we show that overfitting can be minimized when RNA folding scores learnt using a deep neural network are integrated together with Turner’s nearest-neighbor free energy parameters.” | |
| mechanism | thermodynamic regularization Keeping folding scores close to Turner free energies stops the DNN overfitting to training families. — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” | mechanisms[0] |
| sense | family-wise cross-validation - testing on RNA families absent from training, the paper's test of robustness — “We observe that E2Efold almost completely failed to predict secondary structures for unseen families ( F = 0.0361), whereas MXfold2 accurately predicted their secondary structures ( F = 0.628).” | senses[0] |
| relation | RNA secondary structures —operationalised-as→ base-pairs — defined as — “Therefore, instead of conducting such experimental assays, we frequently perform computational prediction of RNA secondary structures, defined as sets of base-pairs with hydrogen bonds between the nucleotides.” — “RNA secondary structures, defined as sets of base-pairs” | concept_edges[0] |
| weight | RNA secondary structures —operationalised-as→ base-pairs — 0.85 — “Therefore, instead of conducting such experimental assays, we frequently perform computational prediction of RNA secondary structures, defined as sets of base-pairs with hydrogen bonds between the nucleotides.” | concept_edges[0] |
| relation | DP —same-as→ dynamic programming — ( ) — “We can efficiently calculate an optimal secondary structure that has the minimum free energy using a dynamic programming (DP) technique, the well-known Zuker algorithm 5 .” — “using a dynamic programming (DP) technique” | concept_edges[1] |
| weight | DP —same-as→ dynamic programming — 0.95 — “We can efficiently calculate an optimal secondary structure that has the minimum free energy using a dynamic programming (DP) technique, the well-known Zuker algorithm 5 .” | concept_edges[1] |
| relation | SCFGs —same-as→ stochastic context-free grammars — ( ) — “Probabilistic generative models such as stochastic context-free grammars (SCFGs) have also been applied to predicting RNA secondary structures.” — “stochastic context-free grammars (SCFGs)” | concept_edges[2] |
| weight | SCFGs —same-as→ stochastic context-free grammars — 0.95 — “Probabilistic generative models such as stochastic context-free grammars (SCFGs) have also been applied to predicting RNA secondary structures.” | concept_edges[2] |
| relation | max-margin framework —same-as→ structured support vector machine — also known as — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” — “the max-margin framework, also known as structured support vector machine” | concept_edges[3] |
| weight | max-margin framework —same-as→ structured support vector machine — 0.9 — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” | concept_edges[3] |
| effect direction | rich parameterization —causes→ overfitting — increases — “However, rich parameterization can easily cause overfitting to the training data, thus preventing robust predictions for a wide range of RNA sequences 15 .” — “rich parameterization can easily cause overfitting” | concept_edges[4] |
| relation | rich parameterization —causes→ overfitting — cause — “However, rich parameterization can easily cause overfitting to the training data, thus preventing robust predictions for a wide range of RNA sequences 15 .” — “rich parameterization can easily cause overfitting to the training data” | concept_edges[4] |
| weight | rich parameterization —causes→ overfitting — 0.8 — “However, rich parameterization can easily cause overfitting to the training data, thus preventing robust predictions for a wide range of RNA sequences 15 .” | concept_edges[4] |
| effect direction | thermodynamic integration —causes→ robustness — increases — “Our experiments showed that thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined, substantially improved the robustness of predictions.” — “substantially improved the robustness of predictions” | concept_edges[5] |
| relation | thermodynamic integration —causes→ robustness — improved — “Our experiments showed that thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined, substantially improved the robustness of predictions.” — “substantially improved the robustness of predictions” | concept_edges[5] |
| weight | thermodynamic integration —causes→ robustness — 0.85 — “Our experiments showed that thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined, substantially improved the robustness of predictions.” | concept_edges[5] |
| relation | thermodynamic regularization —part-of→ thermodynamic integration — consisting of — “Our experiments showed that thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined, substantially improved the robustness of predictions.” — “thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization” | concept_edges[6] |
| weight | thermodynamic regularization —part-of→ thermodynamic integration — 0.9 — “Our experiments showed that thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined, substantially improved the robustness of predictions.” | concept_edges[6] |
| relation | MXfold2 —contrasted-with→ SPOT-RNA — whereas — “Whereas MXfold2 is limited to pseudoknot-free secondary structure prediction owing to its use of the nearest-neighbor model, SPOT-RNA and E2Efold are able to predict pseudoknotted structures because their RNA secondary structure prediction is formulated as multiple binary classification problems that predict whether each pair of nucleotides forms a base pair or not.” — “Whereas MXfold2 is limited to pseudoknot-free secondary structure prediction owing to its use of the nearest-neighbor model, SPOT-RNA and E2Efold are able to predict pseudoknotted structures” | concept_edges[7] |
| weight | MXfold2 —contrasted-with→ SPOT-RNA — 0.85 — “Whereas MXfold2 is limited to pseudoknot-free secondary structure prediction owing to its use of the nearest-neighbor model, SPOT-RNA and E2Efold are able to predict pseudoknotted structures because their RNA secondary structure prediction is formulated as multiple binary classification problems that predict whether each pair of nucleotides forms a base pair or not.” | concept_edges[7] |
| replication | original | |
| limitation | MXfold2 cannot predict pseudoknotted structures because it uses the nearest-neighbor model. — “Whereas MXfold2 is limited to pseudoknot-free secondary structure prediction owing to its use of the nearest-neighbor model, SPOT-RNA and E2Efold are able to predict pseudoknotted structures because their RNA secondary structure prediction is formulated as multiple binary classification problems that predict whether each pair of nucleotides forms a base pair or not.” | limitations[0] |
| limitation | Where neither machine-learning nor thermodynamic methods succeed, MXfold2 also tends to fail. — “However, Fig. 5 suggests that if neither machine learning-based methods nor the thermodynamic-based method can accurately predict the secondary structure, MXfold2 may also be difficult to predict it successfully.” | limitations[1] |
| limitation | Training data are hard to enlarge because empirically determined structures are costly. — “Although the development of sequencing techniques has made it easier to obtain sequence data, it is still not easy to increase the amount of training data because of the significant effort required to obtain those empirically determined secondary structures.” | limitations[2] |
finding · quantity · variable · definition · key-term — grounded in the full text| finding | A DNN-only base model overfits: high accuracy on TestSetA, poor on TestSetB. — “The results in Table 1 and Supplementary Fig. 1 show that the base model achieved high accuracy for TestSetA, but quite poor accuracy for TestSetB, suggesting that the base model might be categorized as being prone to heavy overfitting.” | content_tags.findings[0] |
| finding | E2Efold almost completely fails on unseen families; MXfold2 does not. — “We observe that E2Efold almost completely failed to predict secondary structures for unseen families ( F = 0.0361), whereas MXfold2 accurately predicted their secondary structures ( F = 0.628).” | content_tags.findings[1] |
| finding | Accurate structure prediction does not by itself give reliable thermodynamic stability estimates. — “This suggests that high accuracy in the secondary structure prediction does not directly lead to the reliable evaluation of thermodynamic stability.” | content_tags.findings[2] |
| finding | MXfold2 is about 15 (GPU) and 36 (CPU) times faster than SPOT-RNA. — “In addition, as shown in Supplementary Fig. 3 , MXfold2 is, on average, approximately 15 and 36 times faster than SPOT-RNA when using GPU and CPU, respectively.” | content_tags.findings[3] |
| quantity | F = 0.761 — — “For the predictions on TestSetA, our algorithm, MXfold2, achieved the best accuracy ( F = 0.761), followed by ContextFold ( F = 0.759).” | content_tags.quantities[0] |
| quantity | 803k parameters — — “For example, our algorithm employs 803k parameters, whereas SPOT-RNA and E2Efold use 1746k and 719k parameters, respectively.” | content_tags.quantities[1] |
| variable | thermodynamic regularization — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” | content_tags.variables[0] |
| definition | RNA secondary structures — “Therefore, instead of conducting such experimental assays, we frequently perform computational prediction of RNA secondary structures, defined as sets of base-pairs with hydrogen bonds between the nucleotides.” | content_tags.definitions[0] |
| key-terms | RNA secondary structure, deep neural network, MXfold2, thermodynamic regularization, nearest-neighbor, overfitting, free energy, cross-validation, SPOT-RNA, E2Efold |
| concept | RNA secondary structures — sets of base-pairs with hydrogen bonds between the nucleotides — “Therefore, instead of conducting such experimental assays, we frequently perform computational prediction of RNA secondary structures, defined as sets of base-pairs with hydrogen bonds between the nucleotides.” | concepts[0] |
| concept | base-pairs — not defined in this paper's own words; paired nucleotides — “Therefore, instead of conducting such experimental assays, we frequently perform computational prediction of RNA secondary structures, defined as sets of base-pairs with hydrogen bonds between the nucleotides.” | concepts[1] |
| concept | DP — dynamic programming — “We can efficiently calculate an optimal secondary structure that has the minimum free energy using a dynamic programming (DP) technique, the well-known Zuker algorithm 5 .” | concepts[2] |
| concept | dynamic programming — not defined in this paper's own words; the Zuker-style algorithm that finds the optimal structure — “We can efficiently calculate an optimal secondary structure that has the minimum free energy using a dynamic programming (DP) technique, the well-known Zuker algorithm 5 .” | concepts[3] |
| concept | SCFGs — stochastic context-free grammars — “Probabilistic generative models such as stochastic context-free grammars (SCFGs) have also been applied to predicting RNA secondary structures.” | concepts[4] |
| concept | stochastic context-free grammars — not defined in this paper's own words; probabilistic generative grammars applied to RNA folding — “Probabilistic generative models such as stochastic context-free grammars (SCFGs) have also been applied to predicting RNA secondary structures.” | concepts[5] |
| concept | max-margin framework — structured support vector machine (SSVM) — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” | concepts[6] |
| concept | structured support vector machine — not defined in this paper's own words; the max-margin training framework — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” | concepts[7] |
| concept | rich parameterization — not defined in this paper's own words; models with very many trainable parameters — “However, rich parameterization can easily cause overfitting to the training data, thus preventing robust predictions for a wide range of RNA sequences 15 .” | concepts[8] |
| concept | overfitting — not defined in this paper's own words; fitting the training data at the cost of robust prediction — “However, rich parameterization can easily cause overfitting to the training data, thus preventing robust predictions for a wide range of RNA sequences 15 .” | concepts[9] |
| concept | thermodynamic integration — consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined — “Our experiments showed that thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined, substantially improved the robustness of predictions.” | concepts[10] |
| concept | thermodynamic regularization — prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters — “The deep neural network is trained using the max-margin framework, also known as structured support vector machine (SSVM), to minimize the structured hinge loss function with thermodynamic regularization, which prevents the folding score of the secondary structure from differing significantly from the free energy of the thermodynamic parameters.” | concepts[11] |
| concept | robustness — not defined in this paper's own words; accuracy on sequences structurally unlike the training data — “Our experiments showed that thermodynamic integration, consisting of thermodynamic-integrated folding scores and thermodynamic regularization combined, substantially improved the robustness of predictions.” | concepts[12] |
| concept | MXfold2 — not defined in this paper's own words; the paper's deep-learning folding algorithm — “Whereas MXfold2 is limited to pseudoknot-free secondary structure prediction owing to its use of the nearest-neighbor model, SPOT-RNA and E2Efold are able to predict pseudoknotted structures because their RNA secondary structure prediction is formulated as multiple binary classification problems that predict whether each pair of nucleotides forms a base pair or not.” | concepts[13] |
| concept | SPOT-RNA — not defined in this paper's own words; a DNN-based RNA structure predictor — “Whereas MXfold2 is limited to pseudoknot-free secondary structure prediction owing to its use of the nearest-neighbor model, SPOT-RNA and E2Efold are able to predict pseudoknotted structures because their RNA secondary structure prediction is formulated as multiple binary classification problems that predict whether each pair of nucleotides forms a base pair or not.” | concepts[14] |
Full text — tagged
RNA secondary structureterm·key prediction using deep learning with thermodynamic integrationterm·mth.
Abstract Accurate predictions of RNA secondary structuresterm·dom can help uncover the roles of functional non-coding RNAs. Although machine learning-based models have achieved high performance in terms of prediction accuracy, overfittingterm·mth is a common risk for such highly parameterized models. Here we show that overfittingterm·mth can be minimized when RNA folding scores learnt using a deep neural networkterm·key are integrated together with Turner’s nearest-neighborterm·key free energyterm·key parameters. Training the model with thermodynamic regularizationterm·mth ensures that folding scores and the calculated free energyterm·key are as close as possible. In computational experiments designed for newly discovered non-coding RNAs, our algorithm (MXfold2term·mth) achieves the most robust and accurate predictions of RNA secondary structuresterm·dom without sacrificing computational efficiency compared to several other algorithms. The results suggest that integrating thermodynamic information could help improve the robustnessterm·mth of deep learning-based predictions of RNA secondary structureterm·key.
Subject terms: Structure determination, RNA, Machine learning, Non-coding RNAs Introduction Recent studies have revealed that functional non-coding RNAs (ncRNAs) play essential roles, including transcriptional regulation and guiding modification, participating in various biological processes ranging from development to cell differentiation, with defective functionality being involved in various diseases 1 . Because it is well-known that the functions of ncRNAs are deeply related to their structures rather than their primary sequences, discovering the structures of ncRNAs can elucidate the functions of ncRNAs. However, there are major difficulties in determining RNA tertiary structures through experimental assays such as nuclear magnetic resonance and X-ray crystal structure analysis, because of the high experimental costs and resolution limits on measurements of RNA. Although considerable advances in cryo-electron microscopy research on RNA tertiary structure determination have been achieved in recent years 2 , these limitations have not yet been completely overcome. Therefore, instead of conducting such experimental assays, we frequently perform computational prediction of RNA secondary structuresterm·dom, defined as sets of base-pairsterm·dom with hydrogen bonds between the nucleotides.
The most popular approach for predicting RNA secondary structuresterm·dom is based on thermodynamic models, such as Turner’s nearest-neighborterm·key model 3 , 4 , in which a secondary structure is decomposed into several characteristic substructures, called nearest-neighborterm·key loops, such as hairpin loops, internal loops, bulge loops, base-pairterm·dom stackings, multi-branch loops, and external loops, as shown in Fig. 1 . The free energyterm·key of each nearest-neighborterm·key loop can be calculated by summing the free energyterm·key parameters that characterize the loop. The free energyterm·key parameters have been determined in advance by experimental methods such as optical melting experiments 3 . The free energyterm·key of an entire RNA secondary structureterm·key is calculated by summing the free energyterm·key of the decomposed nearest-neighborterm·key loops. We can efficiently calculate an optimal secondary structure that has the minimum free energyterm·key using a dynamic programmingterm·mth (DPterm·mth) technique, the well-known Zuker algorithm 5 . A number of tools, including Mfold/UNAfold 6 , 7 , RNAfold 8 , 9 , and RNAstructure 10 , 11 , have adopted this approach.
Decomposition of an RNA secondary structureterm·key into nearest-neighborterm·key loops. An RNA secondary structureterm·key can be decomposed into several types of nearest-neighborterm·key loops, including hairpin loops (e.g., bases 11–19), internal loops (e.g., bases 25–29 and 43–47), bulge loops (e.g., bases 4–5 and 57–62), base-pairterm·dom stackings (e.g., bases 23–24 and 48–49), multi-branch loops (e.g., bases 7–9, 21–23, and 49–55), and external loops (e.g., bases 1–2 and 64–66). This diagram was drawn using VARNA 45 .
An alternative approach utilizes machine learning techniques, which train scoring parameters for decomposed substructures from reference structures, rather than experimental techniques. This approach has successfully been adopted by CONTRAfold 12 , 13 , ContextFold 14 , and other methods, and has enabled us to more accurately predict RNA secondary structuresterm·dom. However, rich parameterizationterm·mth can easily cause overfittingterm·mth to the training data, thus preventing robust predictions for a wide range of RNA sequences 15 . Probabilistic generative models such as stochastic context-free grammarsterm·mth (SCFGsterm·mth) have also been applied to predicting RNA secondary structuresterm·dom. Recently, TORNADO 15 implemented an application of SCFGsterm·mth to the nearest-neighborterm·key model, thus achieving performance comparable with that of other machine learning-based methods.
Hybrid methods that combine thermodynamic and machine learning-based approaches to compensate for each other’s shortcomings have been developed. SimFold 16 , 17 more accurately estimates thermodynamic parameters from training data, including RNA sequences and their known secondary structures, as well as the free energyterm·key of their known secondary structures. We previously proposed MXfold 18 , which combines thermodynamic energy parameters and rich-parameterized weight parameters. The model learns more precise parameters for substructures observed in training data and avoids overfittingterm·mth rich-parameterized weight parameters to the training data by resorting to the thermodynamic parameters for assessing previously unobserved substructures.
In recent years, deep learning has made remarkable progress in a wide range of fields, including bioinformatics and life sciences. SPOT-RNAterm·mth 19 and E2Efoldterm·key 20 have been developed for RNA secondary structureterm·key prediction using deep neural networksterm·key (DNNs). Both algorithms formulate RNA secondary structureterm·key prediction as multiple binary classification problems that predict whether each pair of nucleotides forms a base pair or not, using a deep neural networkterm·key trained with a large amount of training data. As with the other machine learning-based methods described above, concerns about overfittingterm·mth caused by rich parameterizationterm·mth remain.
Inspired by MXfold and the DNN-based RNA secondary structureterm·key prediction methods, in this paper, we propose an algorithm for predicting RNA secondary structuresterm·dom using deep learning. Similar to MXfold, we integrate folding scores, which are calculated by a deep neural networkterm·key, with Turner’s nearest-neighborterm·key free energyterm·key parameters. The deep neural networkterm·key is trained using the max-margin frameworkterm·mth with thermodynamic regularizationterm·mth, which is proposed in this paper to make our model robust by maintaining that the folding score calculated by our method and the free energyterm·key calculated by the thermodynamic parameters are as close as possible. To confirm the robustnessterm·mth of our algorithm, we conducted computational experiments with two types of cross-validationterm·key: sequence-wise cross-validationterm·key, with test datasets structurally similar to training datasets, and family-wise cross-validationterm·key, with test datasets structurally dissimilar to training datasets. The resulting method, called MXfold2term·mth, achieved the best prediction accuracy based not only on the sequence-wise cross-validationterm·key, but also on the family-wise cross-validationterm·key, indicating the robustnessterm·mth of our algorithm. In addition, the computational experiment conducted on a dataset with sequence-structure-energy triplets shows a high correlation between folding scores predicted by MXfold2term·mth and the free energyterm·key derived from optical melting experiments.
Results Overview of our algorithm Our algorithm computes four types of folding scores for each pair of nucleotides by using a deep neural networkterm·key, as shown in Fig. 2 . The folding scores are used to calculate scores of the nearest-neighborterm·key loops. Similar to MXfold, we integrate the folding scores, which are calculated using a deep neural networkterm·key, with Turner’s nearest-neighborterm·key free energyterm·key parameters. Then, our algorithm predicts an optimal secondary structure that maximizes the sum of the scores of the nearest-neighborterm·key loops using Zuker-style dynamic programmingterm·mth (DPterm·mth) 5 . The deep neural networkterm·key is trained using the max-margin frameworkterm·mth, also known as structured support vector machineterm·mth (SSVM), to minimize the structured hinge loss function with thermodynamic regularizationterm·mth, which prevents the folding score of the secondary structure from differing significantly from the free energyterm·key of the thermodynamic parameters.
The network structure of our algorithm. It calculates four types of folding scores for each pair of nucleotides ( i , j ) from a given RNA sequence of length L . See the “Methods” section for more details. Act.: the activation function, Norm.: the normalization function, Conv: the convolution layer, d : the dimension of trainable embeddings, N 1 , N 2 , and N 3 : the numbers of times Blocks 1, 2, and 3 repeat, respectively, BiLSTM: bidirectional long short-term memory, MLP: multilayer perceptron, and Concat 2D: concatenation of 1D sequences to 2D matrices.
Effectiveness of two thermodynamic-related techniques To investigate the effectiveness of the thermodynamic-integrated folding scores and thermodynamic regularizationterm·mth, which are the main contribution of this study, we conducted an experiment to compare the cases with and without the thermodynamic-related techniques. We implemented a base model that employs only the deep neural networkterm·key to compute the folding score; that is, it does not use the thermodynamic parameters for the folding score and the thermodynamic regularizationterm·mth. Then, we compared the base model with the use of the thermodynamic-integrated folding score and/or the thermodynamic regularizationterm·mth. For this experiment, we employed TrainSetA, TestSetA, and TestSetB, which have been established by Rivas et al. 15 . Note that TrainSetA and TestSetB are structurally dissimilar, whereas TrainSetA and TestSetA have some structural similarity. All the models were trained with TrainSetA and were then evaluated using TestSetA and TestSetB. The results in Table 1 and Supplementary Fig. 1 show that the base model achieved high accuracy for TestSetA, but quite poor accuracy for TestSetB, suggesting that the base model might be categorized as being prone to heavy overfittingterm·mth. In contrast, the model using the complete thermodynamic-related techniques (the final model) showed much better accuracy for TestSetB compared with the base model, indicating that the two thermodynamic-related techniques enable robust predictions.
Comparison of the accuracy of the secondary structure prediction between the base model and the use of the thermodynamic-related techniques.
Comparison with the existing methods We compared our algorithm with nine available folding algorithms: MXfold version 0.0.2 18 , LinearFold-V (committed on Feb 5, 2020) 21 , CONTRAfold version 2.02 12 , 13 , ContextFold version 1.00 14 , CentroidFold version 0.0.16 22 , 23 , TORNADO version 0.3 15 with ViennaRNA grammar, RNAfold in ViennaRNA package version 2.4.14 8 , 9 , SimFold version 2.1 16 , 17 , and RNAstructure version 6.2 10 , 11 . The parameters of four trainable algorithms, namely, MXfold, CONTRAfold, ContextFold and TORNADO, were trained with TrainSetA. For the other algorithms, default parameters were used.
Figures 3 a and 3 b show PPV-SEN plots of prediction accuracy for TestSetA and TestSetB, respectively. The maximum expected accuracy (MEA)-based secondary structure predictions, TORNADO, CONTRAfold and CentroidFold, can control PPV and SEN through the parameter γ , and therefore, their accuracies are shown as the curves for multiple γ ∈ {2 n ∣ − 5 ≤ n ≤ 10}.
Comparison of prediction accuracy with the existing methods. Positive predictive value (PPV) and sensitivity (SEN) plots for ( a ) TestSetA, ( b ) TestSetB, and ( c ) the combined dataset with TestSetA and TestSetB. ( d ) F -values of the combined test data for all the methods. See Supplementary Tables 2 – 4 and Supplementary Fig. 5 for more details.
For the predictions on TestSetA, our algorithm, MXfold2term·mth, achieved the best accuracy ( F = 0.761), followed by ContextFold ( F = 0.759). While the difference in F -value between MXfold2term·mth and ContextFold was not significant ( p = 0.49, one-sided Wilcoxon singed-rank test), MXfold2term·mth was significantly more accurate than the other methods, including TORNADO ( F = 0.746 at γ = 4.0), MXfold ( F = 0.739), and CONTRAfold ( F = 0.719 at γ = 4.0) ( p < 0.001). On the other hand, for the predictions on TestSetB, which is structurally dissimilar from TrainSetA, we observe that ContextFold achieved the worst accuracy ( F = 0.502), consistent with the possible proneness of ContextFold to overfittingterm·mth, as noted by Rivas et al. 15 . In contrast, we do not observe heavy overfittingterm·mth, as occurred with ContextFold, with MXfold2term·mth ( F = 0.601), which is also significantly more accurate than the other methods, including MXfold ( F = 0.581), CONTRAfold ( F = 0.573 at γ = 4.0), and TORNADO ( F = 0.552 at γ = 4.0) ( p < 0.001).
Figures 3 c and 3 d show the PPV-SEN plot and the distribution of F -values, respectively, on the combined dataset comprising TestSetA and TestSetB. These results indicate that MXfold2term·mth ( F = 0.693) achieved the best accuracy, followed by the trainable methods, namely, MXfold ( F = 0.673), TORNADO ( F = 0.664 at γ = 4.0), CONTRAfold ( F = 0.658 at γ = 4.0), and ContextFold ( F = 0.651), and MXfold2term·mth also outperformed the thermodynamics-based methods ( p < 0.001, one-sided Wilcoxon singed-rank test).
Comparison with DNN-based methods We compared our algorithm, MXfold2term·mth, with two recently published DNN-based methods: SPOT-RNAterm·mth 19 and E2Efoldterm·key 20 . We trained MXfold2term·mth using the available SPOT-RNAterm·mth and E2Efoldterm·key training datasets and tested the prediction accuracy on their test datasets for direct comparison because SPOT-RNAterm·mth does not provide training modules. Furthermore, because their datasets were designed only for sequence-wise cross-validationterm·key, for the purpose of family-wise cross-validationterm·key, we developed the bpRNA-new dataset, which does not contain any families present in their training datasets.
Table 2 and Supplementary Fig. 2 show the comparison of the prediction accuracy among MXfold2term·mth, SPOT-RNAterm·mth, TORNADO, ContextFold, and RNAfold, with MXfold2term·mth and ContextFold trained on the TR0 dataset, which is a subset of the bpRNA-1m dataset 24 . In contrast, SPOT-RNAterm·mth used the pre-trained model trained on the TR0 dataset and an additional dataset for transfer learning. All of these methods were tested using the TS0 dataset, which is also a subset of the bpRNA-1m dataset, for sequence-wise cross-validationterm·key, after which the bpRNA-new dataset was used for family-wise cross-validationterm·key.
Comparison of the accuracy of secondary structure prediction among MXfold2term·mth, SPOT-RNAterm·mth, TORNADO, ContextFold, and RNAfold.
We observe that MXfold2term·mth performs slightly worse in sequence-wise cross-validationterm·key relative to SPOT-RNAterm·mth, but the difference is not significant ( p = 0.31, one-sided Wilcoxon signed-rank test). On the other hand, in terms of family-wise cross-validationterm·key, MXfold2term·mth is significantly more accurate than SPOT-RNAterm·mth ( p < 0.001). In addition, as shown in Supplementary Fig. 3 , MXfold2term·mth is, on average, approximately 15 and 36 times faster than SPOT-RNAterm·mth when using GPU and CPU, respectively.
Figures 4 and 5 show typical examples of successful and failed predictions by MXfold2term·mth. We find in Fig. 4 that our method, which integrates the thermodynamic model, enables accurate predictions even when machine learning-based methods cannot be sufficiently trained. However, Fig. 5 suggests that if neither machine learning-based methods nor the thermodynamic-based method can accurately predict the secondary structure, MXfold2term·mth may also be difficult to predict it successfully.
Comparison of predicted structures of FJ746361 .1_1-144. Here, we employed MXfold2term·mth (PPV = 1.0, Sen = 0.85, F = 0.92), SPOT-RNAterm·mth (PPV = 0.18, Sen = 0.37, F = 0.24), TORNADO (PPV = 0.967, Sen = 0.725, F = 0.829), ContextFold (PPV = 0.32, Sen = 0.23, F = 0.26), and RNAfold (PPV = 0.82, Sen = 0.80, F = 0.81).
Comparison of predicted structures of URS0000D6A925_12908_1-83. Here, we employed MXfold2term·mth (PPV = 0.12, Sen = 0.19, F = 0.155), SPOT-RNAterm·mth (PPV = 0.12, Sen = 0.08, F = 0.09), TORNADO (PPV = 0.0, Sen = 0.0, F = 0.0), ContextFold (PPV = 0.0, Sen = 0.0, F = 0.0), and RNAfold (PPV = 0.0, Sen = 0.0, F = 0.0).
We also compared MXfold2term·mth with E2Efoldterm·key, TORNADO, ContextFold, and RNAfold (Supplementary Table 5 and Supplementary Fig. 4 ), in which E2Efoldterm·key used the pre-trained model trained on a subset of the RNAStrAlign dataset 25 , and MXfold2term·mth and ContextFold were trained on the same dataset. All of these methods were tested on a subset of the ArchiveII dataset 26 for sequence-wise cross-validationterm·key and the bpRNA-new dataset for family-wise cross-validationterm·key. We observe that E2Efoldterm·key almost completely failed to predict secondary structures for unseen families ( F = 0.0361), whereas MXfold2term·mth accurately predicted their secondary structures ( F = 0.628). Furthermore, E2Efoldterm·key trained with TrainSetA also failed to predict secondary structures for TestSetB ( F = 0.0322), which is outside of Fig. 3 . These results suggest that E2Efoldterm·key might be categorized as being prone to heavy overfittingterm·mth.
Correlation with free energyterm·key We investigated the correlation between the free energyterm·key and the predicted folding score using the T-Full dataset 17 , which contains sequence-structure-energy triplets. Table 3 shows the root mean square error (RMSE) and the Spearman’s rank correlation coefficient ( ρ ) between the free energyterm·key of the reference structure and the folding score of the predicted structure for MXfold2term·mth with and without thermodynamic regularizationterm·mth compared with CONTRAfold and RNAfold. The folding score of MXfold2term·mth is highly correlated with the true free energyterm·key, though it is not as high as that of RNAfold. In contrast, even though MXfold2term·mth without thermodynamic regularizationterm·mth is able to predict secondary structures as accurately as MXfold2term·mth with thermodynamic regularizationterm·mth, its folding score is not highly correlated with the true free energyterm·key. This suggests that high accuracy in the secondary structure prediction does not directly lead to the reliable evaluation of thermodynamic stability.
Comparison of the correlation between predicted folding score and the reference free energyterm·key.
Discussion We have proposed a deep learning-based algorithm for predicting RNA secondary structuresterm·dom with thermodynamic integrationterm·mth. Our experiments showed that thermodynamic integrationterm·mth, consisting of thermodynamic-integrated folding scores and thermodynamic regularizationterm·mth combined, substantially improved the robustnessterm·mth of predictions. We conducted benchmarks that compare our algorithm with conventional algorithms, using not only sequence-wise cross-validationterm·key, but also family-wise cross-validationterm·key, which assumes one of the most practical situations for RNA secondary structureterm·key prediction from single sequences. In this comparison, our algorithm achieved the best prediction accuracy without overfittingterm·mth to the training data, as was the case with other rich-parameterized methods such as ContextFold and E2Efoldterm·key. The proposed MXfold2term·mth algorithm should be useful for improving RNA structure modeling, especially for newly discovered RNAs.
Our experiments showed that several methods with rich parameterizationterm·mth achieved excellent performance for test datasets structurally similar to training datasets, but poor performance for test datasets structurally dissimilar from training datasets. Rivas et al. 15 have already pointed out that rich-parameterized methods might easily tend towards overfittingterm·mth. Such secondary structure prediction methods would be useful if we could know a priori whether the sequence we want to predict is a member of a family included in the training data. However, if we have knowledge about the family of the sequences whose structures are to be predicted, homology searches for RNA sequences, such as Infernal 27 , are a better choice. In other words, the lack of prior knowledge about sequences is a major reason the prediction of secondary structures from single sequences is an important goal. In principle, regardless of whether a given sequence is within a family included in a training dataset or not, robust methods should be able to properly predict its secondary structure because its secondary structure is formed under the same physical laws. In general, machine learning-based methods can be expected to improve prediction accuracy by increasing the number of training samples. Although the development of sequencing techniques has made it easier to obtain sequence data, it is still not easy to increase the amount of training data because of the significant effort required to obtain those empirically determined secondary structures. Therefore, it is impractical to predict secondary structures of newly discovered sequences using methods that are not generalizable to a wide variety of sequences, such as those methods that rely on only memorizing sequences and their secondary structures within a large number of neural network parameters. We designed our experiments to predict secondary structures of newly discovered RNAs using test datasets structurally dissimilar from training datasets, e.g., TestSetB for TrainSetA (originally developed by Rivas et al. 15 ) and the bpRNA-new dataset for the TR0 dataset. These conditions are among the most important circumstances for RNA secondary structureterm·key prediction from single sequences, and our method, which is able to achieve high accuracy under these conditions, is useful in practice.
Deep neural network-based methods, including our method, improve prediction accuracy by optimizing a larger number of parameters for training data relative to conventional methods. For example, our algorithm employs 803k parameters, whereas SPOT-RNAterm·mth and E2Efoldterm·key use 1746k and 719k parameters, respectively. Therefore, as shown so far, such rich-parametrized models need to be trained more carefully because they are prone to overfittingterm·mth. SPOT-RNAterm·mth mitigates overfittingterm·mth by building an ensemble of five deep models with transfer learning. To do that, our MXfold2term·mth integrates the folding scores, which are calculated by a deep neural networkterm·key, with Turner’s nearest-neighborterm·key free energyterm·key parameters, and employs training using thermodynamic regularizationterm·mth. Whereas MXfold2term·mth is limited to pseudoknot-free secondary structure prediction owing to its use of the nearest-neighborterm·key model, SPOT-RNAterm·mth and E2Efoldterm·key are able to predict pseudoknotted structures because their RNA secondary structureterm·key prediction is formulated as multiple binary classification problems that predict whether each pair of nucleotides forms a base pair or not.
Estimation of free energies is important for applications other than structure prediction, such as small interfering RNA selection using hybridization thermodynamics 28 . Table 3 indicates that RNAfold and MXfold2term·mth with thermodynamic regularizationterm·mth can calculate folding scores that are highly correlated with true free energyterm·key estimates, at least for sequences for which secondary structures can be predicted with high accuracy. Given that MXfold2term·mth is more accurate in secondary structure prediction compared with RNAfold, as mentioned above, MXfold2term·mth can be expected to evaluate thermodynamic stability with high accuracy for more sequences compared with RNAfold.
Methods Using a deep neural networkterm·key to compute folding scores Figure 2 shows the overview of our deep neural networkterm·key, which is used to calculate four types of folding scores for each pair of nucleotides in a given sequence: helix stacking, helix opening, helix closing, and unpaired region.
The input of the network is an RNA sequence of length L . Each base (A, C, G, or U) is encoded into d -dimensional embeddings using a trainable embedding function.
The L × d -dimensional sequence embedding is entered into 1D convolution blocks (Block 1) repeated N 1 times. Each 1D convolution block consists of two convolutional layers, with kernel sizes of 5 and 3, respectively, and a depth of d . Layer normalization 29 and continuously differentiable exponential linear units (CELU) activation functions 30 are employed. A dropout rate of 50% is used to avoid overfittingterm·mth 31 . To efficiently train deep neural networksterm·key, a residual connection 32 that skips each layer is employed.
The next block (Block 2) contains N 2 layers of the bidirectional long short-term memory (BiLSTM) network 33 . The number of hidden units in each direction of the BiLSTM network is d /2. Each BiLSTM layer is followed by the layer normalization and CELU activation functions. A dropout rate of 50% is also used in BiLSTM layers.
The L × d matrix obtained so far is divided into two matrices of size L × d /2; one of the two matrices expands in the column direction and the other in the row direction, and then the two matrices are concatenated into an L × L × d matrix, as shown in the “Concat 2D” operation in Fig. 2 . This matrix is entered into 2D convolution blocks (Block 3) N 3 times. Each 2D convolution block consists of two convolutional layers, each with a kernel size of 5 × 5 and 3 × 3, respectively, and a depth of d . The layer normalization and CELU functions with a dropout rate of 50% are again utilized.
Finally, four types of folding scores for each pair of i -th and j -th nucleotides are calculated by 3-layer multilayer perceptron (MLP) with h hidden units.
We used the following values for the above hyperparameters: d = 64, N 1 = 4, N 2 = 2, N 3 = 4, and h = 32.
Using DPterm·mth to predicting RNA secondary structureterm·key Similar to traditional folding methods, such as RNAfold and CONTRAfold, our algorithm adopts the nearest-neighborterm·key model 4 , 34 as an RNA secondary structureterm·key model without pseudoknots. As an alternative to the thermodynamic parameters, our algorithm calculates the free energyterm·key of the nearest-neighborterm·key loops using four types of folding scores computed by the deep neural networkterm·key described in the previous section. The four types of the folding scores are helix stacking scores, helix opening scores, helix closing scores, and unpaired region scores (Fig. 2 ). The helix stacking scores characterize the stacking of base-pairsterm·dom that mainly contributes to the thermodynamic stability of RNA structures. The helix opening scores and helix closing scores characterize each end of helix regions, corresponding to the terminal mismatch in Turner’s thermodynamic model. The unpaired region scores characterize unstructured sequences that contribute to the thermodynamic instability of RNA structures. Our algorithm calculates the folding score for a nearest-neighborterm·key loop by summing up the four types of the scores comprising the loop. For example, the folding score of the internal loop consisting of bases 25–29 and 43–47 in Fig. 1 is calculated by summing the helix closing score for (25, 47), the helix opening score for (28, 44), and the unpaired region scores for (26, 28) and (44, 46). Similarly, the folding score of the base-pairterm·dom stacking consisting of 23–24 and 48–49 in Fig. 1 is calculated from the helix stacking score for (23, 49).
Similar to MXfold, we integrate our model with the thermodynamic model by adding the negative value of the Turner free energyterm·key 4 to the folding scores calculated by the deep neural networkterm·key for each nearest-neighborterm·key loop. We define our scoring function of a secondary structure y given an RNA sequence x as where f T ( x , y ) is the contribution of the thermodynamic parameters (i.e., the negative value of the free energyterm·key of y ) and f W ( x , y ) is the contribution of the deep neural networkterm·key, which is calculated as the sum of the folding scores over the decomposed nearest-neighborterm·key loops as described above. The thermodynamic-integrated folding scores defined by Eq. ( 1 ) can also be decomposed into nearest-neighborterm·key loops as each term of Eq. ( 1 ) can be decomposed in the same manner. Therefore, we can efficiently calculate an optimal secondary structure that maximizes the scoring function (Eq. 1 ) using the Zuker-style dynamic programmingterm·mth (DPterm·mth) algorithm 5 . See Supplementary Methods for more details of the Zuker-style DPterm·mth algorithm.
Training parameters with the max-margin frameworkterm·mth To optimize the network parameters λ , we employed a max-margin frameworkterm·mth called the structured support vector machineterm·mth (SSVM) 35 . Given a training dataset D = { ( x ( k ) , y ( k ) ) } k = 1 K , where x ( k ) is the k -th RNA sequence and y ( k ) is the reference secondary structure for the k -th sequence x ( k ) , we aim to find a set of parameters λ that minimizes the objective function where S ( x ) is a set of all possible secondary structures of x . The first term is the structured hinge loss function 35 . Here, Δ ( y , y ^ ) is a margin term of y ^ for y defined as where δ FN and δ FP are tunable hyperparameters that control the trade-off between sensitivity and specificity for learning the parameters. We used δ FN = 0.5 and δ FP = 0.005 by default. The margin term for structured models enables robust predictions by maximizing the margin between f ( x , y ) and f ( x , y ^ ) for y ≠ y ^ . We can calculate the first term of Eq. ( 2 ) using Zuker-style dynamic programmingterm·mth 5 modified by the use of loss-augmented inference 35 . The second and third terms of Eq. ( 2 ) are regularization terms that penalize parameter values that are extreme relative to a postulated distribution. The second term is our proposed thermodynamic regularizationterm·mth, which prevents the folding score of the secondary structure from differing significantly from the free energyterm·key of the thermodynamic parameters. The third term of Eq. ( 2 ) is the ℓ 2 regularization term. We used C 1 = 0.125 and C 2 = 0.01 by default. To minimize the objective function ( 2 ), we employ the Adam optimizer 36 .
Datasets To evaluate our algorithm, we performed computational experiments on several datasets. Supplementary Table 1 shows a summary of the datasets used in our experiments.
The first dataset, which has been carefully established by Rivas et al. 15 , includes TrainSetA, TestSetA, TrainSetB, and TestSetB. TrainSetA and TestSetA were collected from the literature 12 , 16 , 17 , 37 – 39 . TrainSetB and TestSetB, which contain 22 families with 3D structure annotations, were extracted from Rfam 10.0 40 . The sequences in Train/TestSetB share less than 70% sequence identity with the sequences in TrainSetA. We excluded a number of sequences that contain pseudoknotted secondary structures in the original data sources from all four sub-datasets since all algorithms evaluated in this study were designed for predicting RNA secondary structuresterm·dom without pseudoknots. It is important to note that literature-based TrainSetA and Rfam-based TestSetB are structurally dissimilar, whereas TrainSetA and TestSetA have some structural similarity.
To compare our algorithm with SPOT-RNAterm·mth 19 , we employed the same dataset extracted from the bpRNA-1m dataset 24 , which is based on Rfam 12.2 41 with 2588 families. After removing redundant sequences based on sequence identity, the bpRNA-1m dataset was randomly split into three sub-datasets, namely, TR0, VL0, and TS0, for training, validation, and testing, respectively. Note that this partition is not family-wise, but sequence-wise. SPOT-RNAterm·mth performed the initial training and validation using TR0 and VL0, which was followed by testing using TS0 and then transfer learning using other PDB-based datasets.
To confirm the robustnessterm·mth against “unseen” families, we built a dataset that includes families from the most recent Rfam database, Rfam 14.2 42 . Since the release of Rfam 12.2, from which bpRNA-1m is derived, the Rfam project has been actively collecting about 1,500 RNA families, including families detected by newly developed techniques 43 . We first extracted these newly discovered families. Then, as with SPOT-RNAterm·mth, we removed redundant sequences by CD-HIT-EST 44 , with a cutoff threshold of 80%, and discarded sequences whose lengths exceeded 500 bp. Our dataset is referred to as bpRNA-new.
We also compared our algorithm with E2Efoldterm·key 20 , using the same experimental conditions in which the RNAStrAlign dataset 25 and the ArchiveII dataset 26 were employed for training and testing, respectively. Note that this setting is also not family-wise, but sequence-wise, because both datasets contain the same families.
To investigate the correlation between folding score and free energyterm·key, we employed the T-Full dataset, which was compiled by Andronescu et al. 17 . Each sequence in the T-Full dataset has a corresponding reference secondary structure, as well as its free energyterm·key change derived from optical melting experiments.
Performance measure We evaluated the accuracy of predicting RNA secondary structuresterm·dom through the positive predictive value (PPV) and the sensitivity (SEN) of the base pairsterm·dom, defined as where TP is the number of correctly predicted base pairsterm·dom (true positives), FP is the number of incorrectly predicted base pairsterm·dom (false positives), and FN is the number of base pairsterm·dom in the reference structure that were not predicted (false negatives). We also used the F -value as the balanced measure between PPV and SEN, which is defined as their harmonic mean: For the experiment on the T-Full dataset, we evaluated the root mean square error (RMSE) where e i and e ^ i are the free energyterm·key of the reference secondary structure and the negative value of the folding score of the predicted secondary structure for the i -th sequence, respectively. We also evaluated the Spearman’s rank correlation coefficient ( ρ ) between e and e ^ .
Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Supplementary information Supplementary Information Peer Review File Reporting Summary Acknowledgements This work was partially supported by a Grant-in-Aid for Scientific Research (B) (No. 19H04210) and Challenging Research (Exploratory) (No. 19K22897) from the Japan Society for the Promotion of Science (JSPS) to K.S. and a Grant-in-Aid for JSPS Fellows (No. 18J21767) from JSPS to M.A. This work was also supported by a Grant-in-Aid for Scientific Research on Innovative Areas “Frontier Research on Chemical Communications.” The supercomputer system used for this research was provided by the National Institute of Genetics (NIG), Research Organization of Information and Systems (ROIS).
Author contributions K.S. conceived the study, implemented the algorithm, conducted experiments, and drafted the manuscript. K.S. and M.A. discussed the network architectures, designed the experiments, and collected the datasets. K.S. and Y.S. discussed the training algorithm. All authors read, contributed to the discussion of, and approved the final manuscript.
Data availability The datasets used in this study are available at 10.5281/zenodo.4430150.
Code availability The MXfold2term·mth source code is available at https://github.com/keio-bioinformatics/mxfold2term·mth/ , and the MXfold2term·mth web server is available for use at http://www.dna.bio.keio.ac.jp/mxfold2term·mth/ .
Competing interests The authors declare no competing interests.
Footnotes Supplementary information The online version contains supplementary material available at 10.1038/s41467-021-21194-4.
References Associated Data Supplementary Materials Supplementary Information Peer Review File Reporting Summary Data Availability Statement The datasets used in this study are available at 10.5281/zenodo.4430150.
The MXfold2term·mth source code is available at https://github.com/keio-bioinformatics/mxfold2term·mth/ , and the MXfold2term·mth web server is available for use at http://www.dna.bio.keio.ac.jp/mxfold2term·mth/ .