AlignmentFold and AlignmentPartition: Improving the align-then-fold approach for RNA secondary structure prediction

AlignmentFold and AlignmentPartition: Improving the align-then-fold approach for RNA secondary structure prediction | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results AlignmentFold and AlignmentPartition: Improving the align-then-fold approach for RNA secondary structure prediction Abhinav Mittal , View ORCID Profile David H. Mathews doi: https://doi.org/10.1101/2025.07.23.666478 Abhinav Mittal 1 Department of Biochemistry and Biophysics and Center for RNA Biology, University of Rochester Medical Center , Rochester, New York, 14642, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site David H. Mathews 1 Department of Biochemistry and Biophysics and Center for RNA Biology, University of Rochester Medical Center , Rochester, New York, 14642, USA 2 Department of Biostatistics and Computational Biology, University of Rochester Medical Center , Rochester, New York, 14642, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for David H. Mathews For correspondence: david_mathews{at}urmc.rochester.edu Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF ABSTRACT The combination of free energy minimization with sequence covariation analysis is a popular approach for predicting conserved secondary structure using a fixed alignment of homologous RNA sequences. In this work we revisit this approach by developing the new AlignmentFold and AlignmentPartition programs. AlignmentFold and AlignmentPartition predict consensus minimum free energy structures and base pairing probabilities, respectively. We determined new nearest neighbor thermodynamic parameters for gaps and non-canonical base pairs that align with base pairs. These frequently occur for structure prediction with a fixed input alignment. We assessed the impact of alignment prediction method and size on the prediction accuracy. We also assessed the improvement in prediction accuracy conferred by sequence covariation analysis. In contrast to previous solutions, we do not include covariation in the scoring of structure quality and achieve prediction accuracy as good as other tools. AlignmentFold and AlignmentPartition are freely available as part of the RNAstructure software package at https://rna.urmc.rochester.edu/RNAstructure.html . INTRODUCTION RNA sequences play diverse roles in cellular processes ranging from gene expression to anti-viral defense ( Cech and Steitz 2014 ). Consequently, they emerged as drug targets, drugs, and vaccines ( Damase et al. 2021 ). While new classes of non-coding RNAs (ncRNAs) with their functions continue to come to light ( Flynn et al. 2021 ), a better understanding of RNA structure and folding is expected to provide further insights ( Spitale and Incarnato 2022 ). The set of canonical base pairs, including Watson-Crick-Franklin AU and GC base pairs and wobble GU base pairs, constitute the secondary structure for an RNA. For functional RNAs, the secondary structure tends to remain evolutionarily conserved for the function to be conserved. Improved RNA secondary structure prediction is key to discovering new functional RNAs ( Gao et al. 2021 ; Fu et al. 2015 ; Gruber et al. 2010 ; Zhang et al. 2023 ) and offers a better understanding of the known ones ( Hagey et al. 2022 ; Li et al. 2021 ; Szutkowska et al. 2022 ). Comparative sequence analysis is the gold standard for RNA secondary structure determination in the absence of a solved 3D structure ( Gutell et al. 2002 ). However, the availability of a large number of homologous sequences is a prerequisite for RNA secondary structure determination. Base pairs are identified by finding covariations in an alignment that is refined to reflect the structure. These covariations, for example an AU base pair in one sequence that replaces a GC base pair in another sequence, provide evidence for conserved pairs. Covariations occur frequently across evolution because RNA structure is conserved although sequence is not. The process of comparative analysis is semi-automated, but still requires considerable manual effort to deduce covariation and iteratively refine the predicted structure ( Gutell 2014 ). One important reason for the difficulty in automating this process is that the covariations that reveal the conserved structure complicate the alignment of the sequences based on nucleotide identity alone. The nearest neighbor thermodynamic model (a.k.a. Turner Rules) can be used to estimate the Gibbs free energy change of folding for a given secondary structure (ΔG°)( Xia et al. 1998 ; Mathews et al. 2004 ; Mittal et al. 2024 ). Dynamic programming algorithms using this model can predict the structure for a given sequence with reasonable accuracy by estimating the lowest free energy structure ( Mathews et al. 2010 ). Dynamic programing algorithms are also used to estimate base pairing probabilities with partition function calculations ( McCaskill 1990 ). The pairing probabilities offer a more nuanced view of the secondary structure, especially for sequences that sample multiple secondary structures ( Mathews 2006 ). The base pairing probabilities also provide information about the confidence of predicted pairs ( Mathews 2004 ). Algorithms have been developed to automate comparative sequence analysis by applying structure prediction to a set of homologous sequences to predict the consensus structure that is the conserved secondary structure common to all sequences. ‘Align-then-fold’ has been one of the popular approaches to determine the consensus secondary structure using a set of homologous sequences ( Bernhart et al. 2008 ; Hofacker et al. 2002 ; Knudsen and Hein 1999 , 2003; Rivas 2020 ; Seeman et al. 2008 ). In this approach, a fixed input alignment is provided by the user. The thermodynamic nearest neighbor model along with covariation analysis is used to predict structures where base pairs are common for all sequences in the alignment. Dynamic programming is used to search for the best scoring consensus structure or the consensus base pairing probabilities. These approaches are popular because they work well when the sequences are similar and the computational cost is lower than approaches that simultaneously predict structure and the sequence alignment. The structure prediction accuracy is improved over single sequence structure prediction because structures with conserved base pairs are scored more favorably by the nearest neighbor model and the covariation terms. As summarized in Figure 1 , we report two new programs: AlignmentFold and AlignmentPartition. These take the align-then-fold approach to determine the consensus minimum free energy structure and base pairing probabilities, respectively. For this work, we developed principled free energy change parameters that account for gaps and non-canonical base pairing that arise during simultaneous folding of sequences ( Mittal et al. 2024 ). We take advantage of the ability of RNAstructure to use an extended nucleotide alphabet, including the gap symbol, to use these new parameters ( Kierzek et al. 2022 ). Download figure Open in new tab Figure 1. AlignmentFold and AlignmentPartition use an extended nearest neighbor parameter set to predict a consensus secondary structure from a fixed alignment of homologous RNAs. We derived principled parameters to include the gap symbol. In this work, we benchmarked structure prediction methods using manually curated alignments and predicted alignments. The predicted alignments test expected performance of these programs before the structure of an RNA family is known; the manually curated alignments are informed by structure. We benchmarked our new programs against popular available programs, RNAalifold ( Hofacker et al. 2002 ; Bernhart et al. 2008 ), PETfold ( Seeman et al. 2008 ) and CaCoFold ( Rivas 2020 ) using 5, 10, and 20 sequence sets designed to mimic the early stages of structure determination for newly discovered ncRNAs. In contrast to our expectations, we found that including terms that account for covariation analysis brings only modest improvement at best in prediction accuracy. An increase in the number of sequences in the input might bring marginal improvement in accuracy. We also tested the accuracy when using alignments generated by several sequence alignment programs. The alignments generated by MAFFT X-INS-I ( Katoh and Toh 2008 ), which considers structure when aligning the sequences, lead to more accurate consensus structure prediction. RESULTS AND DISCUSSION Folding Model AlignmentFold and AlignmentPartition find the lowest free energy secondary structure and partition function, respectively, by mapping predicted base pairs to all the sequences in a sequence alignment. The partition function estimates the pairing probabilities for base pairs in the conserved structure. The programs calculate the mean folding free energy change across the sequences for a given structure. The alignment gap symbol is treated as a nucleotide in an extended alphabet of bases. Base pairs other than AU, GC and, GU occur when the structure is mapped to the individual sequences in the alignment. These include non-canonical pairs like “GG” and gapped pairs like “-A” ( Figure 2A ). We developed nearest neighbor parameters to handle these cases, following the 2004 Turner Rules ( Mathews et al. 2004 ; Mittal et al. 2024 ). This takes advantage of the ability to use extended alphabets of nucleotides ( Kierzek et al. 2022 ). We also modified the single sequence folding software in the RNAstructure package to work with multiple sequence alignments. This provides suboptimal structure prediction (alternative low free energy structure prediction), simplifies maintenance and simplifies the incorporation of new features. The details of the implementation are provided in the METHODS section. Download figure Open in new tab Figure 2. Appearances of alignment gaps and non-canonical pairs. Panel A shows a hairpin-stem-loop structure mapped to individual sequences in an alignment. This mapping can result in gaps appearing in base paired positions, for example the -A base pair that appears in the first sequence. The mappping results in a GG base pair in the second sequence. Panel B illustrates that gaps are removed from hairpin loops (and also internal loops) before free energy estimation. In this example, the six-membered hairpin loop contains two gaps that are removed to make a tetraloop. Parameter Optimization To optimize four adjustable parameters (introduced below), we adopted a set of training RNA families that are separate from families used to test performance ( Tan et al. 2017 ; Szikszai et al. 2022 ; Lu et al. 2009 ; Liu et al. 2010 ; Rivas et al. 2012 ). Consistent with our prior work, we used 5S ribosomal RNA (Eubacteria subfamily), group I intron (IC1 subfamily), tmRNA, and tRNA families for training and 16S rRNA (Alphaproteobacteria subfamily), SRP RNA (Protozoan subfamily), RNase P RNA (bacterial type A subfamily) and telomerase RNA families for testing, all of which were drawn from our RNAStralign data set ( Tan et al. 2017 ). For determining a consensus secondary structure from a fixed alignment, a criterion to determine positions that can base pair is required. For the alignment, positions are columns in the alignment, represented by nucleotides or gaps in each sequence. For AlignmentFold, we used a simple base pairing criterion where positions (i,j) can base pair if the fraction of sequences that can form a canonical base pair exceeds a cutoff. A higher, more restrictive cutoff leads to a higher positive predictive value (PPV; fraction of predicted pairs that are correct), lower sensitivity (fraction of known pairs correctly predicted), lower F 1 score (the harmonic mean of sensitivity and PPV), and substantially improved run-time performance (Supplementary Table S1). To strike a balance between sensitivity, specificity, and runtime, we selected 0.3 as the default cutoff. A less restrictive base pairing cutoff can be chosen by the user to allow a greater number of pairs of positions in alignment to be eligible for base pairing. This is more resilient to alignment errors, but substantially increases the run time due to the expanded consensus secondary structure space. The maximum size of an internal loop is set to 30 unpaired bases by default for single sequence folding in the RNAstructure package. For AlignmentFold, this was increased to 40 to accommodate the plausible increase in the size of internal loops due to the presence of gaps in the alignment. Increasing the maximum size of internal loops leads to improved sensitivity (Supplementary Table S2). It also leads to higher runtime due to the expanded consensus secondary structure space. For AlignmentPartition, we used the same parameters, i.e. 0.3 as the base pairing cutoff and 40 as the maximum size of internal loop. This is important because the estimated pair probabilities reflect the same landscape of possible structures that free energy minimization structure prediction considers. The ProbKnot ( Bellaousov and Mathews 2010 ; Zhang et al. 2019 ) and MaxExpect ( Lu et al. 2009 ) tools in the RNAstructure package predict structures from base pairing probabilities. We adapted these to use the output of AlignmentPartition to predict conserved structures. ProbKnot is capable of predicting pseudoknots ( Liu et al. 2010 ), although AlignmentFold and MaxExpect are not. We found 4.0 to be the optimal value of the gamma parameter for MaxExpect on the training set (Supplementary Table S3) and 0.1 to be optimal threshold for ProbKnot (Supplementary Table S4). Accuracy Benchmarks We benchmarked AlignmentFold, AlignmentPartition (generating structures with MaxExpect or ProbKnot), RNAalifold ( Hofacker et al. 2002 ; Bernhart et al. 2008 ), PETfold ( Seeman et al. 2008 ) and CaCoFold ( Rivas 2020 ) on alignments estimated using Clustal Omega ( Sievers et al. 2011 ), MAFFT X-INS-I ( Katoh and Toh 2008 ) and TurboFold II ( Tan et al. 2017 ) and also on the manually curated alignments. It is important to benchmark align-then-fold structure prediction tools using predicted alignments because the accepted alignments were generated using the structure for guidance ( Malik et al. 2024 ). RNA secondary structure prediction from an alignment of homologous RNAs presents a chicken-and-egg problem where we need a high-quality sequence alignment to arrive at accurate prediction of conserved structure and vice-versa ( Seetin and Mathews 2012 ). Also, we used alignments of 5, 10, and 20 sequences to represent the early stages of structure determination when a limited number of homologous RNAs are available for a newly discovered RNA family. The best performance for structure prediction accuracy was obtained using curated alignments, as expected because those are informed by the conserved structure (Supplementary Table S5). The best structure prediction accuracies for predicted alignments were obtained using MAFFT X-INS-i ( Figure 3 and Supplementary Table S5). Prediction accuracy on sequence alignments generated using TurboFold II was similar (within 3%; Supplementary Tables S5-S9), except for SRP RNA where MAFFT X-INS-i is substantially better (Supplementary Table S8). Structure prediction accuracy is worst with alignments generated by Clustal Omega, which estimates alignments without considering the structure and relies solely on sequence identity. Also, PETfold performs substantially better (> 10%) on MAFFT X-INS-i-aligned RNase P sequences. Both MAFFT X-INS-i and TurboFold II consider secondary structure information generated through thermodynamic nearest neighbor-based computational approaches. Download figure Open in new tab Figure 3. F 1 score for RNAalifold (default), RNAalifold (RIBOSUM), RNAalifold (−-cfactor=0), PETfold, PETfold (−p 1), AlignmentFold, AlignmentProbability MaxExpect (AP-MaxExpect), AlignmentProbability ProbKnot (AP-ProbKnot), and TurboFold II on 200 sets of 5, 10 and 20 sequences each aligned using MAFFT X-INS-I except TurboFold II, which takes unaligned sequences as input. Plotted uncertainty is ± the standard deviation. The average F 1 score of single sequence folding with RNAstructure Fold (lowest free energy structure prediction) has been plotted using a dotted line. F 1 scores which are significantly different (p < 0.05) than AP-ProbKnot are marked by *. The accuracies of most align-then-fold tools are similar, although there are idiosyncrasies observed across families. AlignmentFold, AlignmentPartition (ProbKnot and MaxExpect), PETfold and RNAalifold offer better prediction accuracy as measured by F 1 score than single sequence structure prediction by Fold for 16S rRNA, RNase P and telomerase RNA ( Figure 3 ). PPV and sensitivity are available as Supplementary Figures S1 and S2. For SRP RNA, these tools can offer superior PPV but substantially poorer sensitivity with structurally informed alignments (Supplementary Table S10). Except for CaCoFold, improvement in prediction accuracy in terms of F 1 score with increase in the number of sequences is not guaranteed (Supplementary Table S5). The AlignmentPartition-ProbKnot structure predictions outperform all other methods including PETfold, RNAalifold, and CaCoFold for SRP RNA and telomerase RNA across all alignment modalities (Supplementary Figures S1 and S2 and Supplementary Tables S8 and S9). AlignmentPartition-MaxExpect predictions are more accurate than AlignmentFold and AlignmentPartition-ProbKnot for 16S rRNA ( Figure 3 and Supplementary Table S6). ‘Align-and-fold’ algorithms, unlike ‘align-then-fold’ algorithms, take unaligned homologous sequences as input and return structures specific to each sequence in input along with an alignment. These algorithms tend to have better prediction accuracy than align-then-fold algorithms ( Tan et al. 2017 ). TurboFold II is one such ‘align-and-fold’ algorithm. TurboFold II has better prediction accuracy than PETfold, RNAalifold, and AlignmentFold when input to these algorithms is an estimated alignment (Supplementary Table S11). The only exception being the 16S rRNA family where PETfold outperforms TurboFold II, and other align-then-fold methods also offer comparable accuracy to TurboFold II (Supplementary Tables S6 and S11). This is consistent with prior benchmarks where it was noted that the pairwise sequence identities of the 16S rRNA sequences is higher than that of other families ( Tan et al. 2017 )(Supplementary Figure S3). Comparison of Methods In this work, we benchmarked RNAalifold version 2.4. It by default uses covariation-based scoring as described in 2002 ( Hofacker et al. 2002 ). It uses a generalized formula to calculate the contribution of the covariation to the final score. RNAalifold with the -r option uses specialized RIBOSUM-like matrices for calculating the contribution of covariation to the final score as described in 2008 ( Bernhart et al. 2008 ). It adapts to the levels of sequence identity in the alignment in that it stabilizes the formation of base pairs with compensatory mutations more generously when the level of sequence identity is high. The contribution of the covariation score to the total score often exceeds 50% (Supplementary Table S12). RNAalifold allows the user to tune the contribution of covariation-based scoring vis-à-vis thermodynamics for consensus structure prediction. A user can set it with ‘--cfactor’ option. By default, it is set to 1. We set it to 0 to benchmark RNAalifold without covariation-based scoring. Typically, RNAalifold with RIBOSUM-like scoring (−r flag) performs the best for the RNAalifold methods. It is closely followed by RNAalifold without covariation-based scoring, which in turn performs better than RNAalifold with default settings ( Figure 4 and Supplementary Table S5). Download figure Open in new tab Figure 4. F 1 score for RNAalifold (−-cfactor=0), RNAalifold (default), and RNAalifold (RIBOSUM) on 200 sets of 5, 10 and 20 sequences each aligned using MAFFT X-INS-i. Uncertainties are plotted as ± the standard deviation. F 1 scores which are significantly different (p < 0.05) than RNAalifold (−-cfactor=0) are marked by *. The F 1 accuracy of RNAalifold (default) tends to decrease with increasing number of sequences in the input alignment (Supplementary Table S5). There are some exceptions to this trend. For example, RNAalifold (default) outperforms RNAalifold with no covariation-based scoring and RNAalifold with RIBOSUM-like scoring on the 16S rRNA family across all alignment methods (Supplementary Table S6). This is particularly interesting because RIBOSUM-like matrices were derived from analysis of 16S rRNA sequences ( Bernhart et al. 2008 ). PETfold provides a maximum expected accuracy structure by optimally combining pairing probabilities determined using RNAfold and Pfold ( Knudsen and Hein 2003 ). While RNAfold relies on thermodynamic parameters to find the minimum free energy structure, Pfold draws from an explicit evolutionary model and SCFG-based probabilistic model for secondary structures. PETfold was optimized over curated alignments from Rfam.( Griffiths-Jones et al. 2005 ) The default parameter setting limits the influence of the evolutionary model on structure prediction and emphasizes the thermodynamics. PETfold removes columns with more than 25% gaps from the alignment before structure prediction and introduces a gap in their positions in the final consensus structure. It calculates “expected overlap” for a given consensus structure over all the sequences individually and avoids the issue of non-canonical pairing that arises during simultaneous folding of the sequences. PETfold with the -p option set at 1 relies solely on base pairing probabilities generated with RNAfold and does not consider probabilities generated by Pfold. Therefore, PETfold with -p 1 option does not consider covariation information. PETfold with default settings tends to perform slightly (∼1%) better in F 1 score than PETfold with the -p 1 option (Supplementary Table S5; MAFFT X-INS-i alignments). However, there is an exception to this trend. PETfold without covariation information outperforms PETfold with default settings as well as all other methods on 16S rRNA alignments of size 5, 10, and 20 sequences across all alignment methods (Supplementary Table S6). But PETfold with default settings again performs better on large 16S rRNA alignments of 100 sequences ( Figure 5 and Supplementary S13). PETfold with default settings amongst all methods offers the best F 1 performance for 100 RNase P RNA sequences aligned with MAFFT X-INS-i or Clustal Omega (Supplementary Table S13 and Figure 5 ). PETfold offers suboptimal consensus structure prediction as well. Download figure Open in new tab Figure 5. F 1 score for RNAalifold (default), RNAalifold (RIBOSUM), RNAalifold (−-cfactor=0), PETfold, PETfold (−p 1), AlignmentFold, AlignmentProbaility MaxExpect, AlignmentProbability ProbKnot, and CaCoFold on 10 sets of 100 sequences each aligned using MAFFT X-INS-i. CaCoFold (Cascade variation/covariation Constrained Folding algorithm) takes a three-fold approach to consensus structure prediction ( Rivas 2020 ). First, it identifies pairs of positions in alignment for which covariation is statistically significant when considering phylogeny. Second, it identifies pairs of positions for which there is negative evidence of base pairing. These positions vary without significant covariation with other positions and are therefore forbidden from base pairing. Finally, folding is performed with probabilistic algorithms in a way that incorporates the above information about base pairing into one structure. Other methods benchmarked here consistently outperform CaCoFold on estimated alignments. CaCoFold, however, offers highest prediction accuracy for the curated alignment of RNase P RNA of 20 and 100 sequences (S7 and S13, Supplementary Information). These results suggest that CaCoFold requires a larger set of sequences to identify significant covariations. CaCoFold can predict pseudoknots like AlignmentPartition-ProbKnot and in contrast to RNAalifold, PETfold, AlignmentFold, and AlignmentPartition MaxExpect, which cannot predict pseudoknots. Time Benchmarks The runtime complexity of AlignmentFold is O(HN 3 ) where H is the number of sequences in the alignment and N is the length of alignment. This is comparable to O(N 3 ) of single sequence folding where N is the length of the sequence. Supplementary Table S14 shows the runtime performance of AlignmentFold and AlignmentPartition. CONCLUSIONS These benchmarks highlight the idiosyncratic performance landscape of align-then-fold methods across alignment prediction methods and sequence numbers. It also highlights the importance of benchmarking these tools on estimated alignments and across a range of number of sequences. An important finding of these benchmarks on up to 20 sequences is that there is limited improvement offered by the inclusion of covariation information in structure prediction, especially for estimated alignments. They suggest that the accuracy is a result of a complex interplay of multiple factors like base pairing acceptance criteria, handling of gaps, and underlying folding models in secondary structure prediction from the fixed alignment. This work brings forth the advantages and limitations of align-then-fold tools vis-à-vis single sequence folding and align-and-fold approaches like TurboFold II. AlignmentFold and AlignmentPartition offer an intuitive and transparent approach to consensus structure prediction. By using recursions previously developed for single sequence folding for processing multiple sequence alignments, it demonstrates an innovative approach towards software development for structure prediction. AlignmentFold and AlignmentPartition are incorporated in the RNAstructure software package, which is available freely under the GNU Public License at https://rna.urmc.rochester.edu . METHODS Software Implementation The existing recursions in the RNAstructure package for single sequence folding were used for AlignmentFold and AlignmentPartition ( Reuter and Mathews 2010 ). This was implemented with C++ function pointers. Function calls that calculate free energy increments for stacks and loops for a single sequence were replaced with function pointers. The pointers indicate functions that calculate these energies for multiple sequences when the input is a multiple sequence alignment and for a single sequence when the input is a single sequence. The energy model is therefore the full Turner 2004 model including the coaxial stacking of helices in multibranch and exterior loops ( Mittal et al. 2024 ; Ward et al. 2019 ). This approach will facilitate updates to the software, for example, as the free energy model undergoes further improvements ( Zuber et al. 2022 ). Improving run time performance For a given alignment, a 2D matrix stores information about whether the positions i,j in the alignment are eligible for base pairing as per the criteria. Another set of 2D matrices store information about the number of gaps between any two positions for each sequence in the alignment. For each sequence in the alignment, a version of the sequence without gaps is also stored. A 1-D vector corresponding to each sequence in the alignment serves as a map for base positions from a sequence in alignment with gaps to the gapless sequence. All the above-mentioned entities are initialized at the beginning of execution, and these allow rapid look up during run time leading to faster execution. Benchmarking run time performance For benchmarking run time performance, AlignmentFold and AlignmentPartition were compiled using GCC 4.8 and run on Red Hat Enterprise Linux 7 Maipo (x86_64) with an Intel Xeon CPU E5-2695 v2 @ 2.40GHz and 4 GB memory. Determining Structure Prediction Accuracy The accuracy of RNA secondary structure prediction is calculated in terms of PPV (positive predictive value), sensitivity, and F 1 score. PPV is the fraction of the predicted base pairs that also occur in the accepted structure. Sensitivity is the fraction of base pairs in the accepted structure that are correctly predicted. F 1 score is harmonic mean of sensitivity and PPV ( Mathews 2019 ). Predicted base pairs are considered correct if a nucleotide either on 5’ or 3’ end of the correct base is off by one position ( Mathews et al. 1999 ). For instance, a predicted base pair (i, j) is correct if base pair (i, j), (i±1, j) or (i, j±1) exists in the accepted structure. This is due to uncertainty in secondary structure determination by comparative sequence analysis ( Gutell et al. 2002 ) and thermodynamic fluctuation in local structures ( Fu et al. 2014 ; Znosko et al. 2002 ; Woodson and Crothers 1987 ). The scorer program of RNAstructure was used to calculate PPV and sensitivity. Benchmark Data The RNAStralign ( Tan et al. 2017 ) database was divided and used for determining hyperparameters and for accuracy benchmarking. RNAStralign is aggregated from expertly curated online databases of RNA secondary structure and sequence alignment. Sequences and their structures in RNAStralign are categorized into homologous families. The families included are 5S ribosomal RNA ( Szymanski et al. 2002 ), Group I intron ( Zhou et al. 2008 ), tmRNA ( Zwieb et al. 2003 ), tRNA ( Jühling et al. 2009 ), 16S ribosomal RNA ( Cannone et al. 2002 ), Signal Recognition Particle (SRP) RNA ( Rosenblad et al. 2003 ), RNase P RNA ( Brown 1999 ), and telomerase RNA ( Nawrocki et al. 2015 ). 16S rRNA, 5S rRNA, RNase P RNA, SRP RNA and group I intron families are further categorized into subfamilies in line with the original databases. Optimizing base pairing cutoff parameter In AlignmentFold and AlignmentPartition, base pairs are only permitted when the fraction of sequences that can form canonical base pairs exceeds a pairing cutoff. To determine the optimal base pairing cutoff, 250 groups of input sequences, comprising five, ten and twenty homolog sets, were randomly selected from the 5S ribosomal RNA (Eubacteria subfamily), group I intron (IC1 subfamily), tmRNA, and tRNA families of the RNAStralign database. These sequence sets were aligned using MAFFT (X-INS-i)( Katoh and Toh 2008 ). A consensus structure was determined for each of the alignments using AlignmentFold with varied base pairing cutoff. Cutoffs 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7 and 0.8 were benchmarked (Supplementary Table S1). The consensus structure was mapped back to each sequence in the alignment to determine the predicted structure for each sequence. Only AU, GC and GU base pairs were retained during mapping. Accuracy was calculated for each sequence individually. 0.3 was selected as the default base pairing cutoff for both AlignmentFold and AlignmentPartition. Optimizing maximum size of internal loop In AlignmentFold and AlignmentPartition, there is a maximum number of unpaired nucleotides for internal and bulge loops, which is adjustable. Secondary structure space is not explored for internal loops of greater size for structure prediction. To determine the default size, alignments were generated as described for optimizing the base pairing cutoff above. A consensus structure was determined for each of the alignments using AlignmentFold with varied maximum internal loop sizes. Sizes 30, 35, 40, 45, 50, 55, 60 and 80 were benchmarked (Supplementary Table S2). The consensus structure was mapped back to each sequence in the alignment. Only AU, GC and GU base pairs were retained during mapping. Accuracy was calculated for each sequence individually. 40 was selected as the default maximum internal loop size for both AlignmentFold and AlignmentPartition. Optimizing gamma for MaxExpect MaxExpect uses a partition function save file as input. This file contains the information needed for estimating base pairing probabilities ( Lu et al. 2009 ). MaxExpect infers a structure by maximizing the sum of the base-pairing and single stranded nucleotide probabilities. Pairing probabilities are weighted by factor γ. where P bp ( i,j )is the base pairing probability for nucleotides at positions i and j and P ss ( k ) is the probability of remaining single-stranded (unpaired) for nucleotide at position k . The sums are over all base pairs (BPs) and single stranded nucleotides (SS) in structure s . To determine the optimal gamma, alignments were generated from the training set as described above for optimizing the base pairing cutoff. Partition function save files were generated for alignments using AlignmentPartition. Maximum expected accuracy structures were generated for varied values of gamma: 0.5, 1, 2, 4, 5, 6, 7, and 8 (Supplementary Table S3). The structures were mapped back to each sequence in the alignment. Only AU, GC and GU base pairs were retained during mapping. Accuracy was calculated for each sequence individually. 4 was selected as the default γ for structure determination using MaxExpect. Optimizing threshold for ProbKnot ProbKnot uses a partition function save file as input and returns a structure with the most probable base pairs. Each base pair ( i,j ) in this structure has probability p ( i,j ) that is highest amongst the competing pairs( Zhang et al. 2019 ). Also, p ( i,j ) is greater than or equal to a threshold probability. To determine the default threshold, alignments were generated from the training set as described earlier. Structures were generated with varied values of threshold: 0.0, 0.1, 0.2, 0.3, 0.4, 0.6, and 0.8 (Supplementary Table S4). The structures were mapped back to each sequence in the alignment. Only AU, GC and GU base pairs were retained during mapping. Accuracy was calculated for each sequence individually. 0.1 was selected as optimal threshold. Isolated base pairs In single sequence folding, a base pair (i,j) is considered isolated if base pairs (i-1,j+1) and (i+1,j-1) cannot base pair. Single sequence folding in the RNAstructure package by default does not allow these isolated base pairs, but users have an option to allow them. In AlignmentFold and AlignmentPartition, base pairing between positions (i,j) is considered isolated if positions (i,j) fulfill the base pairing criterion stated above but neither (i+1,j-1) nor (i-1,j+1) positions do so. However, an exception is made and pairing is allowed for positions (i,j) when there exist positions (m,n) such that m can pair to n, i<m<n< j and all the intervening positions (i+1,i+2…m-2, m-1, n+1,n+2…j-2,j-1) have at least 50% nucleotides as gaps. Similarly, base pairing is allowed for (i,j) when there exist positions (m,n) such that m can pair to n, mj and all the intervening positions (m+1,m+2…i-2, i-1, j+1,j+2…n-2,n-1) have at least 50% nucleotides as gaps. Accuracy Benchmarks The command line options used for benchmarking each program are provided in Table 1 . View this table: View inline View popup Download powerpoint Table 1. Commands for methods of sequence alignment and structure prediction benchmarked in this work. For benchmarking, 200 groups comprising of five, ten and twenty homologous sequences were randomly selected with replacement from the 16S rRNA (Alphaproteobacteria subfamily), SRP RNA (Protozoan subfamily), RNase P RNA (bacterial type A subfamily) and telomerase RNA. For SRP RNA, sequences shorter than 200 nucleotides were excluded because their structures are not consistent with those of longer sequences ( Tan et al. 2017 ; Rosenblad et al. 2003 ). Similarly, for 16S rRNA sequences shorter than 1400 nucleotides were excluded. The benchmarks were therefore performed on families other than those used to train the pairing cutoff parameter ( Szikszai et al. 2022 ; Rivas et al. 2012 ; Lu et al. 2009 ). For benchmarks, the sequences were aligned using a multiple sequence alignment software or extracted from curated alignment and then consensus structure prediction was performed. The consensus structure was mapped back to each of the individual sequences and the accuracy determined. All methods were benchmarked on the same groups of sequences. As a control, single sequence folding was performed on each individual sequence using Fold (single sequence free energy minimization) and the accuracy calculated. Although a sequence can appear multiple times in sequence alignments, the accuracy of a sequence is averaged only once for the Fold calculations. Benchmarking was also performed on 10 groups of 100 homologous sequences randomly selected with replacement from 16S rRNA (Alphaproteobacteria subfamily) and RNase P RNA (bacterial type A subfamily). The benchmarking was limited to sequences aligned using MAFFT X-INS-i and those extracted from curated alignments. Significance testing To assess the statistical significance of the differences in sensitivity, PPV and F 1 , two-tailed, paired t-tests were performed using the ‘scipy.stats.ttest_rel’ function of SciPy v1.2.3 ( Virtanen et al. 2020 ). The type I error rate was set to 0.05. Thermodynamic parameters for gaps and base pairs other than AU, GC, and GU AlignmentFold maps the structure to all the sequences in the alignment while searching for the optimal structure. It treats gaps as a nucleotide. Base pairs other than AU, GC and, GU occur when the structure is mapped to the individual sequences in the alignment. We developed new nearest neighbor parameters to handle cases where the non-canonical pairs are predicted as base pairs or when gap symbols appear in pairs. In general, we average Turner 2004 parameters to estimate a free energy change for the resulting structure with the gap removed. This is intended to balance the accuracy of folding free energy estimates with computational efficiency. We used the extended alphabet feature of RNAstructure to account for gaps, “-” as an additional nucleotide ( Kierzek et al. 2022 ). In the RNAstructure package, all the files containing the energy parameters reside in the ‘data_tables’ directory. The files containing the extended parameter set start with ‘msa’. Gaps are removed from hairpin loops and internal loops before free energy evaluation in AlignmentFold (See Figure 2B ), following the strategy used in RNAalifold ( Bernhart et al. 2008 ). If the unpaired number of nucleotides in the hairpin loop are fewer than 3, which is forbidden in single sequence folding, we use 8 kcal/mol as the loop closure penalty. This approximates the 5.4 kcal/mol cost of closing a hairpin of 3 unpaired nucleotides, plus the average cost of opening a base pair stack (2.1 kcal/mol). For a bulge or an internal loop, if the number of unpaired nucleotides is zero after removal of gaps in the loop region, we treat the motif as a stack of base pairs. Helix stack parameters are specified in msa.stack.dg. When a gap-gap pair is followed by another gap-gap pair, a non-canonical pair, or a canonical pair, the stack is 0 kcal/mol. When a gap-nucleotide base pair is followed by a non-canonical pair, a canonical pair, or a gap-gap pair, 0.6 kcal/mol, the penalty used per nucleotide asymmetry in internal loops is used. When a gap-nucleotide base pair is followed by another gap-nucleotide base pair where the gaps appear in opposite strands, 0 kcal/mol is used. When a gap-nucleotide base pair if followed by another gap-nucleotide base pair where the gaps appear on the same strand, 1.2 kcal/mol is used because this is twice the asymmetry penalty. When a non-canonical base pair follows another non-canonical base pair, 0.2 kcal/mol is used. This is equivalent to free energy change incurred by adding two unpaired nucleotides to an internal loop. For stacks consisting of one non-canonical base pair and one canonical base pair, the free energy change was estimated by linear regression from RNA tables for 1×1 internal loops and terminal mismatches. We treated the 1×1 internal loop table (rna.int11.dg) as combinations of two stacks consisting of one non-canonical base pair and one canonical base pair. Similarly, we treated the terminal mismatches (rna.tstack.dg) as a stack of a pair and a non-canonical pair. With these data, linear regression estimated the sequence dependent free energy contribution of the stack of a canonical pair followed by a non-canonical pair. These stacking parameters are also used for the coaxial stacking of adjacent helices with no intervening mismatch as contained in msa.coaxial.dg ( Walter and Turner 1994 ). For dangling ends (msa.dangle.dg), when the base pair adjacent to the dangling base is non-canonical or contains one or more gaps, 0 kcal/mol is used. Also, when the dangling base is a gap, 0 kcal/mol is used. Terminal mismatch stacks are specified in msa.tstack.dg for exterior loops, in msa.tstackm.dg for multibranch loops, and in msa.tstackcoax for coaxial stacks of adjacent helices with an intervening mismatch ( Tyagi and Mathews 2007 ; Kim et al. 1996 ; Walter et al. 1994 ). Following the 2004 nearest neighbor parameters, these sequence-dependent terminal mismatches values are identical in these different structural contexts ( Mittal et al. 2024 ). When the terminal base pair is non-canonical pair or contains one or two gaps, 0 kcal/mol is used. If the terminal base pair is canonical and one of the mismatching bases is a gap, then the motif is treated as a dangling end and corresponding free energy from rna.dangle.dg is used. When both nucleotides in the terminal mismatch are gaps, 0 kcal/mol is used. Internal loops of size 1×n (n>2), 2×3, or i×j (i,j >2) can be stabilized by first mismatches at each end (msa.tstacki.dg, msa.tstack1n.dg, and msa.tstack23.dg). When non-canonical base pairing or gaps are present in first mismatches, 0 kcal/mol is used. When one or both the mismatching bases are gaps, 0 kcal/mol is used as well. First mismatches in hairpin loops are specified in msa.tstackh.dg. When the terminal base pair is non-canonical or contains gaps, 0 kcal/mol is used. If one of the mismatched bases is a gap, then the motif is treated as a dangling end and corresponding free energy from rna.dangle.dg is used. If both nucleotides in the mismatch are gaps, then 0 kcal/mol is used. For the sequence-independent coaxial stack of two adjacent helices with an intervening mismatch (msa.coaxstack.dg), 0 kcal/mol is used when non-canonical base pair is followed by a mismatch or one or more bases in the terminal base pair or the mismatching bases are gaps. In these cases, the two helices would be beyond the single mismatch allowed for coaxial stacking. Internal loops of size 2×2 are specified in msa.int22.dg. When one or more gaps appear in closing base pairs of 2×2 internal loops, then they are treated as “wildcards” (i.e. they are treated as though they could represent any of the canonical nucleotides A, C, G, or U) and the mean free energy change across all sequence possibilities for the gap is used. When one or more non-canonical base pairs close a 2×2 internal loop, we treat this as a larger internal loop. We assume a canonical base pair adjacent to the non-canonical pair. For a single non-canonical pair closure, the free energy is the mean of canonical closures at that position with additional penalty of 0.2 kcal/mol, which would be incurred by a 2×2 internal loop becoming a 3×3 internal loop. When both the closing base pairs for the 2×2 internal loop are non-canonical, the free energy change is the mean of canonical closures with twice the 0.2 kcal/mol penalty. When gaps and non-canonical base pairing appear in closing base pairs for 2×2 loops, we treat gaps as “wildcards” and assume presence of a canonical base pair adjacent to the non-canonical pair. We use the mean free energy increment of all relevant possible sequences and canonical closures, with additional penalty of 0.2 kcal/mol. For all these cases, if the free energy change calculated is less than 0 kcal/mol, we use a minimum free energy change of 0 kcal/mol instead. Internal loops of size 2×1 are specified in msa.int21.dg. When one or more gaps appear in closing base pairs of 2×1 internal loops, they are treated as “wildcards” and the mean free energy change across all sequence possibilities for the “wildcard” is used. When a non-canonical base pair closes 2×1 internal loop, we treat this as a larger internal loop. We assume a canonical base pair adjacent to the non-canonical pair. The free energy change is the mean of canonical closures at that position, with an additional penalty of 0.2 kcal/mol, which would be incurred by 2×1 loop becoming 3×2 loop. When both the closing base pairs for 2×1 internal loop are non-canonical, we make similar assumptions. The free energy change is the mean of canonical closures with twice the penalty. When gaps and non-canonical base pairing appear in closing base pairs for 2×1 loops, we treat gaps as “wildcards” and assume presence of a canonical base pair adjacent to the non-canonical pair. We take the mean of relevant possibilities and canonical closures, with an additional penalty of 0.2 kcal/mol. For all these cases, if free energy change calculated is less than 0 kcal/mol, we use 0 kcal/mol instead. Internal loops of size 1×1 are parameterized in msa.int11.dg. When one or more gaps appear in closing base pairs of 1×1 internal loops, hey are treated as “wildcards” and the mean free energy across all sequence possibilities is used. When a non-canonical base pair closes a 1×1 internal loop, we treat this as larger internal loop. We assume a canonical base pair adjacent to the non-canonical pair. The free energy is the mean of canonical closures at that position, with additional penalty of 0.2 kcal/mol (which would be incurred by 1×1 loop becoming 2×2 loop). When both the closing base pairs for 1×1 internal loop are non-canonical, we make similar assumptions. The free energy change is the mean of canonical closures with twice the penalty. When gaps and non-canonical base pairing appear in closing base pairs for 1×1 loops, we treat gaps as “wildcards” and assume the presence of a canonical base pair adjacent to the non-canonical pair. We take mean of relevant possibilities and canonical closures, with an additional penalty of 0.2 kcal/mol. For all these cases, if free energy change calculated is less than 0 kcal/mol, we use 0 kcal/mol instead. The specific hairpin loop sequence tables, including the tetraloop table (msa.tloop.dg), the triloop table (msa.triloop.dg), and the hexaloop table (msa.hexaloop.dg) remain identical to their tables in the 2004 nearest neighbor parameters. msa.miscloop.dg remains the same as rna.miscloop.dg. This includes the multibranch loop parameters and the helix end penalties. The length-dependent free energy penalties for closing hairpin, internal, and bulge loops (contained in msa.loop.dg) are sequence independent and remain the same as those in 2004 nearest neighbor parameters. ACKNOWLEDGEMENTS This work was supported by United States National Institutes of Health grant NIH R35 GM145283 to D.H.M. We thank Liang Huang and Apoorv Malik (Oregon State University) for discussions. Funder Information Declared National Institutes of Health, https://ror.org/01cwqze88 , R35GM145283 Footnotes https://rna.urmc.rochester.edu/RNAstructure.html References ↵ Bellaousov S , Mathews DH . 2010 . ProbKnot: Fast prediction of RNA secondary structure including pseudoknots . RNA 16 : 1870 – 1880 . doi: 10.1261/rna.2125310 . OpenUrl Abstract / FREE Full Text ↵ Bernhart SH , Hofacker IL , Will S , Gruber AR , Stadler PF . 2008 . RNAalifold: Improved consensus structure prediction for RNA alignments . BMC Bioinformatics 9 : 1 – 13 . doi: 10.1186/1471-2105-9-474 . OpenUrl CrossRef PubMed ↵ Brown JW . 1999 . The Ribonuclease P Database . Nucleic Acids Res 27 : 314 – 314 . doi: 10.1093/NAR/27.1.314 OpenUrl CrossRef PubMed Web of Science ↵ Cannone JJ , Subramanian S , Schnare MN , Collett JR , D’Souza LM , Du Y , Feng B , Lin N , Madabusi LV. , Müller KM , et al. 2002 . The comparative RNA web (CRW) site: an online database of comparative sequence and structure information for ribosomal, intron, and other RNAs . BMC Bioinformatics 3 . doi: 10.1186/1471-2105-3-2 . OpenUrl CrossRef PubMed ↵ Cech TR , Steitz JA . 2014 . The Noncoding RNA Revolution—Trashing Old Rules to Forge New Ones . Cell 157 : 77 – 94 . doi: 10.1016/J.CELL.2014.03.008 . OpenUrl CrossRef PubMed Web of Science ↵ Damase TR , Sukhovershin R , Boada C , Taraballi F , Pettigrew RI , Cooke JP . 2021 . The Limitless Future of RNA Therapeutics . Front Bioeng Biotechnol 9 : 1 – 24 . doi: 10.3389/fbioe.2021.628137 . OpenUrl CrossRef ↵ Flynn RA , Pedram K , Malaker SA , Batista PJ , Smith BAH , Johnson AG , George BM , Majzoub K , Villalta PW , Carette JE , et al. 2021 . Small RNAs are modified with N-glycans and displayed on the surface of living cells . Cell 184 : 3109 - 3124.e22 . doi: 10.1016/j.cell.2021.04.023 . OpenUrl CrossRef ↵ Fu Y , Sharma G , Mathews DH . 2014 . Dynalign II: common secondary structure prediction for RNA homologs with domain insertions . Nucleic Acids Res 42 : 13939 – 13948 . doi: 10.1093/NAR/GKU1172 . OpenUrl CrossRef PubMed ↵ Fu Y , Xu ZZ , Lu ZJ , Zhao S , Mathews DH . 2015 . Discovery of Novel ncRNA Sequences in Multiple Genome Alignments on the Basis of Conserved and Stable Secondary Structures . PLoS One 10 : e0130200 . doi: 10.1371/journal.pone.0130200 . OpenUrl CrossRef PubMed ↵ Gao W , Jones TA , Rivas E. 2021 . Discovery of 17 conserved structural RNAs in fungi . Nucleic Acids Res 49 : 6128 – 6143 . doi: 10.1093/nar/gkab355 . OpenUrl CrossRef PubMed ↵ Griffiths-Jones S , Moxon S , Marshall M , Khanna A , Eddy SR , Bateman A. 2005 . Rfam: annotating non-coding RNAs in complete genomes . Nucleic Acids Res 33 : D121 – D124 . doi: 10.1093/nar/gki081 . OpenUrl CrossRef PubMed Web of Science ↵ Gruber AR , Findeiß S , Washietl S , Hofacker IL , Stadler PF . 2010 . RNAz 2.0: improved noncoding RNA detection . Pac Symp Biocomput 69–79 . doi: 10.1142/9789814295291_0009 . OpenUrl CrossRef ↵ Gutell RR . 2014 . Ten lessons with Carl Woese about RNA and comparative analysis . RNA Biol 11 : 254 – 272 . doi: 10.4161/rna.28718 . OpenUrl CrossRef PubMed ↵ Gutell RR , Lee JC , Cannone JJ . 2002 . The accuracy of ribosomal RNA comparative structure models . Curr Opin Struct Biol 12 : 301 – 310 . doi: 10.1016/S0959-440X(02)00339-1 . OpenUrl CrossRef PubMed Web of Science ↵ Hagey RJ , Elazar M , Pham EA , Tian S , Ben-Avi L , Bernardin-Souibgui C , Yee MF , Moreira FR , Rabinovitch MV , Meganck RM , et al. 2022 . Programmable antivirals targeting critical conserved viral RNA secondary structures from influenza A virus and SARS-CoV-2 . Nat Med . doi: 10.1038/s41591-022-01908-x . OpenUrl CrossRef PubMed ↵ Hofacker IL , Fekete M , Stadler PF . 2002 . Secondary structure prediction for aligned RNA sequences . J Mol Biol 319 : 1059 – 1066 . doi: 10.1016/S0022-2836(02)00308-X . OpenUrl CrossRef PubMed Web of Science ↵ Jühling F , Mörl M , Hartmann RK , Sprinzl M , Stadler PF , Pütz J. 2009 . tRNAdb 2009: compilation of tRNA sequences and tRNA genes . Nucleic Acids Res 37 . doi: 10.1093/NAR/GKN772 . OpenUrl CrossRef ↵ Katoh K , Toh H. 2008 . Improved accuracy of multiple ncRNA alignment by incorporating structural information into a MAFFT-based framework . BMC Bioinformatics 9 : 1 – 13 . doi: 10.1186/1471-2105-9-212 . OpenUrl CrossRef PubMed ↵ Kierzek E , Zhang X , Watson RM , Kennedy SD , Szabat M , Kierzek R , Mathews DH . 2022 . Secondary structure prediction for RNA sequences including N6-methyladenosine . Nat Commun 2022 131 13 : 1 – 10 . doi: 10.1038/s41467-022-28817-4 . OpenUrl CrossRef PubMed ↵ Kim J , Walter AE , Turner DH . 1996 . Thermodynamics of coaxially stacked helixes with GA and CC mismatches . Biochemistry 35 : 13753 – 13761 . doi: 10.1021/bi960913z . OpenUrl CrossRef PubMed ↵ Knudsen B , Hein J. 2003 . Pfold: RNA secondary structure prediction using stochastic context-free grammars . Nucleic Acids Res 31 : 3423 . doi: 10.1093/NAR/GKG614 . OpenUrl CrossRef PubMed Web of Science ↵ Knudsen B , Hein J. 1999 . RNA secondary structure prediction using stochastic context-free grammars and evolutionary history . Bioinformatics 15 : 446 – 454 . doi: 10.1093/bioinformatics/15.6.446 . OpenUrl CrossRef PubMed Web of Science ↵ Li S , Zhang H , Zhang L , Liu K , Liu B , Mathews DH , Huang L. 2021 . LinearTurboFold: Linear-time global prediction of conserved structures for RNA homologs with applications to SARS-CoV-2 . Proc Natl Acad Sci U S A 118 . doi: 10.1073/pnas.2116269118 . OpenUrl Abstract / FREE Full Text ↵ Liu B , Mathews DH , Turner DH . 2010 . RNA pseudoknots: folding and finding . F1000 Biol Rep 2 . doi: 10.3410/B2-8 . OpenUrl CrossRef ↵ Lu ZJ , Gloor JW , Mathews DH . 2009 . Improved RNA secondary structure prediction by maximizing expected pair accuracy . RNA 15 : 1805 – 1813 . doi: 10.1261/rna.1643609 . OpenUrl Abstract / FREE Full Text ↵ Malik A , Zhang L , Gautam M , Dai N , Li S , Zhang H , Mathews DH , Huang L. 2024 . LinearAlifold: Linear-time consensus structure prediction for RNA alignments . J Mol Biol 168694. doi: 10.1016/j.jmb.2024.168694 . OpenUrl CrossRef ↵ Mathews DH . 2019 . How to benchmark RNA secondary structure prediction accuracy . Methods 162–163 : 60 – 67 . doi: 10.1016/J.YMETH.2019.04.003 . OpenUrl CrossRef ↵ Mathews DH . 2006 . Revolutions in RNA Secondary Structure Prediction . J Mol Biol 359 : 526 – 532 . doi: 10.1016/J.JMB.2006.01.067 . OpenUrl CrossRef PubMed Web of Science ↵ Mathews DH . 2004 . Using an RNA secondary structure partition function to determine confidence in base pairs predicted by free energy minimization . RNA 10 : 1178 . doi: 10.1261/RNA.7650904 . OpenUrl Abstract / FREE Full Text ↵ Mathews DH , Disney MD , Childs JL , Schroeder SJ , Zuker M , Turner DH . 2004 . Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure . Proc Natl Acad Sci U S A 101 : 7287 – 7292 . doi: 10.1073/PNAS.0401799101 . OpenUrl Abstract / FREE Full Text ↵ Mathews DH , Moss WN , Turner DH . 2010 . Folding and Finding RNA Secondary Structure . Cold Spring Harb Perspect Biol 2 . doi: 10.1101/cshperspect.a003665 . OpenUrl Abstract / FREE Full Text ↵ Mathews DH , Sabina J , Zuker M , Turner DH . 1999 . Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure . J Mol Biol 288 : 911 – 940 . doi: 10.1006/JMBI.1999.2700 . OpenUrl CrossRef PubMed Web of Science ↵ McCaskill JS . 1990 . The equilibrium partition function and base pair binding probabilities for RNA secondary structure . Biopolymers 29 : 1105 – 1119 . doi: 10.1002/bip.360290621 . OpenUrl CrossRef PubMed Web of Science ↵ Mittal A , Turner DH , Mathews DH . 2024 . NNDB: An Expanded Database of Nearest Neighbor Parameters for Predicting Stability of Nucleic Acid Secondary Structures . J Mol Biol . doi: 10.1016/J.JMB.2024.168549 . OpenUrl CrossRef ↵ Nawrocki EP , Burge SW , Bateman A , Daub J , Eberhardt RY , Eddy SR , Floden EW , Gardner PP , Jones TA , Tate J , et al. 2015 . Rfam 12.0: updates to the RNA families database . Nucleic Acids Res 43 : D130 – D137 . doi: 10.1093/NAR/GKU1063 . OpenUrl CrossRef PubMed ↵ Reuter JS , Mathews DH . 2010 . RNAstructure: Software for RNA secondary structure prediction and analysis . BMC Bioinformatics 11 : 1 – 9 . doi: 10.1186/1471-2105-11-129 . OpenUrl CrossRef PubMed ↵ Rivas E. 2020 . RNA structure prediction using positive and negative evolutionary information . PLOS Comput Biol 16 : e1008387 . doi: 10.1371/JOURNAL.PCBI.1008387 . OpenUrl CrossRef PubMed ↵ Rivas E , Lang R , Eddy SR . 2012 . A range of complex probabilistic models for RNA secondary structure prediction that includes the nearest-neighbor model and more . RNA 18 : 193 – 212 . doi: 10.1261/RNA.030049.111 . OpenUrl Abstract / FREE Full Text ↵ Rosenblad MA , Gorodkin J , Knudsen B , Zweib C , Samuelsson T. 2003 . SRPDB: Signal Recognition Particle Database . Nucleic Acids Res 31 : 363 – 364 . doi: 10.1093/NAR/GKG107 . OpenUrl CrossRef PubMed Web of Science ↵ Seeman SE , Gorodkin J , Backofen R. 2008 . Unifying evolutionary and thermodynamic information for RNA folding of multiple alignments . Nucleic Acids Res 36 : 6355 – 6362 . doi: 10.1093/NAR/GKN544 . OpenUrl CrossRef PubMed Web of Science ↵ Seetin MG , Mathews DH . 2012 . RNA structure prediction: An overview of methods . Methods Mol Biol 905 : 99 – 122 . doi: 10.1007/978-1-61779-949-5_8 . OpenUrl CrossRef PubMed ↵ Sievers F , Wilm A , Dineen D , Gibson TJ , Karplus K , Li W , Lopez R , McWilliam H , Remmert M , Söding J , et al. 2011 . Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega . Mol Syst Biol 7 : 539 . doi: 10.1038/msb.2011.75 . OpenUrl CrossRef PubMed Web of Science ↵ Spitale RC , Incarnato D. 2022 . Probing the dynamic RNA structurome and its functions . Nat Rev Genet 2022 243 24 : 178 – 196 . doi: 10.1038/s41576-022-00546-w . OpenUrl CrossRef PubMed ↵ Szikszai M , Wise M , Datta A , Ward M , Mathews DH . 2022 . Deep learning models for RNA secondary structure prediction (probably) do not generalize across families . Bioinformatics 38 : 3892 – 3899 . doi: 10.1093/BIOINFORMATICS/BTAC415 . OpenUrl CrossRef PubMed ↵ Szutkowska B , Wieczorek K , Kierzek R , Zmora P , Peterson JM , Moss WN , Mathews DH , Kierzek E. 2022 . Secondary Structure of Influenza A Virus Genomic Segment 8 RNA Folded in a Cellular Environment . Int J Mol Sci 23 : 2452 . doi: 10.3390/IJMS23052452 . OpenUrl CrossRef PubMed ↵ Szymanski M , Barciszewska MZ , Erdmann VA , Barciszewski J. 2002 . 5S Ribosomal RNA Database . Nucleic Acids Res 30 : 176 – 178 . doi: 10.1093/NAR/30.1.176 . OpenUrl CrossRef PubMed Web of Science ↵ Tan Z , Fu Y , Sharma G , Mathews DH . 2017 . TurboFold II: RNA structural alignment and secondary structure prediction informed by multiple homologs . Nucleic Acids Res 45 : 11570 – 11581 . doi: 10.1093/NAR/GKX815 . OpenUrl CrossRef PubMed ↵ Tyagi R , Mathews DH . 2007 . Predicting helical coaxial stacking in RNA multibranch loops . RNA 13 : 939 – 951 . doi: 10.1261/RNA.305307 . OpenUrl Abstract / FREE Full Text ↵ Virtanen P , Gommers R , Oliphant TE , Haberland M , Reddy T , Cournapeau D , Burovski E , Peterson P , Weckesser W , Bright J , et al. 2020 . SciPy 1.0: fundamental algorithms for scientific computing in Python . Nat Methods 2020 173 17 : 261 – 272 . doi: 10.1038/s41592-019-0686-2 . OpenUrl CrossRef PubMed ↵ Walter AE , Turner DH . 1994 . Sequence Dependence of Stability for Coaxial Stacking of RNA Helixes with Watson-Crick Base Paired Interfaces . Biochemistry 33 : 12715 – 12719 . doi: 10.1021/bi00208a024 . OpenUrl CrossRef PubMed ↵ Walter AE , Turner DH , Kim J , Lyttle MH , Müller P , Mathews DH , Zuker M. 1994 . Coaxial stacking of helixes enhances binding of oligoribonucleotides and improves predictions of RNA folding . Proc Natl Acad Sci 91 : 9218 – 9222 . doi: 10.1073/PNAS.91.20.9218 . OpenUrl Abstract / FREE Full Text ↵ Ward M , Sun H , Datta A , Wise M , Mathews DH . 2019 . Determining parameters for non-linear models of multi-loop free energy change . Bioinformatics 35 : 4298 – 4306 . doi: 10.1093/bioinformatics/btz222 . OpenUrl CrossRef PubMed ↵ Woodson SA , Crothers DM . 1987 . Proton nuclear magnetic resonance studies on bulge-containing DNA oligonucleotides from a mutational hot-spot sequence . Biochemistry 26 : 904 – 912 . doi: 10.1021/BI00377A035 . OpenUrl CrossRef PubMed Web of Science ↵ Xia T , SantaLucia J , Burkard ME , Kierzek R , Schroeder SJ , Jiao X , Cox C , Turner DH . 1998 . Thermodynamic parameters for an expanded nearest-neighbor model for formation of RNA duplexes with Watson -Crick base pairs . Biochemistry 37 : 14719 – 14735 . doi: 10.1021/bi9809425 . OpenUrl CrossRef PubMed ↵ Zhang C , Zhang Y , Pyle AM . 2023 . rMSA: A Sequence Search and Alignment Algorithm to Improve RNA Structure Modeling . J Mol Biol 435 : 167904 . doi: 10.1016/J.JMB.2022.167904 . OpenUrl CrossRef PubMed ↵ Zhang L , Zhang H , Mathews DH , Huang L. 2019 . ThreshKnot: Thresholded ProbKnot for Improved RNA Secondary Structure Prediction . arXiv Biomol . doi: 10.48550/arXiv.1912.12796 . OpenUrl CrossRef ↵ Zhou Y , Lu C , Wu QJ , Wang Y , Sun ZT , Deng JC , Zhang Y. 2008 . GISSD: Group I Intron Sequence and Structure Database . Nucleic Acids Res 36 . doi: 10.1093/NAR/GKM766 . OpenUrl CrossRef ↵ Znosko BM , Silvestri SB , Volkman H , Boswell B , Serra MJ . 2002 . Thermodynamic parameters for an expanded nearest-neighbor model for the formation of RNA duplexes with single nucleotide bulges . Biochemistry 41 : 10406 – 10417 . doi: 10.1021/BI025781Q . OpenUrl CrossRef PubMed Web of Science ↵ Zuber J , Schroeder SJ , Sun H , Turner DH , Mathews DH . 2022 . Nearest neighbor rules for RNA helix folding thermodynamics: improved end effects . Nucleic Acids Res 50 : 5251 – 5262 . doi: 10.1093/NAR/GKAC261 . OpenUrl CrossRef PubMed ↵ Zwieb C , Gorodkin J , Knudsen B , Burks J , Wower J. 2003 . tmRDB (tmRNA database) . Nucleic Acids Res 31 : 446 – 447 . doi: 10.1093/NAR/GKG019 . OpenUrl CrossRef PubMed Web of Science View the discussion thread. Back to top Previous Next Posted July 28, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following AlignmentFold and AlignmentPartition: Improving the align-then-fold approach for RNA secondary structure prediction Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share AlignmentFold and AlignmentPartition: Improving the align-then-fold approach for RNA secondary structure prediction Abhinav Mittal , David H. Mathews bioRxiv 2025.07.23.666478; doi: https://doi.org/10.1101/2025.07.23.666478 Share This Article: Copy Citation Tools AlignmentFold and AlignmentPartition: Improving the align-then-fold approach for RNA secondary structure prediction Abhinav Mittal , David H. Mathews bioRxiv 2025.07.23.666478; doi: https://doi.org/10.1101/2025.07.23.666478 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13895) Bioinformatics (41953) Biophysics (21456) Cancer Biology (18595) Cell Biology (25521) Clinical Trials (138) Developmental Biology (13381) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24323) Genetics (15612) Genomics (22511) Immunology (17738) Microbiology (40401) Molecular Biology (17184) Neuroscience (88623) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7644) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)