In this research, we designed a deep-learning workflow that integrates knowledge of different thermodynamic structure models with constrained optimization to tackle the RNA secondary structure folding problem. The designed integration algorithm is called DEBFold (deep ensemble boosting and folding). The DEBFold workflow comprises two stages to convert the given RNA sequences into their potential secondary structures (Figure 1): (1) structure location folding probability estimation. DEBFold was designed to take an RNA sequence into a deep-learning network that helps to integrate knowledge from various thermal models. Then, it generates the folding constraints similar to the SHAPE experiment results. (2) SHAPE-like constrained optimization. In the second stage, free energy minimization constrained by the predicted structure location folding probability will help provide the final structure prediction. These two stages are elaborated on in the following subsections.

Model training in data-based learning approaches requires the data diversity in the training set to mimic the degree of variations in real problems.³⁸ Researchers have experimented that a large number of sequences from a small number of structural families can easily lead to model overfitting.²³ The possible reason for the situation is that structurally similar sequences provide limited variations for deep models to fetch the inherent patterns that map sequences to the corresponding structures. In order to prepare a ground-truth data set that incorporates a high structure diversity, we downloaded the data deposited in bpRNA-1m³⁹ and retrieved the sequences with known Rfam families. In bpRNA-1m, RNAs curated from Rfam 12.2⁴⁰ were selected and processed. In total, 43,269 valid RNA sequences from 2125 families were found.

In order to tune and evaluate DEBFold in a structural family independent manner, the ground-truth data set is divided and sampled into two parts: the training-validation and test sets. Since RNAs in the same family are structurally similar,²³ the data processing was carried out in a family-wise manner. We resorted to 25-fold cross-validation in DEBFold Stage 1 model training due to the scarcity of known distinct RNA structural families in the collected ground-truth data set. However, the average of the cross-validation results provides only an optimistic estimation of the model performance. An extra test set should be set aside as totally clean data for a fair performance evaluation. Therefore, among the gathered 2125 families, 44 (2%, or 1/50 due to 25-fold cross-validation, of the total 2125 families) of them were reserved in the test set. These 44 families were selected by stratified sampling that preserves the length distribution of all 2125 RNA families. In the calculation of the length distribution of RNA families, the longest sequence within a structural family was chosen as the representative. These 44 families were also ensured to be in distinct Rfam clans from the rest of the remaining families. The remaining 2081 families were used as the training-validation folds in cross-validation. We sampled at most 5 sequences from the 2081 families to be included in the training-validation set (8782 sequences in total). Finally, we took one sequence from each of the 44 test families for testing, resulting in a test set consisting of 44 RNAs. This test set is termed TestSetα in this research. We further confirmed the sequence similarity between the training-validation and test RNA sequences using CD-HIT.⁴¹ Under the threshold of 0.85, TestSetα contains no RNA with sequence similarity to any member of the training-validation set. In order to have more extensive testing, we further prepared two additional independent test sets (TestSetβ and TestSetγ) as described in the following two subsections.

To evaluate the RNA secondary structure prediction performance, we resorted to the structure F1 score calculation. The precision, recall, and F1 score for a predicted structure are computed as follows^3,81516where TP, FP, and FN stand for the numbers of true-positive, false-positive, and false–negative pairs, respectively. For a predicted structure and its corresponding real structure, TP counts the number of correctly predicted pairings, FP is the number of predicted pairings that do not appear in the real structure, and FN considers the number of real base pairings missed by the prediction. The final F1 score can balance the trade-off between the FP and FN values and is thus used as the overall evaluation metric.

We first evaluated the generalization prediction performance of DEBFold and compared its performance with classical thermodynamics-based prediction tools on prepared TestSetα, TestSetβ, and TestSetγ. Nineteen prediction tools that were still publicly available were compared in this section: (1) tools from single-sequence thermodynamics models: RNAfold,²⁵ RNALfold,⁴⁴ Fold,²⁹ MaxExpect,²⁶ RNAprob,²⁸ RME (both the PARS-model and the DMS-model),⁴⁵ PKNOTS,⁴⁶ IPknot,⁶ HotKnots,⁴⁷ IterativeHFold,⁴⁸ ShapeKnots,⁴⁹ and ProbKnot;²⁷ (2) tools based on homologous sequence alignment and folding: RNAalifold,¹⁴ TurboFold,²⁹ SPARSE,⁵⁰ aliFreeFold,¹⁵ LocARNA,⁵¹ comRNA,⁵² and MXSCARNA.⁵³ Algorithms that fail to provide executable codes were not included in this comparison. Default parameters suggested by the authors were adopted in these tools for a fair comparison. For IterativeHFold, we set RNAfold as the restriction structure generator for it. For tools that require multiple sequence alignment results, sequences in the same Rfam family were taken. In some Rfam families, the sequences are nearly identical, even though they were collected from different species. These nearly identical sequences within the same family were used in tools that required homologous alignments if no other homologous sequence was available for that family. If no sequence is available in the same RNA structural family of the given sequence, the required homologues were searched against the RNAcentral database⁵⁴ using BLAST as suggested in our previous work.³

The results of different tools on TestSetα, TestSetβ, and TestSetγ are summarized in Table 1. In Table 1, the performance of each tool was summarized by its median rank and overall median F1 score on the three test sets. Among single-sequence input tools [RNAfold, RNALfold, Fold, MaxExpect, RNAprob, RME (both the PARS-model and the DMS-model), PKNOTS, IPknot, HotKnots, IterativeHFold, ShapeKnots, and ProbKnot], DEBFold achieves the overall top-one median F1 score rank on the three test sets. The improvement is supposed to be from the integration of various existing thermal knowledge followed by deep learning boosting. Compared with tools that utilize homologous sequences (RNAalifold, TurboFold, SPARSE, aliFreeFold, LocARNA, comRNA, and MXSCARNA), DEBFold still has a better overall median F1 score rank. Combing these comparisons, DEBFold provides cutting-edge performance over existing tools on the test sets. Notice that in TestSetα, tools that rely on sequence alignment results generally perform better than single-sequence input tools since the RNA sequences from Rfam 12.2 mostly come with good family-wise homologous sequences. This advantage may diminish in many real-world use cases, such as in TestSetγ. DEBFold overcomes this restriction by using a deep learning model that integrates the results of many single-sequence structure prediction tools. It is worthwhile noting that DEBFold also boosts the prediction F1 score over the 6 predictions (RNAfold, IPknot, MaxExpect, ProbKnot, RNAprob, and Fold) integrated in the input encoding tensor of the deep network. In summary, DEBFold is concluded to outperform existing classical structure prediction tools with less input requirements.

Recently, more and more deep-learning-based structure prediction tools were developed to help predict RNA secondary structures. We also sought to compare DEBFold to these tools. To the best of our knowledge, existing deep-learning-based RNA secondary structure prediction tools include the following: SPOT-RNA,²¹ SPOT-RNA2,⁵⁵ MXfold2,²² GCNFold,⁵⁶ UFold,⁵⁷ REDFold,⁵⁸ and e2eFold.⁵⁹ Since SPOT-RNA and SPOT-RNA2 released by the authors do not encompass any retraining codes and the tools themselves incorporate the whole bpRNA-1m data set in its model training process, independent test sets other than TestSetα, which was prepared from the bpRNA-1m Rfam 12.2 part, should be used to ensure fair comparison between DEBFold and SPOT-RNA/SPOT-RNA2. For this purpose, we downloaded the bpRNA-new data set from the MXfold2 work and prepared TestSetβ based on Rfam 14.4 (excluding those already in Rfam 12.2) for fairly evaluating the performance of DEBFold with existing deep-learning-based prediction tools. In TestSetβ, both the length bias and data contamination issues were eliminated. Details of the TestSetβ preparation steps can be found in the “Contamination-Free Family-Wise Independent Test Set for Evaluating Existing Deep-Learning-Based Structure Prediction Tools” section. Another PDB-derived independent test set, TestSetγ, was also prepared for this purpose. The preparation steps of TestSetγ can be found in the “PDB-Derived Source-Independent Test Set” section.

The final test results of DEBFold and the existing deep-learning-based prediction tools on TestSetα, TestSetβ, and TestSetγ are listed in Table 2. For SPOT-RNA, SPOT-RNA2, and UFold, TestSetα results were not reported due to complete data contamination caused by using the whole bpRNA-1m data set in the training process. REDFold has potential slightly optimistic evaluation results on TestSetα and TestSetβ since a small portion (122 families, whose details were not clearly reported by the original paper) of the Rfam 14.4 families were included in the training process. MXfold2 also included a small portion of the Rfam data (22 Rfam families from Rfam 10.0 and 151 sequences from the S-151Rfam data set, whose details were not clearly reported by the original paper), which might result in potential slightly optimistic results for TestSetα and TestSetβ. We calculated the median performance ranking and median F1 score of each tool on suitable test sets in the comparison. For tools that are suitable to be evaluated on all three test sets (DEBFold, MXfold2, REDFold, GCNfold, and e2efold), DEBFold has superior median F1 score results. Compared with SPOT-RNA, UFold, and SPOT-RNA2, DEBFold shows the top overall median F1 score rank. Among these tools, SPOT-RNA was first pretrained on the bpRNA-1m data set. Then, these pretrained models were fine-tuned and aggregated on the PDB data set. In other words, both the RNA structural families from bpRNA-1m and the PDB-identified RNA sequences were considered in the tool. Because of the final fine-tuning process, the final published SPOT-RNA tool somewhat favors 3D-structure-derived RNA base pairings. However, in PDB-derived TestSetγ, DEBFold still outperforms SPOT-RNA. Therefore, we conclude that DEBFold can provide better results than SPOT-RNA even on the 3D-structure-derived base pairs. It is worth noting that SPOT-RNA2 can achieve good results by considering the homologous sequence information. Nonetheless, SPOT-RNA2 only gets the same median performance rank as DEBFold while demanding more than 7500 times longer execution time than DEBFold (see Table 5). In view of the overall consideration of prediction efficiency and accuracy, DEBFold still has better results than SPOT-RNA2. Overall, the performance improvements in DEBFold are supposed to be attributed to a novel two-stage pipeline that utilizes a 1D structure representation and deep network integration of various thermal models. These analyses indicate that DEBFold can achieve state-of-the-art generalization structure prediction performance over that of currently available deep-learning-based prediction tools.

The second step of DEBFold feeds the computed structure location folding probabilities into thermodynamically constrained optimization algorithms to obtain the final RNA secondary structure prediction. Previous sections show that DEBFold can have a better structure prediction over the integrated constitutional prediction results by using more accurate folding scores as the optimization soft constraints. We next evaluated the robustness of DEBFold against different constrained optimization algorithms. In DEBFold, Fold was the adopted algorithm for the thermodynamically constrained optimization of the final predicted structures. Besides Fold, RNAfold and ShapeKnots can also help perform thermodynamics-constrained optimization based on soft constraints. These three constrained optimization algorithms were separately combined with DEBFold Stage I (DEBFold-Fold, DEBFold-RNAfold, and DEBFold-ShapeKnots) and then evaluated on TestSetα, TestSetβ, and TestSetγ. The comparison results are summarized in Table 3. As shown in Table 3, DEBFold achieves nearly identical median F1 score results (within one percent) on the three test sets when different thermodynamics-constrained optimization tools are used. From this comparison, it is suggested that DEBFold can provide accurate folding scores that help boost the final structure prediction and is robust against different thermodynamics-constrained optimization tools.

In DEBFold Stage 1, the base locations where the structure pairings occur are labeled. The higher performance of the labeling process usually leads to better constrained optimization results. We compared the structure location labeling performance of DEBFold Stage 1 with the simple consensus model and another existing technique called GRASP.⁶⁰ In DEBFold Stage 1, the probability of a base location being involved in the final structure pairings is predicted by a deep convolutional network. In previous research, Ke et al. proposed using the XGBoost model on the flattened one-hot encoding tensor of a sequence window to get the pairing probability of an RNA base location and implemented this method as a tool called GRASP.⁶⁰ In both DEBFold Stage 1 and GRASP, a fair threshold of 0.5 was adopted when the existence of base pairings was labeled for different RNA base locations. On the other hand, the simple consensus model labels the existence of the base pairing for an RNA base location if all 6 tools (RNAfold,²⁵ IPknot,⁶ MaxExpect,²⁶ ProbKnot,²⁷ RNAProb,²⁸ and Fold²⁹) agree on the pairing existence. In order to evaluate the location labeling results, we resort to the confusion matrix generated from the labeling results of a specific tool for each RNA. For structure location labeling in each RNA, TP is the number of correctly labeled pairing locations, FN represents the number of missed pairing locations, FP shows the number of nonpairing locations mistakenly marked to be paired, and TN counts the number of correctly labeled nonpairing locations. Based on the defined confusion matrix, the recall, precision, and F1 score values for the structure location labeling problem are defined similarly to those in eqs 15 and 16. We calculated the evaluation metrics for each RNA on the three prepared test sets (TestSetα, TestSetβ, and TestSetγ) and obtained the median F1 score of all of the RNAs in each test set. As summarized in Table 4, DEBFold Stage 1 achieves F1 scores superior to those of both the consensus model and GRASP on all three test sets. By this comparison, we conclude that DEBFold Stage 1 is better than the simple consensus model and the existing tool GRASP in labeling the existence of the base pairings for RNA base locations.

Due to the advancement of sequencing technology, lots of RNA sequences are found, leading to the need for large-scale structure investigation directly on the sequences. To fulfill this need, the execution time required for each prediction tool should also be inspected in addition to structure prediction accuracy. We evaluated and compared the execution time of DEBFold and different structure prediction tools on the collected 147 PDB-derived sequences using a workstation with Intel i7 CPU cores and 128GB RAM. The results are summarized in Table 5. On these 147 RNA sequences, most single-input folding algorithms used less than a minute to process the predictions. Since DEBFold integrates the outcomes of six different single-input folding prediction results, it took around 1 min to finish all predictions. For algorithms that consider homologous sequences (RNAalifold,¹⁴ TurboFold,²⁹ SPARSE,⁵⁰ aliFreeFold,¹⁵ LocARNA,⁵¹ comRNA,⁵² and MXSCARNA⁵³), around an additional 10 min was demanded to get the homologous sequences using BLAST. Therefore, DEBFold has better time efficiency than these homology-based tools and achieves better performance. In deep-learning-based predictions, SPOT-RNA required about 3 times more execution time than did DEBFold. Among all tools, SPOT-RNA2 had very low execution efficiency and needed more than 142 h, which is more than 7500 times as long as the execution time of DEBFold, to predict merely the structures of the 15 TestSetγ sequences selected by stratified sampling from all 147 RNAs. More execution time is required to finish all 147 RNAs. The high execution time of SPOT-RNA2 suggests that SPOT-RNA2 is probably not ready for a large-scale RNA structure investigation. Moreover, SPOT-RNA2 requires more than 2TB of disk space to store the information needed for its algorithm. In summary, DEBFold not only achieves state-of-the-art prediction accuracy but also retains an acceptable execution time for large-scale RNA sequence investigation.

In this work, we developed a deep-learning model that can integrate the prediction results from different single-sequence folding tools to more accurately predict RNA secondary structures. DEBFold is a two-stage pipeline for RNA secondary structure prediction. In the first stage, DEBFold identifies the locations where structural pairings occur. Then, the pairing probabilities for the base locations of a given sequence are transformed into the SHAPE-like score tensor. Based on the SHAPE-like score tensor, the structure prediction bearing the minimum free energy under these soft constraints is given by Fold. In RNA structure prediction, pseudoknots (intertwined helices on a plane) and noncanonical pairings (the base pairs formed by hydrogen bonding differing from the patterns of standard Watson–Crick rules⁶¹) present significant challenges for prediction algorithms. Since DEBFold Stage 1 does not distinguish pseudoknot, canonical, and noncanonical RNA pairings, these pairings are all potentially considered in this stage. It is possible to adopt ShapeKnots or other minimum free energy structure prediction optimizers in DEBFold Stage 2 if these challenging types of base pairings are to be considered in the final optimization process. However, since there are insufficient noncanonical or pseudoknot pairings in the ground-truth data set, DEBFold currently mainly focuses on canonical base pairings.

Previously, Szikszai et al.²³ pointed out that deep-learning models for RNA secondary structure prediction can easily overfit a data set if a sequence-based cross-validation and a test process are adopted. This type of overfitting occurs even when a massive number of sequences from only a few families are provided. Szikszai et al. argued that the chief cause for the phenomenon lies within the structural similarities among RNA sequences in the same family. Based on their research, we adopted a family-wise approach in developing DEBFold and further prepared test sets that are more suitable for evaluating this problem in existing deep-learning-based prediction tools. Although DEBFold overcomes the overfitting problem and provides state-of-the-art performance, inherent limitations remain. Currently, the known RNA families are far smaller than the known sequences in the community. Despite the large number of RNAs with known sequences available, most of them belong to the same structural families. Because of the limited RNA structural families available (only 2606 different families in Rfam 14.4), the performance of single-sequence structure prediction using deep learning is largely restricted. The limited RNA structural families also constrain the prediction accuracy of long noncoding RNAs (lncRNAs) since the overfitting problem boils down to the lack of structural family diversity. Currently, popular data sets contain scarce numbers of structural families for long RNAs, inhibiting the structural prediction of lncRNAs. The current version of DEBFold is suitable only for sequences up to 512 bps. Although some of the existing deep-learning-based tools claim to be able to deal with lncRNAs, the moderate performance observed in Table 2 indicates that the originally claimed prediction accuracy may still need more attention. Notice that DEBFold performs better than classical folding algorithms that incorporate multiple-sequence alignments. Therefore, whether multiple-sequence alignments can better help suggest structure modules in the DEBFold pipeline requires further study.