Hierarchical analysis of RNA secondary structures with pseudoknots based on sections

TL;DR

A new method for predicting RNA pseudoknots by analyzing local sections improves computational efficiency and reveals insights into RNA folding dynamics.

Contribution

A hierarchical approach for pseudoknot prediction that reduces computational cost and reveals RNA folding follows sequential co-transcriptional dynamics.

Findings

Over 90% of pseudoknots occur in the top 3% of section pairs ranked by minimum free energy gain.

The method achieves sensitivity >0.9 and positive predictive value >0.8 for 2-section pseudoknots.

3-section pseudoknots show asymmetric behavior with early-formed connections being predictable but later ones not.

Abstract

Predicting RNA structures containing pseudoknots remains computationally challenging due to high processing costs and complexity. While standard methods for pseudoknot prediction require O(N6) time complexity, we present a hierarchical approach that significantly reduces computational cost while maintaining prediction accuracy. Our method analyzes RNA structures by dividing them into contiguous regions of unpaired bases (“sections”) derived from known secondary structures. We examine pseudoknot interactions between sections using a nearest-neighbor energy model with dynamic programming. Our algorithm scales as O(n2ℓ4), offering substantial computational advantages over existing global prediction methods. Analysis of 726 transfer messenger RNA and 454 Ribonuclease P RNA sequences reveals that biologically relevant pseudoknots are highly concentrated among section pairs with large minimum free energy (MFE) gain. Over 90% of connected section pairs appear within just the top 3% of section pairs ranked by MFE gain. For 2-clusters, our method achieves high prediction accuracy with sensitivity exceeding 0.9 and positive predictive value above 0.8. For 3-clusters, we discovered asymmetric behavior where “former” section pairs (formed early in the sequence) are predicted accurately, while “latter” section pairs do not follow local energy predictions. This asymmetry suggests that complex pseudoknot formation follows sequential co-transcriptional folding rather than global energy minimization, providing insights into RNA folding dynamics.

RNA molecules fold into structures to perform biological functions. However, predicting complex RNA structures known as “pseudoknots” is computationally expensive. Current methods often attempt to calculate the entire structure simultaneously, which requires significant computational resources. In this paper, we introduce a hierarchical approach that simplifies pseudoknot prediction. We break the RNA sequence into smaller “sections” of unpaired bases and calculate the energy required for these sections to bind locally, rather than solving for the global structure. Our analysis shows that strong local interactions are favored by biology; with over 90% of pseudoknots occurring within the top 3% of the most energetically favorable section pairs. This finding allows us to focus computational effort on the small subset of interactions that are most likely to form pseudoknots, rather than testing every possible combination. Our method achieves >90% sensitivity for simple 2-section pseudoknots. However, for complex 3-section pseudoknots, only early-forming connections are predictable. This reveals that RNA does not simply fold into the most stable structure. Instead, folding is sequential, with earlier regions establishing interactions that constrain the final structure before synthesis of the later regions.

Figures

Graphical representations of RNA structures and section-based pseudoknot classifications.Yellow lines: secondary structure base pairs; red lines: pseudoknot base pairs; blue regions: sections (unpaired nucleotides); black regions: non-sections. a) RNA structure in stem-loop and linear representations. b) Stem-loop (left) and linear (right) representations of the same structure. Blue arrows indicate sections (contiguous unpaired regions). The black arrow indicates a non-section where nucleotides participate in secondary structure pairing. c) Topologically distinct pseudoknot structures of both 2-clusters and 3-clusters.

Energy model calculation examples demonstrating local versus global structural considerations.a) Pseudoknot base pair structure. The yellow bands represent sets of base pairs belonging to a secondary structure, while the red band represents a set of consecutive pseudoknots. Bottom right shows detailed base pair structure in the pseudoknot. b) Pseudoknot connecting sections across multi-loops. The newly formed loop that includes the pseudoknot contains all sections from the connected multi-loops, requiring global rather than local energy calculations.

Implementation flowchart of the RNA pseudoknot structure prediction algorithm.The flowchart illustrates the hierarchical approach of our method, starting with input parsing, followed by section identification from the secondary structure. Energy calculations determine the minimum free energy (MFE) between section pairs. Structure prediction then branches into 2-cluster and 3-cluster analysis paths, with the latter using weighted energy minimization (w=0.8) to balance former and latter section pair contributions. The output includes predicted structures, base pair assignments, and accuracy metrics.

Distribution of MFE gain for section pairs in RNA databases.a) Histogram showing exponential decay of MFE gain frequency for all section pairs in tmRNA and RNase P RNA sequences from RNAstrand database. b) Simple probabilistic model explaining the observed exponential decay, assuming independent base pairing probability and ∼1.0 kcal/mol energy contribution per base pair.

Analysis of pseudoknot formation likelihood based on energetic predictions.a) Connecting probability as a function of MFE gain for tmRNA and RNase P RNA sequences, defined as the ratio of actually connected section pairs to all section pairs with a given MFE value. b) Representative example of a non-connecting section pair with substantial MFE gain (80 kcal/mol), illustrating how conformational entropy losses can outweigh local energy gains in long sections (∼104 and ∼48 nucleotides).

Optimization of weighting parameter for 3-cluster pseudoknot structure prediction.Prediction accuracy metrics (sensitivity and positive predictive value) as a function of weight parameter w applied to the free energy contribution of latter section pairs in the weighted energy minimization approach. Optimal performance occurs at w = 0.8, indicating preferential weighting of former section pair interactions while maintaining contribution from latter pairs in type-1 3-cluster configurations.

Equations

Full text

Introduction

Ribonucleic acid (RNA) plays fundamental roles across all life forms, carrying genetic information, regulating gene expression, and performing various catalytic functions. Non-coding RNA, which does not translate into proteins, plays important roles in many catalytic and regulatory cellular processes [1–4]. The biological functions of non-coding RNA have strong connections with its molecular structure, particularly pseudoknots [5–8]. Thus, predicting and understanding the structure of RNA for any given base sequence is important in biology, medicine, pharmacy, and other related fields. A fundamental challenge in RNA structure prediction is defining which interactions constitute the secondary structure versus pseudoknots. RNA secondary structure has multiple definitions in literature [9–11]. In this work, we define RNA secondary structure as nested base pairs that can be represented without crossing lines, while pseudoknots are non-planar interactions that create crossing patterns [12–16]. While no universal criterion exists for distinguishing secondary structure from pseudoknot base pairs [17,18], our definition enables computational tractability, aligns with standard structural notations (e.g. dot-bracket format), and provides a framework for our hiearchical method of identifying unpaired sections and predicting pseudoknot interactions between them. RNA secondary structure forms the basic frame of an RNA structure, typically exhibiting greater thermodynamic stability than pseudoknot interactions. Several reliable algorithms have been developed to predict RNA secondary structure for any given base sequence, such as dynamic programming algorithm (DPA) [19–24] and stochastic methods [25–27]. DPA is one of the most widely used methods and has been applied in many different areas [28,29]. DPA accurately predicts RNA secondary structures by minimizing the free energy contribution under a local energy model in O(L^3^) time for an RNA sequence of length L.

Predicting RNA structures that include pseudoknots remains challenging due to computational complexity and high processing costs. The most straightforward approach would be to DPA to find minimum free energy structures across all possible base pairings. However, this approach is NP-complete even with simple energy models that only consider individual base pair contributions [30–34]. In response, algorithms have been developed that focuses on specific subsets of pseudoknotted structures [35–40]. These methods improve efficiency by eliminating physically impossible configurations early in the analysis. For example, a dynamic programming algorithm achieving O(N^6^) time complexity and O(N^4^) space complexity was introduced [10], providing the first method to fold optimal pseudoknotted RNAs using standard thermodynamic models. Though widely used, this approach becomes computationally intensive for longer RNA sequences. A later study showed that simpler dynamic programming techniques could achieve O(N^5^) time complexity [41], offering improved performance but lacking the ability to handle complex pseudoknot configurations like kissing-hairpins, which are commonly found in RNA structures. Alternative approaches have emerged, including stochastic methods [42] and novel computational strategies [43,44]. While these techniques show promise in specific scenarios, comparative analyses [45] reveal accuracy limitations, particularly when dealing with complex pseudoknot structures [46].

In order to overcome the trade-off between the accuracy of the model and the computational cost, there is now an increasing focus on hierarchical folding method [40,47,48]. The hierarchical folding method uses a two-step approach: first predicting the pseudoknot-free secondary structure, then calculating pseudoknot base pairs based on this foundation. A proposed DPA variant [49] finds pseudoknotted structures in O(N^3^) time and O(N^2^) space by applying this hierarchical principle, matching the time efficiency of algorithms that only predict secondary structures. This approach can handle complex configurations like kissing hairpins, which traditionally required O(N^6^) time to process. However, comparative studies show these methods face significant challenges with accuracy, particularly for complex pseudoknots in longer sequences [36,45,50]. These limitations appear largely due to the fact that pseudoknot structures are often long-range and thus affected by the many possible 3D RNA conformations, which are difficult to account for in current models.

In this paper, we analyze local pseudoknot interactions by examining how different “sections” of an RNA (derived from known secondary structures) are connected via pseudoknots. Using a nearest-neighbor energy model and dynamic programming, we calculate the minimum free energy of potential pseudoknot formations between unpaired regions. Our approach focuses strictly on local energy contributions between sections and excludes interactions between different pseudoknots.

Building on this section-wise interaction framework, we develop a computationally efficient method that scales as $[eqn]$ , where n represents the number of sections and $[eqn]$ represents their typical length. This hierarchical approach provides further understanding into how local structural elements contribute to global RNA configurations.

Methods

Definitions of terms

We define the terminologies used in the research of RNA structure.

Primary structure of an RNA refer to its chemical sequence.Secondary structure of an RNA is the set of all local nested pairing of complementary bases that can be represented in a planar graph without crossing lines. Even pairings between distant bases (such as initial and final bases) are considered part of the secondary structure if they maintain this planar property. Common secondary structure motif includes hairpin loops, hairpin stems, helical duplexes, bulges, internal loops, and multi-loops. Internal loops consist of unpaired nucleotides connecting exactly two helical segments $[eqn]$ , whereas multibranch loops are junction points where three or more helices converge $[eqn]$ .Pseudoknots are base pairings between secondary motifs that cannot be represented without crossing lines in a planar path. These crossing interactions are tetiary structures. Pseudoknots are defined by these crossing patterns within the global RNA topology, not as isolated elements. The non-planarity emerges from the relationship between these interactions and the secondary structure, reflecting the topological complexity characteristic of pseudoknots.Pseudoknot base pair is a base pair that participates in a non-nested pairing arrangement, creating crossing interactions when represented in a planar diagram. These pairs are part of the global topological structure where bases i,j form a pair and bases k,l form another pair such that i < k < j < l, resulting in crossing lines in a 2D representation.Pseudoknotted secondary structure is defined as the set of all base pairs, including pseudoknots. It is specifically used when the set of base pairs includes pseudoknots.

The definitions above are explained graphically in Fig 1A. Next, we will introduce terminologies to be used in our analysis.

A section is a contiguous subset of bases in the RNA sequence that are not involved in base pairing, as determined from the known structures in our dataset. Specifically, if the ith base and jth base (i < j) are not paired in the secondary structure, they belong to the same section if all bases between them (from position $[eqn]$ to j – 1) are also not involved in base pairing. Each section represents a single-stranded region that may potentially form pseudoknot interactions with other sections. This definition relies on having accurate secondary structure information, which in our study comes from established RNA structure databases.Section pair refers to a pair of sections. The ith section and jth section are connected if there exists a pseudoknot that connects a base in ith section and another in jth section. See Fig 1B for a graphical description of sections.N-cluster is a cluster consisting of N sections connected by pseudoknots. The number of topologically distinct pseudoknot structures increases rapidly once $[eqn]$ . We will only consider 2-clusters and 3-clusters in this work. The different types of 3-clusters are shown in Fig 1C.Pseudoknot order is a metric for RNA structural complexity defined by the minimum number of base pair decompositions to create a nested structure [9]. In our work, 2-clusters correspond to pseudoknot order 1 (1 decomposition; e.g. H/K-type) while 3 clusters correspond to pseudoknot order 2 (2 decompositions; L/M-type).

Graphical representations of RNA structures and section-based pseudoknot classifications.Yellow lines: secondary structure base pairs; red lines: pseudoknot base pairs; blue regions: sections (unpaired nucleotides); black regions: non-sections. a) RNA structure in stem-loop and linear representations. b) Stem-loop (left) and linear (right) representations of the same structure. Blue arrows indicate sections (contiguous unpaired regions). The black arrow indicates a non-section where nucleotides participate in secondary structure pairing. c) Topologically distinct pseudoknot structures of both 2-clusters and 3-clusters.

Energy model

We use the nearest-neighbor energy model which is based on the chemical and thermodynamic properties of the RNA chains. Nearest-neighbor energy model is widely used to predict secondary structure in conjunction with DPA [51–58]. In the nearest-neighbor energy model, values of free energy contribution are assigned to each closed loop (e.g. stackings, bulge loops, interior loops, multi-loops) in RNA secondary structure, and the total free energy of the structure is calculated as the sum of all the free energy contribution from the loops [59]. The parameters and detailed rules in the energy model used in our analysis are based on mfold version 3.6 [20,21]. While these parameters were derived from experiments in salt-free conditions, we apply them as an approximation for pseudoknot prediction in the absence of comprehensive energy tables calibrated for varying salt concentrations.

Formally, the energy model is defined as

$$[eqn]$$

Using Eq 1, we calculate the free energy contribution of the pseudoknot between two sections in the example shown (in linear representation) in Fig 2. According to the conventional rule of the nearest-neighbor energy model, the free energy contribution of the pseudoknot in Fig 2A is given by:

$$[eqn]$$

Energy model calculation examples demonstrating local versus global structural considerations.a) Pseudoknot base pair structure. The yellow bands represent sets of base pairs belonging to a secondary structure, while the red band represents a set of consecutive pseudoknots. Bottom right shows detailed base pair structure in the pseudoknot. b) Pseudoknot connecting sections across multi-loops. The newly formed loop that includes the pseudoknot contains all sections from the connected multi-loops, requiring global rather than local energy calculations.

However, in this work, the contribution of loops that include base pairs in secondary structures are neglected, and only loops that include pseudoknot base pairs are considered. Hence, the free energy contribution of the psuedoknot in Fig 2 is calculated as

$$[eqn]$$

Our energy model uses two main nearest-neighbour thermodynamic parameters: stacking energies and loop energies. Stacking energies represent free energy contributions of stacked pairs formed by consecutive base pairs, with values ranging from –3.4 kcal/mol (most stable) to +1.5 kcal/mol (least stable). Loop energies compose of multiple parameter arrays such as destabilising values for interior loops (4.1 to 7.4 kcal/mol for lengths 1–30), bulge loops energies (3.9 to 6.7 kcal/mol), mismatch energy parameters (–2.7 to +0.2 kcal/mol) [21]. We exclude free energy contributions from loops containing secondary structure base pairs to focus on local interactions between sections under analysis. This exclusion is essential for our hierarchical approach. For example, when sections within our target pair reside in a multi-loop (Fig 2B), calculating their local pseudoknot structure would necessitate predicting pseudoknot formations across other section pairs simultaneously, forcing expansion from local to global structure prediction and defeating the purpose of our section-based approach. Moreover, global base pair structures predicted by DPA are known to be inaccurate even with additional energy parameters for pseudoknot motifs [49], likely due to inability to account for 3D conformation effects [60,61]. Rather than incorporating 3D conformational effects into DPA, we focus on predicting local structure using only local information.

Our energy model further simplifies entropic contributions by treating loop contributions as purely additive. This additive approximation, independent analysis of section pairs, and loop exclusion achieve $[eqn]$ computational cost while maintaining focus on local interactions. The empirical validity of these approximations is demonstrated in our Results section. For exact thermodynamic parameter values, refer to our code. For exact thermodynamic parameter values used in this study, refer to our Github.

While our approach for energy calculations shares mathematical similarities with RNA hybridization algorithms through nearest-neighbor dynamic programming models, there are key differences in our implementation. Our method specifically addresses intramolecular interactions between sections within a single RNA molecule, unlike hybridization algorithms that focus on intermolecular binding [62,63]. Our approach also incorporates topological constraints specific to pseudoknot formation within an existing RNA secondary structure. For 3-clusters analysis, we apply weighted energy minimization (as shown in Fig 6, where w = 0.8 proved optimal), a consideration not typically present in standard hybridization algorithms [64–66]. Furthermore, our method distinguishes between different pseudoknot topologies (Fig 1C) and prevents crossing pseudoknot interactions within section pairs. We apply these techniques to analyse local pseudoknot interactions within known unpaired regions, which helps demonstrate how RNA structures can be examined in a hierarchical manner.

Algorithm

Our algorithm determines the MFE structure between section pairs by considering all possible structures where base pairs can form pseudoknots with the secondary structure, constrained so that pseudoknot base pairs within the same section pair cannot cross each other. This allows pseudoknots between different sections while preventing nested pseudoknots within the same section pair. Moreover, pairing between two sections from the same loop is prohibited because such base pairs would be nested rather than pseudoknotted. Sections are identified by analyzing the loop regions from the predicted secondary structure, where each continuous stretch of unpaired nucleotides within a loop forms a distinct section. The algorithm systematically identifies all sections by scanning the RNA sequence and grouping consecutive unpaired bases that belong to the same loop structure. Fig 3 shows the flowchart of our proposed algorithm.

Implementation flowchart of the RNA pseudoknot structure prediction algorithm.The flowchart illustrates the hierarchical approach of our method, starting with input parsing, followed by section identification from the secondary structure. Energy calculations determine the minimum free energy (MFE) between section pairs. Structure prediction then branches into 2-cluster and 3-cluster analysis paths, with the latter using weighted energy minimization (w=0.8) to balance former and latter section pair contributions. The output includes predicted structures, base pair assignments, and accuracy metrics.

We denote the two sections in the given section pair as “section1” and “section2”, with length $[eqn]$ and $[eqn]$ respectively. We define a matrix C whose entry C(i,j) ( $[eqn]$ ) represents the minimum value of free energy contribution over candidate structures in which ith base in section1 and jth base in section2 are paired and $[eqn]$ th base in section1 or $[eqn]$ th base in section2 are not connected when either $[eqn]$ or $[eqn]$ . We set $[eqn]$ when the ith base in section1 and the jth base in section2 cannot form a base pair. By definition, the MFE between the section-pairs that we are calculating is the minimum entry of the C matrix. Thus, the problem of finding the MFE turns into the calculation of matrix C. The entries of matrix C can be calculated using the following recursion relation.

$$[eqn]$$

After calculating the C matrix (see Algorithm 1), we determine the MFE structure between section pairs by selecting the binding sites through dynamic programming that minimizes free energy, where the algorithm identifies the minimum entry in the C matrix, representing the optimal base pairing configuration. Base pairs are formed only between complementary bases (AU, UA, CG, GC, GU, UG). This restriction to canonical pairs (Watson-Crick and wobble) is due to the availability of thermodynamic parameters in mfold 3.6. Analysis of our datasets shows canonical pairs constitute over 99% of pseudoknot pairs in transfer messenger RNA (tmRNA) and over 90% in ribonuclease P RNA (RNase P RNA). For sections to interact, they must be from different loops, contain complementary bases, and have negative free energy contribution due to their interaction. For 2-clusters, we use the MFE structure between the section pair. For 3-clusters, we apply weighted energy minimization: (energy of former pair) + w × (energy of latter pair), where $[eqn]$ . Sections are ordered by their position in the RNA sequence (5’ to 3’ direction), with the former section pair connecting sections appearing earlier and the latter section pair connecting sections appearing later. This ordering reflects co-transcriptional folding, where 5’ regions can form interactions before 3’ regions are fully synthesized [67,68]. We ensure consistency by preventing bases already paired in former sections from forming pairs in latter sections. This hierarchical approach constructs the global structure by first identifying favorable local interactions, then combining them while maintaining structural constraints.

Algorithm 1 Dynamic programming algorithm for MFE calculation between section pairs.

$[eqn]$ , $[eqn]$ : section lengths;  s1, s2: starting positions;  C: DP matrix;  $[eqn]$ : possible base pairs;  mismatches: mismatch types between bases

1: function SEGMENT_MFE(segments, loops, bases, $[eqn]$ , $[eqn]$ , seg1, seg2)

2:   $[eqn]$ get_section_lengths(segments, $[eqn]$ , seg1, seg2)

3:   $[eqn]$ get_section_starting_points(segments, $[eqn]$ , seg1, seg2)

4:   if $[eqn]$ = $[eqn]$ then return $[eqn]$

5:   end if          $[eqn]$   Sections in same loop cannot form pseudoknots

6:   Initialize $[eqn]$ and $[eqn]$

7:   for $[eqn]$ to l1−1 do

8:    for $[eqn]$ to l2−1 do

9:     $[eqn]$ pair_type( $[eqn]$ , $[eqn]$ )

10:     $[eqn]$ mismatch_type( $[eqn]$ , $[eqn]$ )

11:    end for

12:   end for

13:   Initialize $[eqn]$ and set $[eqn]$

14:   for $[eqn]$ downto 0 do

15:    for $[eqn]$ to l2−1 do

16:     if $[eqn]$ then

17:      $[eqn]$           $[eqn]$   i,j cannot pair

18:     else

19:      $[eqn]$           $[eqn]$   Case: i,j is the first bond

20:      if i<l1−1 and j > 0 then          $[eqn]$   Stacking

21:       $[eqn]$ $[eqn]$ , C[i + 1][j−1] +  $[eqn]$ , $[eqn]$

22:      end if

23:      if i<l1−1 then          $[eqn]$   Bulge loop case 1

24:       for $[eqn]$ to j do

25:        $[eqn]$ $[eqn]$ , C[i + 1][j−d] +  $[eqn]$ , $[eqn]$ , mismatches[i][j− d + 1], $[eqn]$

26:       end for

27:      end if

28:      if j > 0 then          $[eqn]$   Bulge loop case 2

29:       for $[eqn]$ to l1−i do

30:        $[eqn]$ ( $[eqn]$ , C[i + d][j−1] +  $[eqn]$ , $[eqn]$ , mismatches[i +  d−1][j], mismatches[i + 1][j−1]))

31:       end for

32:      end if

33:      for all $[eqn]$ do   $[eqn]$   Interior loop

34:       $[eqn]$ ( $[eqn]$ , $[eqn]$ $[eqn]$ , $[eqn]$ , $[eqn]$ , mismatches[i + 1][j−1]))

35:      end for

36:      $[eqn]$

37:      $[eqn]$

38:     end if

39:    end for

40:   end for

41:   return MFE

42: end function

We prove the computational complexity of our algorithm is O(n^2^l^4^). If n is the number of sections and l the typical length of one section, we compute the MFE for all possible section pairs for n sections, yielding $[eqn]$ . For each section pair with lengths l1 and l2, we construct matrix C of size $[eqn]$ . For each entry C(i,j), there are three cases: stacking with O(1) complexity, bulge loop with O(l1) and O(l2) complexity respectively, interior loop with $[eqn]$ complexity due to two nested loops of length l1 and l2. Thus, the worst-case time complexity for each C(i,j) entry is $[eqn]$ or O(l^4^) assuming $[eqn]$ . Thus, for all O(n^2^) section pairs, the total computational cost is O(n^2^l^4^). For typical RNA molecules in our dataset, we observed $[eqn]$ , where N is the total RNA length. This relationship suggests our algorithm’s complexity in terms of total RNA length is approximately O(N^2^).

Our algorithm was implemented in C and compiled with gcc optimization flags. Runtime analysis on our complete dataset (1,180 sequences, 102–1,331 nt) showed a total execution time of 7.54 seconds, with mean runtime of 6.39 ms per sequence (median: 5.35 ms, IQR: 3.98–6.65 ms). Only 6.8% of sequences were classified as statistical outliers beyond 1.5×IQR, demonstrating consistent performance suitable for large-scale RNA database analysis.

Results and discussion

In this work, we analyzed 726 tmRNA sequences and 454 RNase P RNA sequences from the RNAstrand database [69]. These RNA families were selected because they represent a comprehensive and diverse dataset of biologically relevant structures with well-documented pseudoknots. Rather than selecting a small subset of examples, we analyzed the complete set of these RNA families available in the database to ensure a thorough and unbiased evaluation of our approach. Both families consist of relatively long RNA sequences with many pseudoknots, making them ideal test cases for our section-based prediction method.

Basic properties of investigated RNA sequences

Table 1 summarises key characteristics of tmRNA and RNase P RNA sequences from the RNAstrand database. Analysis of the RNAstrand database reveals structural similarities between the two RNA classes (Table 1). Mean sequence lengths are 368 nt (tmRNA) and 333 nt (RNase P RNA), with ∼30 sections per sequence and mean section lengths of 8 nt and 5 nt respectively. tmRNA sequences show greater variability in their measurements. The increase in standard deviation is due to sequences exceeding 1000 bases, with each sequence containing a single extended section that forms no pseudoknot connections with other sections. These outlier sequences do not impact our analysis of pseudoknot interactions between section pairs. The substantial length of these RNA sequences presents a challenge, as conventional prediction methods struggle to accurately determine the structure of such extended sequences [36].

Table 1: Structural characteristics of RNA sequences analyzed from RNAstrand database.Section length refers to the number of unpaired nucleotides within each contiguous unpaired region.

Predicting the connected section pairs based on MFE

Fig 4A shows the distribution of MFE gain (absolute value of MFE) for all section pairs. Frequency decays exponentially with increasing MFE gain. The exponential decay follows from two assumptions. Firstly, the probability of forming each additional base pair between two sections is independent of existing structure (“probability of base pairing”). Secondly, each base pair contributes approximately 1.0 kcal/mol to the free energy, representing the average stacking energy in our model. Under these assumptions, the histogram in Fig 4 fits an exponential function with parameter a:

$$[eqn]$$

Distribution of MFE gain for section pairs in RNA databases.a) Histogram showing exponential decay of MFE gain frequency for all section pairs in tmRNA and RNase P RNA sequences from RNAstrand database. b) Simple probabilistic model explaining the observed exponential decay, assuming independent base pairing probability and ∼1.0 kcal/mol energy contribution per base pair.

The probability of base pairing can be calculated as

$$[eqn]$$

Curve fitting analysis reveals base pairing probabilities of 0.608 for tmRNA and 0.632 for RNase P RNA. These values reflect fundamental RNA base pairing constraints. With six possible pairing types (AU, UA, CG, GC, UG, GU), the probability of two adjacent bases forming a new base pair is 3/8 (0.375). The probability of at least one of two adjacent base pairs forming next to an existing pair is 1−(5/8)^2^ = 39/64 (0.609), matching our observed pairing probabilities.

Our results show that MFE structures between section pairs primarily consist of stacking interactions and small internal loops of length 2. The exponential decay in the distribution reflects the rarity of section pairs with large MFE gain. Specifically, only 5.5% of pseudoknots in tmRNA and 2.8% in RNase P RNA have MFE gain exceeding 5.0 kcal/mol. Even fewer exceed 10.0 kcal/mol: just 0.47% in tmRNA and 0.23% in RNase P RNA.

We investigated the connecting probability for both tmRNA and RNase P RNA in our database (Fig 5). Connecting probability is defined as the ratio between connected section pairs and all section pairs with a given MFE:

$$[eqn]$$

Analysis of pseudoknot formation likelihood based on energetic predictions.a) Connecting probability as a function of MFE gain for tmRNA and RNase P RNA sequences, defined as the ratio of actually connected section pairs to all section pairs with a given MFE value. b) Representative example of a non-connecting section pair with substantial MFE gain (80 kcal/mol), illustrating how conformational entropy losses can outweigh local energy gains in long sections (∼104 and ∼48 nucleotides).

Connecting probability correlates with high MFE gain as expected, since larger energy gains stabilize structures. However, we observed non-connecting section pairs with unexpectedly large MFE gains. These sections have large total MFE gain but small MFE gain per base pair compared to genuinely connected pairs. The loss of conformational entropy, which our nearest-neighbor energy model ignores, outweighs the stabilizing effect of each individual pair.

Connected section pairs concentrate among those with large MFE gains, which represent a small fraction of all section pairs. We sorted all section pairs by MFE gain and examined the cumulative distribution of connected pairs. We find that more than 90% of connected section pairs in tmRNA are included in the top 10^4^ section pairs (3%) with largest maximum free energy gains (out of all the $[eqn]$ possible section pairs). Over 90% of connected section pairs fall within the top 3% by MFE gain for tmRNA (10^4^ of $[eqn]$ pairs) and top 1% for RNase P RNA (1,500 of 1.8×10^5^ pairs). Therefore, it is sufficient to only consider a small fraction of section pairs with large MFE gain in order to predict RNA structures with pseudoknots. On average, we only need to consider ∼14 section pairs per tmRNA sequence (10,000/726) or ∼3 section pairs per RNase P RNA sequence (1,500/454) to account for 90% of pseudoknots.

Comparison between MFE structures and real structures for connected section pairs

Connected section pairs exclusively form either 2-clusters or 3-clusters, with no observed instances of larger cluster formations (Table 2). Structures with pseudoknot order 1(2-clusters: 2337 in tmRNA, 543 in RNase P RNA) substantially outnumber structures with pseudoknot order 2 (3-clusters: 110 in tmRNA, 1 in RNase P RNA). This distribution is consistent with hierarchical folding, where lower-order pseudoknots form before higher-order structures. We compared MFE predictions with actual free energies (real FE) for connected 2-cluster section pairs. While the proportion of 2-clusters where real FE = MFE varies considerably between RNA classes (68.7% for tmRNA versus 39.9% for RNase P RNA), >80% of all 2-clusters demonstrate free energy gain >0.8 × MFE gain.

Table 2: The result of analyses based on section pairs.

We compared MFE-predicted base pair arrangements with actual structures for 2-clusters. For each cluster, we calculated sensitivity (TPR) and positive predictive value (PPV) using true positives (TP), false positives (FP), and false negatives (FN). When multiple MFE structures existed, we averaged these metrics:

$$[eqn]$$

With sensitivity >0.9 and PPV >0.8 for both RNA classes, our algorithm outperforms global predictions methods, which typically achieves a sensitivity score of 0.5–0.7 and PPV <0.7 for complex pseudoknotted structures [45].

It is important to acknowledge that our analysis relies on knowledge of the true secondary structure from which we derive the unpaired sections. This differs from the more challenging problem of de novo pseudoknot prediction without prior structural information. While our results demonstrate that local energetic considerations can effectively predict pseudoknot formations between known unpaired regions, additional research would be needed to determine whether similar accuracy could be achieved when starting with predicted rather than true secondary structures.

3-clusters exhibit several topologically distinct configurations (Fig 1C), with 110 type-1 3-clusters in tmRNA and only one type-2 3-cluster in RNase P RNA. Given limited type-2 samples, we analyze only type-1 3-clusters, designating section1-section2 (see Fig 1C) connections as “former pairs” and section2-section3 as “latter pairs”.

We analysed the 110 type-1 3-clusters in tmRNA by comparing MFE calculations with actual structures for both pair types, treating them independently (similar to 2-cluster analysis). For former section pairs, approximately 50% have real structures identical to predicted MFE structures, while over 80% have free energy contributions lower than 0.8×MFE (Table 3). These results closely resemble our 2-cluster prediction accuracy. In contrast, latter section pairs show markedly different behavior, with only a small fraction having free energy contributions approximating their MFE. Former section pairs demonstrate prediction accuracy comparable to 2-clusters, while latter section pairs exhibit particularly low PPV values. Across both pair types, sensitivity consistently exceeds PPV, indicating our method tends to predict excess base pairs. This over-prediction is expected, as our approach treats these connections as independent 2-clusters despite their presence in more complex 3-cluster structures.

Table 3: The comparison of MFE and free energy contribution of real structure for connected section pairs in 3-clusters.

The failure suggests that former section pairs form independently, whereas latter section pair formations are strongly influenced by pre-existing base-pair connections. This asymmetry appears to reflect the fundamental biology of RNA synthesis, which proceeds directionally from 5’ to 3’ along the sequence. During transcription, former section pairs likely establish their interactions before latter sections are even fully synthesised, as supported by previous research on co-transcriptional folding dynamics [70,71]. The consistency of this pattern across our dataset suggests it represents a general feature of RNA folding rather than a family-specific characteristic. These findings challenge the assumption that RNA structures in cells form solely through global free energy minimization [72–75]. Instead, our results suggest that the dynamic and sequential nature of RNA synthesis/folding plays a crucial role in determining the final structural configuration, particularly for complex pseudoknot arrangements. This sequential dependency within pseudoknot order 2 structures demonstrates that even within a single higher-order pseudoknot, simpler interactions must establish before complex interactions can form.

To evaluate whether our method’s accuracy depends on sequence length, we analyzed all RNA sequences from our dataset. For each sequence, we compared predicted connections (negative MFE values) against actual connections in known structures. We calculated the true positive rate (TPR) and positive predictive value (PPV) for each sequence and examined their correlation with sequence length using Pearson correlation analysis.

For TPR versus sequence length, r = −0.026 $[eqn]$ , while for PPV versus sequence length, r = 0.046 $[eqn]$ . Both metrics show negligible correlation with sequence length $[eqn]$ with high p-values (p > 0.05), confirming no significant dependency. Under the null hypothesis of no correlation, these results support the conclusion that prediction accuracy remains consistent across the tested sequence length range.

The prediction accuracy independence is due to our hierarchical approach. Our algorithm calculates minimum free energy locally between section pairs rather than globally across the entire RNA structure. Consequently, prediction accuracy is determined by local section-pair interactions and section length $[eqn]$ rather than total RNA length L, consistent with our $[eqn]$ computational complexity. For a given section length, computational cost scales with the number of section pairs, not with total sequence length, ensuring that prediction accuracy remains stable regardless of overall RNA size.

The comparison between MFE structures and real structures for 3-clusters

Our previous analysis treated section pairs independently, ignoring that bases in section2 cannot simultaneously pair with both section1 and section3. We now predict complete 3-cluster structures while enforcing this section pair dependency with two prediction strategies:

Equal weighting: Minimize the combined free energy contribution from both former and latter section pairs, treating them with equal importance.Sequential determination: First identify the minimum free energy structure of the former section pair, then determine the latter section pair structure while ensuring no bases in section2 are reused.

These approaches represent special cases of a weighted method that minimizes:

$$[eqn]$$

where 0 < w < 1. We implemented an energy minimisation algorithm that evaluates all possible base pairings between sections (Fig 6), calculating MFE contributions for each former-latter pair combination using nearest-neighbor thermodynamic parameters from mfold 3.6. The algorithm considers stacking interactions, bulge loops, and interior loops while preventing bases already paired in the former pair from participating in latter pair interactions. We determined optimal w by testing values from 0 to 1, calculating weighted energies, predicted base pairs, and accuracy metrics (sensitivity and PPV) for each. Optimal performance occurs at w = 0.8, indicating preferential weighting of former pairs while maintaining latter pair contributions. However, 3-cluster latter pair prediction remains poor across all w values, suggesting their formation is influenced by global structure or 3D conformation beyond our local energy model’s scope.

Optimization of weighting parameter for 3-cluster pseudoknot structure prediction.Prediction accuracy metrics (sensitivity and positive predictive value) as a function of weight parameter w applied to the free energy contribution of latter section pairs in the weighted energy minimization approach. Optimal performance occurs at w = 0.8, indicating preferential weighting of former section pair interactions while maintaining contribution from latter pairs in type-1 3-cluster configurations.

A limitation of our approach is the use of energy parameters derived from salt-free experimental conditions (mfold 3.6) to predict structures that typically form in the presence of salt. This approximation was necessary due to the lack of comprehensive energy tables for RNA in various salt concent rations. Despite this limitation, our method demonstrates strong predictive power for 2-clusters, suggesting that local pseudoknot formation may be reasonably approximated by these parameters even without explicitly modeling ionic effects. Future work could significantly improve prediction accuracy by incorporating salt-dependent energy parameters as they become available.

Our method is successful despite using additive entropy approximation for pseudoknots. Although previous work has shown that entropy contributions in pseudoknotted structures are not strictly additive [76,77], our results for 2-clusters (sensitivity >0.90, PPV >0.80) suggest this simplification works well within our hierarchical section-based framework. The high accuracy indicates that for common pseudoknot configurations, local energy minimization with simplified entropy terms captures the essential physics governing structure formation. The method’s limitations with latter section pairs in 3-clusters likely indicate where more complex entropic considerations become necessary. Nonetheless, this simplification allows for structural prediction of longer RNA sequences while maintaining good accuracy for the most common pseudoknot configurations.

Conclusion

We have presented a section-based approach to RNA pseudoknot analysis that employs mfold nearest-neighbor energy model with DPA. This method offers both computational efficiency and prediction accuracy advantages over traditional approaches. Our algorithm scales as $[eqn]$ , where n represents the number of sections and $[eqn]$ represents the average section length. Assuming section length remains independent of total RNA sequence length, the computational complexity reduces to O(L^2^) for an RNA of length L. A key finding of our research is the highly concentrated distribution of biologically relevant pseudoknots among section pairs with large MFE gains. By focusing on just the top 3% of section pairs ranked by MFE gain, we can capture approximately 90% of all connected section pairs. This observation enables substantial reductions in the search space required for structure prediction.

Our approach predicts pseudoknot order 1 structures (2-clusters) with sensitivity above 0.90 and positive predictive value exceeding 0.80. These results significantly outperform conventional methods that attempt to predict global pseudoknotted structures in a single step. This indicates that once unpaired regions are identified, 2-cluster pseudoknot formations between them appear to be primarily determined by local information. This observation suggests potential value in hierarchical approaches to RNA structure prediction, though such approaches would need to address the challenge of accurately identifying relevant unpaired regions and resolving potential competition between pseudoknot and regular secondary structure formation. Our work provides evidence for the effectiveness of local energy calculations specifically for the subproblem of pseudoknot prediction between known unpaired regions. However, significant limitations exist for pseudoknot order 2 structures (3-clusters). Nearly all 3-clusters conform to the type-1 configuration (Fig 1C). While our approach can predict former section pairs in these 3-clusters with accuracy comparable to 2-clusters, it performs poorly when predicting latter section pairs. The inadequacy of our local energy model for latter section pairs suggests that these pseudoknots are significantly influenced by factors beyond local interactions, such as the dynamic process of RNA folding and the global 3D RNA conformation. These findings highlight the limitations of purely thermodynamic approaches to complex RNA structure prediction and suggest that accurate modeling of 3-clusters may require incorporating kinetic folding pathways and long-range structural constraints.