← Về thư mục
📄 / / profiles / r-and-d / raw_evidence / OpenAI4S / skills / bioskills / bio-rna-structure-secondary-structure-prediction / SKILL.md

name: bio-rna-structure-secondary-structure-prediction description: Predicts RNA secondary structure with ViennaRNA, treating the Boltzmann ensemble (partition function, base-pair probabilities, centroid, MEA, stochastic samples) as the object rather than a single MFE fold. Covers consensus folding from alignments (RNAalifold), SHAPE-constrained folding, RNA-RNA interaction (RNAcofold/RNAduplex/RNAup), local and linear-time methods for long RNA, and pseudoknot-aware tools. Use when folding an RNA and choosing between MFE, centroid, MEA, or ensemble sampling; judging whether a single structure is well-defined; folding long RNAs where a global MFE is meaningless; handling suspected pseudoknots; or weighing thermodynamic versus comparative versus deep-learning prediction. origin: openai4s category: bioskills/rna-structure metadata: tool_type: cli primary_tool: ViennaRNA third_party: name: GPTomics/bioSkills repository: https://github.com/GPTomics/bioSkills commit: d91ed3d563019e649dc854c56ccd62551359488a license: MIT

Version Compatibility

Reference examples tested with: ViennaRNA 2.6+, matplotlib 3.8+, numpy 1.26+

Before using code patterns, verify installed versions match. If versions differ: - CLI: <tool> --version then <tool> --help to confirm flags - Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Secondary Structure Prediction

"Predict the secondary structure of my RNA sequence" -> Compute base pairs under a nearest-neighbor thermodynamic model, but report the Boltzmann ENSEMBLE (partition function, base-pair probabilities, centroid/MEA, per-base confidence), not a single fold. - CLI: RNAfold -p for single-sequence ensemble folding - CLI: RNAalifold for consensus structure from an alignment - CLI: RNAcofold / RNAduplex / RNAup for RNA-RNA interaction - Python: import RNA (fold_compound) for scripted ensemble analysis

The governing principle: the MFE is one sample from an ensemble, not "the structure"

The minimum free energy (MFE) structure is the single lowest-energy fold, but every possible structure has probability proportional to exp(-G/RT). The partition function (McCaskill 1990) sums over that whole Boltzmann ensemble and yields the probability of each base pair, not just one fold. The MFE is frequently NOT the biologically relevant structure: riboswitches sit at a poised switch between two folds, many RNAs are genuinely dynamic, and the functional fold is often slightly suboptimal. Report the MFE alone and the uncertainty is hidden.

Three consequences shape every decision below: - Trust the ENSEMBLE, not the point estimate. Use partition-function quantities (base-pair probability > 0.9, low positional entropy, low ensemble diversity) as the per-pair and per-base confidence. A low MFE with a diffuse base-pair-probability matrix means the single structure is not trustworthy. - Accuracy is modest and length-dependent. Single-sequence thermodynamic folding recovers ~70-73% of base pairs for short, well-behaved RNAs (tRNA, 5S rRNA) and degrades sharply beyond ~700 nt. A single global MFE of a multi-kilobase mRNA or viral genome is close to meaningless -> use local (RNAplfold), linear-time (LinearFold/LinearPartition), comparative, or probing-restrained methods. - ViennaRNA folds only nested canonical + wobble (G-U) pairs. Pseudoknots are silently excluded (general pseudoknot prediction is NP-hard); non-canonical pairs (sheared G-A, base triples, the Leontis-Westhof geometric families) are invisible. A confidently wrong nested fold is the failure mode for frameshift elements, riboswitch aptamers, and viral UTRs.

Which ViennaRNA program for which question

Question / input Tool Why
Fold one sequence, want structure + confidence RNAfold -p MFE + partition function (centroid/MEA/diversity)
Ensemble free energy only, no dot plot (speed) RNAfold -p0 skips base-pair probabilities, ~50% faster
Consensus structure from an alignment of homologs RNAalifold thermodynamics + covariation
Two RNAs that dimerize (full model + concentrations) RNAcofold (-c) intra+inter pairs, equilibrium species
Fast two-strand hybridization screen RNAduplex inter-molecular pairs only, no internal structure
sRNA/miRNA-target where site accessibility matters RNAup opening energy + hybridization (correct for buried sites)
Sample alternative conformations RNAsubopt -p N / fc.pbacktrack(n) Boltzmann sampling
Long mRNA: local pairing / target accessibility RNAplfold windowed pair + unpaired probabilities
Long RNA/genome: find local structured elements RNALfold locally stable structures, bounded span
Sequence longer than a few kb LinearFold / LinearPartition O(n), avoids a meaningless global O(n^3) fold
Pseudoknot suspected IPknot / ProbKnot / Knotty nested folders structurally cannot
Have SHAPE/DMS reactivities RNAfold --shape (Deigan) restrain folding with experimental data

Which "answer" to report

Goal Answer Note
Quick single estimate MFE over-calls weak pairs; not for long/low-complexity RNA
Conservative, ensemble-representative structure centroid minimum total base-pair distance to ensemble; fewer false pairs, can under-pair
Best single structure (esp. with probing data) MEA (tune gamma) high gamma -> more pairs/recall, low gamma -> fewer/precision
Per-pair / per-base confidence base-pair probability (>0.9) + positional entropy structure-agnostic confidence track
Conformational switching / multiple states stochastic sampling cluster the samples into populations
Is one structure well-defined? ensemble diversity (low) + ensemble defect length-relative, not an absolute cutoff

RNAfold: single-sequence ensemble folding

# MFE only
echo "GGGCUAUUAGCUCAGUUGGUUAGAGCGCACCCCUGAUAAGGGUGAGGUCGCUGAUUCGAAUUCAGCAUAGCCCA" | RNAfold --noPS

# Ensemble: partition function + base-pair probabilities + centroid + MEA + ensemble diversity
# --noPS suppresses the *_dp.ps / *_ss.ps PostScript files RNAfold writes to the CWD by default
echo ">myRNA" > rna.fa && echo "GGGCUAUUAGCUCAGUUGGUUAGAGCGCACC" >> rna.fa
RNAfold -p --MEA --noLP --noPS < rna.fa

Key flags (verified against the current manpage):

Option Effect
-p partition function + base-pair-probability matrix (unlocks centroid/MEA/diversity)
-p0 ensemble free energy ONLY, no base-pair probabilities (faster)
--MEA[=gamma] maximum-expected-accuracy structure (default gamma 1.0); --MEA implies -p
-d2 dangling-end model (default); use -d0 for comparative/alignment folding to avoid dangle artifacts
-d3 also allow coaxial stacking of adjacent helices in multiloops (MFE folding only; the partition function -p ignores -d3 and falls back to -d2, so ensemble quantities do not reflect it)
--noLP forbid lonely (isolated) base pairs; standard for well-folded RNA
--maxBPspan N cap base-pair span; crude knob for long sequences
-T 37 folding temperature in Celsius (default 37)
--shape FILE / --shapeMethod SHAPE-directed folding (see structure-probing)
-g allow G-quadruplex formation (default off; turn on for G-rich sequences where G4s compete with canonical pairing)
--noPS suppress PostScript drawings (always set in scripts to avoid CWD clutter)

Centroid, MEA, sampling, and per-base confidence (Python)

fc.pf() MUST be called before bpp(), centroid(), MEA(), pbacktrack(), positional_entropy(), ensemble_defect(), or mean_bp_distance() -- they all read the partition-function matrices and silently return empty/garbage otherwise.

import RNA

seq = 'GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA'
fc = RNA.fold_compound(seq)

mfe_struct, mfe = fc.mfe()
_, ensemble_g = fc.pf()                 # partition function; ensemble G <= MFE always

centroid, _ = fc.centroid()             # conservative, fewer false-positive pairs
mea_struct, _ = fc.MEA()                # MEA(gamma); gamma>1 favors pairing (recall), gamma<1 precision
diversity = fc.mean_bp_distance()       # ensemble diversity: low = well-defined (read RELATIVE to length)
defect = fc.ensemble_defect(mfe_struct) # expected wrongly-paired positions of this structure vs ensemble
entropy = fc.positional_entropy()       # per-base Shannon entropy: low = confidently paired-or-unpaired

Sample alternative conformations from the Boltzmann ensemble (riboswitches, bistable RNA) with the CLI RNAsubopt -p N (N stochastic samples) or the Python fc.pbacktrack(N) after fc.pf(); cluster the samples to find conformational populations. fc.pbacktrack requires a ViennaRNA build with stochastic backtracking enabled -- if it returns nothing, use RNAsubopt -p N.

Decision rule: one well-defined structure -> centroid or MEA; report stability/confidence -> partition-function quantities (base-pair probabilities + positional entropy); conformational switching -> stochastic sampling; quick single estimate -> MFE (with the caveats above).

Constrained and SHAPE-directed folding

import RNA

seq = 'GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA'

# Hard constraints: force positions unpaired or paired
fc = RNA.fold_compound(seq)
fc.hc_add_up(35, RNA.CONSTRAINT_CONTEXT_ALL_LOOPS)   # 1-indexed; force position 35 unpaired
fc.hc_add_bp(1, 72, RNA.CONSTRAINT_CONTEXT_ALL_LOOPS)
constrained, c_mfe = fc.mfe()

# Soft SHAPE pseudo-energy restraint (Deigan model). The reactivity vector is 1-INDEXED:
# prepend -999 as a placeholder for index 0; -999 elsewhere means "no data" (NOT zero reactivity).
# m=1.8, b=-0.6 is the standard SHAPE pair (Hajdin 2013, the ViennaRNA default), not Deigan 2009's own m=2.6/b=-0.8.
fc2 = RNA.fold_compound(seq)
reactivities = [-999.0] + [0.1, 0.05, 0.8, 0.9] + [-999.0] * (len(seq) - 4)
fc2.sc_add_SHAPE_deigan(reactivities, 1.8, -0.6)
shape_struct, shape_mfe = fc2.mfe()   # this energy INCLUDES the SHAPE pseudo-energy; do NOT compare it to the unrestrained MFE

The energy returned after a SHAPE restraint folds in the pseudo-energy bonus, so it is not on the same scale as an unconstrained MFE -- compare the STRUCTURES (base-pair distance, SHAPE agreement), not the two energy numbers. See structure-probing for obtaining reactivities and for the SHAPE-vs-DMS parameter choice.

Comparative (consensus) folding: homologs beat a single sequence

Evolution conserves structure while sequence drifts, so a compensatory substitution (an A-U in one species becoming G-C at the same two columns) is direct evidence of a real pair that thermodynamics alone cannot see. RNAalifold folds an alignment with a combined thermodynamic + covariation score.

# Consensus structure; format (Stockholm/Clustal/FASTA) is auto-detected, alignment is positional.
# --ribosum_scoring improves covariation detection; -d0 avoids dangle artifacts at gapped columns.
RNAalifold --ribosum_scoring -d0 -p --noPS alignment.sto
Option Effect
--cfactor covariation weight (default 1.0; lower leans on thermodynamics)
--nfactor penalty for sequences that cannot form the consensus pair (default 1.0)
--ribosum_scoring use RIBOSUM covariation matrices (recommended)
-p consensus partition function + base-pair probabilities

RNAalifold accuracy depends entirely on alignment quality and real covariation: near-identical sequences carry no covariation signal and it degrades toward noisy single-sequence folding. RNAalifold assumes a FIXED, correct alignment; when homologs cannot be aligned reliably, TurboFold II co-estimates alignment AND structure across the sequences jointly (Tan et al. 2017) and is the better choice. A predicted consensus structure is a HYPOTHESIS until covariation is statistically validated -- test it with R-scape (see covariation-analysis), which found no significant covariation support for the proposed HOTAIR/Xist/SRA lncRNA structures.

Long RNA: do not fold one global structure

For a multi-kilobase mRNA, lncRNA, or viral genome, a single global O(n^3) MFE both over-pairs across long ranges and is slow, and folding is co-transcriptional and local in reality.

# Windowed local pairing + per-position UNPAIRED (accessibility) probabilities
RNAplfold -W 200 -L 150 -u 30 < long_rna.fa      # -W window, -L max base-pair span, -u accessibility region length

# Scan for locally stable structured elements with a bounded span
RNALfold -L 150 < long_rna.fa

# Linear-time approximate MFE (LinearFold) and partition function (LinearPartition), if installed
echo "GGGAAACCC..." | linearfold
echo "GGGAAACCC..." | linearpartition

LinearFold's 5'->3' beam search can match or improve accuracy versus the exact cubic algorithm on long RNAs (the exact global model is not more correct when a single structure is not meaningful), besides being far faster.

Pseudoknots: ViennaRNA cannot, by construction

The standard dynamic programming forbids crossing pairs, and general pseudoknot prediction is NP-hard (Lyngso & Pedersen 2000) -- RNAfold/RNAalifold silently return the best NESTED structure. Suspect a pseudoknot for tmRNA, telomerase RNA, RNase P, many riboswitch aptamers (SAM-II, preQ1), -1 ribosomal frameshift elements, IRES, and group I/II intron cores.

Tool Class / method Note
IPknot integer programming over base-pair probabilities fast, broad class, the pragmatic default
ProbKnot MEA assembly from McCaskill probabilities (RNAstructure) any topology, fastest/most scalable
Knotty MFE over the broad CCJ class more complex crossing topologies
pknotsRG MFE over restricted simple recursive pseudoknots, O(n^4) narrower class

Pseudoknot prediction is substantially less accurate and more expensive than nested folding -- treat any predicted pseudoknot as a hypothesis to corroborate with a second tool, covariation, or probing.

RNA-RNA interaction: pick by the binding question

Tool Models Use when
RNAcofold both intramolecular AND intermolecular pairs; -c gives equilibrium concentrations full dimerization model
RNAduplex inter-molecular pairs only (no internal structure), fast first-pass target screen
RNAup opening (accessibility) energy + hybridization energy sRNA/miRNA-target where the site may be buried in structure (the physically correct choice)

The two strands are concatenated with & (RNAcofold/RNAduplex); RNAup takes the two sequences on separate lines.

# Full dimer model (intra + inter pairs); -p for the heterodimer partition function
echo "GCGCGCAUAU&AUAUGCGCGC" | RNAcofold -p --noPS
# With -c, RNAcofold reads the two monomer concentrations and reports equilibrium fractions of the
# five species (AB, AA, BB, A, B) -- use it to ask how much dimer actually forms, not just whether it is favorable.

# Fast inter-molecular-only hybridization screen (no internal structure)
echo "GCGCGCAUAU&AUAUGCGCGC" | RNAduplex

# Accessibility-corrected sRNA/miRNA-target binding: opening energy + hybridization (-b includes both)
RNAup -b < two_sequences.fa

Thermodynamics vs deep learning: DL is not a default for novel RNA

Deep-learning predictors (SPOT-RNA, UFold, E2Efold) report high accuracy ON FAMILIES SEEN IN TRAINING, but under family-fold cross-validation that removes train/test homology their accuracy collapses to at or below the thermodynamic baseline (Szikszai et al. 2022); the apparent gains are intra-family memorization, and benchmark sets are ~55% rRNA / >90% rRNA+tRNA (Flamm et al. 2022). For a genuinely novel RNA (unseen Rfam family), no single-sequence method (DL or thermodynamic) is reliable -- the robust evidence is covariation (R-scape) and experimental probing. If using DL, prefer the thermodynamics-integrated hybrid MXfold2 over end-to-end nets; never cite intra-family accuracy as proof of de-novo performance.

"Is this more structured than random?" (z-score vs shuffled controls, RNAz, randfold): shuffles MUST preserve DINUCLEOTIDE composition (Altschul-Erikson), because MFE is dominated by GC content and base stacking, a dinucleotide property -- a mononucleotide shuffle inflates significance and makes almost anything look stable. For an alignment, RNAz combines a dinucleotide-controlled z-score with a structure conservation index (SCI = consensus MFE / mean single-sequence MFE; ~1 = a conserved structure), but reads Clustal/MAF, not Stockholm. A negative z-score means "more stable than random," NOT "this structure is correct"; covariation is the stronger evidence standard.

Method-class selection (the big fork)

Situation Recommended Avoid as default
Single novel RNA, no homologs thermodynamic ensemble (RNAfold -p / RNAstructure) pure end-to-end DL
Aligned homologs with covariation RNAalifold + R-scape validation single-sequence MFE
Homologs but no trusted alignment TurboFold II (joint alignment + structure) align-then-RNAalifold on a poor alignment
Have SHAPE/DMS reactivities probing-restrained folding (Deigan/Zarringhalam) unrestrained MFE
Long mRNA / transcriptome scale RNAplfold / LinearFold+LinearPartition global O(n^3) MFE
Pseudoknot biology IPknot/ProbKnot/Knotty + cross-check RNAfold (cannot)
Willing to use DL MXfold2 (thermodynamics-integrated) E2Efold/UFold on unseen families

When competing methods or parameters are in play, verify current behavior against the installed tool's --help and the latest docs before trusting a number.

Getting the structure out and drawing it

Dot-bracket is the default text form; CT and BPSEQ are the interchange formats downstream tools (ProbKnot, IPknot, RNAstructure) read and write. To DRAW a structure, use forna (web), R2DT (template-based standard layouts for known families), or VARNA; RNAfold's own *_ss.ps PostScript drawing is exactly what --noPS suppresses, so drop --noPS when the built-in diagram is wanted. The example renders the base-pair probability dot plot with matplotlib.

Common Errors

Symptom Cause Fix
AttributeError: module 'RNA' has no attribute 'sequence_shuffle' no such ViennaRNA function use a dinucleotide-preserving shuffle (ushuffle, esl-shuffle -d) for z-scores
bpp()/centroid()/MEA() return empty or garbage fc.pf() not called first call fc.pf() before any ensemble quantity
*_dp.ps / *_ss.ps files appearing in the working directory RNAfold/RNAalifold write PostScript by default pass --noPS (and run in a scratch dir)
Long mRNA gives one improbable global fold global MFE is meaningless past ~700 nt use RNAplfold / LinearFold / LinearPartition
Predicted structure has a pseudoknot the tool "missed" RNAfold cannot represent crossing pairs use IPknot / ProbKnot / Knotty
Consensus structure looks confident but is wrong RNAalifold trusts the alignment; no real covariation validate with R-scape; check alignment quality
RNAalifold --aln alignment.sto treated as input flag --aln is an OUTPUT (annotated PostScript) flag pass the alignment positionally: RNAalifold alignment.sto
SHAPE-constrained fold barely changes reactivity vector mis-indexed or zeros where data is missing vector is 1-indexed (prepend -999); use -999 for no-data, not 0

References