---

name: "viennarna-structure-prediction" description: "Predict RNA secondary structure, MFE folding, base-pair probabilities, RNA-RNA interactions via ViennaRNA Python bindings. Pipeline: sequence → MFE → partition function and pair-probability matrix → dot-bracket → duplex. Use for siRNA/sgRNA targeting, ribozyme design, RNA accessibility. Use RNAfold CLI for batch use without Python." license: "MIT"

ViennaRNA Structure Prediction

Overview

ViennaRNA is the gold-standard toolkit for RNA secondary structure prediction based on thermodynamic nearest-neighbor parameters. It predicts the minimum free energy (MFE) structure and dot-bracket notation for a given RNA sequence, computes the full partition function to obtain base pair probabilities, and models RNA-RNA interactions via co-folding and duplex prediction. The Python bindings (import RNA) expose the full ViennaRNA C library with sequence-level and fold-compound APIs. Command-line programs (RNAfold, RNAalifold, RNAduplex) are also available and demonstrated here.

When to Use

• Predicting the minimum free energy secondary structure of an RNA sequence (mRNA, lncRNA, miRNA precursor, aptamer)
• Computing base pair probability matrices to assess structural uncertainty and identify well-defined stem-loops
• Designing or evaluating siRNA accessibility by folding the target mRNA region and checking for double-stranded structure
• Assessing sgRNA targeting efficiency by predicting guide RNA secondary structure that may reduce on-target activity
• Modeling RNA-RNA interactions (co-folding or duplex prediction) for miRNA-target binding or antisense oligonucleotide design
• Calculating folding free energies for a set of sequences to compare thermodynamic stability
• Use mfold (web server) or RNAstructure instead when you need Mfold algorithm predictions specifically or need the Efold partition function; ViennaRNA uses the Turner 2004 nearest-neighbor parameters and is the standard for research-grade thermodynamic prediction

Prerequisites

• Python packages: ViennaRNA (Python bindings), matplotlib, numpy
• Data requirements: RNA sequences as strings (ACGU alphabet; T is auto-converted to U by ViennaRNA)
• Environment: Python 3.8+; conda installation strongly recommended (handles C library dependencies)

# Install via conda (recommended)
conda install -c conda-forge -c bioconda viennarna

# Verify installation
python -c "import RNA; print(RNA.__version__)"
# 2.6.4

# Install additional Python dependencies
pip install matplotlib numpy pandas

# Optional: verify CLI tools are available
RNAfold --version
# RNAfold 2.6.4

Quick Start

import RNA

# Predict MFE structure for an RNA sequence
sequence = "GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA"
structure, mfe = RNA.fold(sequence)

print(f"Sequence:  {sequence}")
print(f"Structure: {structure}")
print(f"MFE:       {mfe:.2f} kcal/mol")
# Sequence:  GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA
# Structure: (((((((..((((........)))).(((((.......))))).....(((((.......))))))))))))....
# MFE:       -31.30 kcal/mol

Workflow

Step 1: Sequence Preparation and MFE Folding

Load an RNA sequence and compute its minimum free energy secondary structure using RNA.fold(). Validate the input and inspect the dot-bracket output.

import RNA

def prepare_sequence(seq: str) -> str:
    """Normalize sequence: uppercase, replace T→U, validate alphabet."""
    seq = seq.upper().replace("T", "U").strip()
    invalid = set(seq) - set("ACGUNX")
    if invalid:
        raise ValueError(f"Invalid characters in sequence: {invalid}")
    return seq

# E. coli tRNA-Phe (GenBank: M10217)
raw_seq = "GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA"
sequence = prepare_sequence(raw_seq)

structure, mfe = RNA.fold(sequence)

print(f"Sequence length: {len(sequence)} nt")
print(f"Structure:       {structure}")
print(f"MFE:             {mfe:.2f} kcal/mol")

# Validate: structure length must equal sequence length
assert len(structure) == len(sequence), "Structure and sequence length mismatch"

# Count stems (paired bases)
n_paired   = structure.count("(") + structure.count(")")
n_unpaired = structure.count(".")
print(f"Paired bases: {n_paired}  |  Unpaired bases: {n_unpaired}")
print(f"Stem fraction: {n_paired/len(sequence):.2f}")

Step 2: Create a Fold Compound for Advanced Analysis

The RNA.fold_compound object is the central API for partition function, base pair probabilities, and constrained folding.

import RNA

sequence = "GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA"

# Create fold compound (wraps the sequence with model parameters)
fc = RNA.fold_compound(sequence)

# Compute MFE structure via the fold compound API
structure, mfe = fc.mfe()
print(f"MFE structure: {structure}")
print(f"MFE:           {mfe:.2f} kcal/mol")

# Evaluate free energy of an alternative structure
alt_structure = "." * len(sequence)   # fully unfolded
energy = fc.eval_structure(alt_structure)
print(f"Fully unfolded energy: {energy:.2f} kcal/mol")
print(f"Folding stabilization:  {energy - mfe:.2f} kcal/mol")

Step 3: Partition Function and Base Pair Probabilities

Compute the thermodynamic partition function to obtain ensemble-level base pair probabilities. High-probability pairs indicate well-defined structural elements.

import RNA
import numpy as np

sequence = "GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA"
n = len(sequence)

fc = RNA.fold_compound(sequence)

# Step 1: MFE folding (required before pf for proper initialization)
structure_mfe, mfe = fc.mfe()

# Step 2: Rescale Boltzmann factors for numerical stability (optional but recommended)
fc.exp_params_rescale(mfe)

# Step 3: Compute partition function
structure_pf, gibbs_free_energy = fc.pf()
print(f"Gibbs free energy (ensemble): {gibbs_free_energy:.2f} kcal/mol")
print(f"MFE structure:  {structure_mfe}")
print(f"Centroid (pf):  {structure_pf}")

# Step 4: Retrieve base pair probability matrix
bpp = fc.bpp()   # returns (n+1)x(n+1) matrix; 1-indexed

# Convert to 0-indexed numpy array for analysis
probs = np.zeros((n, n))
for i in range(1, n + 1):
    for j in range(i + 1, n + 1):
        if bpp[i][j] > 0.0:
            probs[i - 1][j - 1] = bpp[i][j]
            probs[j - 1][i - 1] = bpp[i][j]

# Identify high-confidence pairs (p > 0.9)
high_conf = [(i, j, probs[i, j]) for i in range(n) for j in range(i + 1, n) if probs[i, j] > 0.9]
print(f"\nHigh-confidence base pairs (p > 0.9): {len(high_conf)}")
for i, j, p in high_conf[:5]:
    print(f"  {sequence[i]}{i+1} — {sequence[j]}{j+1}: p={p:.3f}")

Step 4: Visualize Base Pair Probability Matrix

Plot the base pair probability matrix as a heatmap to visualize structural regions.

import RNA
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

sequence = "GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA"
n = len(sequence)

fc = RNA.fold_compound(sequence)
structure_mfe, mfe = fc.mfe()
fc.exp_params_rescale(mfe)
fc.pf()
bpp = fc.bpp()

# Build matrix
probs = np.zeros((n, n))
for i in range(1, n + 1):
    for j in range(i + 1, n + 1):
        if bpp[i][j] > 0.0:
            probs[i - 1][j - 1] = bpp[i][j]
            probs[j - 1][i - 1] = bpp[i][j]

fig, ax = plt.subplots(figsize=(8, 7))
im = ax.imshow(probs, cmap="hot_r", vmin=0, vmax=1, origin="upper", aspect="equal")
plt.colorbar(im, ax=ax, label="Base pair probability")

ax.set_xlabel("Nucleotide position")
ax.set_ylabel("Nucleotide position")
ax.set_title(f"Base Pair Probability Matrix\n(n={n} nt, MFE={mfe:.2f} kcal/mol)")
plt.tight_layout()
plt.savefig("bpp_matrix.png", dpi=150, bbox_inches="tight")
print("Saved: bpp_matrix.png")

Step 5: RNA-RNA Duplex Prediction (Co-folding)

Use RNA.cofold() to predict the interaction between two RNA sequences by concatenating them with an & separator.

import RNA

# miRNA (hsa-miR-21-5p) and its target sequence in mRNA (PTEN 3'UTR region)
mirna_seq  = "UAGCUUAUCAGACUGAUGUUGA"
target_seq = "UCAACAUCAGUCUGAUAAGCUA"   # approximate complementary target

# Co-fold: concatenate with & separator
cofold_seq = mirna_seq + "&" + target_seq
structure, mfe = RNA.cofold(cofold_seq)

print(f"miRNA:          {mirna_seq}")
print(f"Target:         {target_seq}")
print(f"Co-fold MFE:    {mfe:.2f} kcal/mol")

# Parse the structure — & is retained in output
n1, n2 = len(mirna_seq), len(target_seq)
struct_mirna  = structure[:n1]
struct_ampersand = structure[n1]
struct_target = structure[n1 + 1:]

print(f"miRNA structure:  {struct_mirna}")
print(f"Target structure: {struct_target}")

paired_in_duplex = struct_mirna.count("(") + struct_mirna.count(")")
print(f"Bases paired across the duplex: {paired_in_duplex}")

Step 6: RNA Accessibility Analysis for siRNA Design

Compute the accessibility of a target region within a longer mRNA sequence — critical for siRNA and antisense oligonucleotide efficiency.

import RNA
import numpy as np

def compute_accessibility(mrna_seq: str, window: int = 40) -> list:
    """
    Compute per-position probability of being unpaired (accessible) using
    a sliding-window approach on the partition function.
    Returns list of (position, accessibility) tuples.
    """
    n = len(mrna_seq)
    fc = RNA.fold_compound(mrna_seq)
    _, mfe = fc.mfe()
    fc.exp_params_rescale(mfe)
    fc.pf()
    bpp = fc.bpp()

    # Probability of being paired at each position
    p_paired = np.zeros(n)
    for i in range(1, n + 1):
        for j in range(1, n + 1):
            if i != j:
                p = bpp[min(i,j)][max(i,j)]
                p_paired[i - 1] += p

    p_unpaired = 1.0 - np.clip(p_paired, 0, 1)
    return p_unpaired

# Example: 80-nt mRNA segment with a known accessible region
mrna = "AUGCUAGCUAGCUAGCUAUGCUAGCUAGCUUUUUUUUUUUUAUGCUAGCUAGCUAGCUAGCUAGCUAGCUAGCUAGC"
p_unpaired = compute_accessibility(mrna)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 3))
ax.bar(range(1, len(mrna) + 1), p_unpaired, color="#2166ac", alpha=0.8)
ax.axhline(0.5, color="red", lw=1, ls="--", label="50% unpaired")
ax.set_xlabel("Position (nt)")
ax.set_ylabel("P(unpaired)")
ax.set_title("RNA Accessibility Profile")
ax.legend()
plt.tight_layout()
plt.savefig("rna_accessibility.png", dpi=150, bbox_inches="tight")
print("Saved: rna_accessibility.png")

# Top 5 most accessible positions (siRNA target candidates)
best = sorted(enumerate(p_unpaired, 1), key=lambda x: -x[1])[:5]
print("\nMost accessible positions:")
for pos, prob in best:
    print(f"  Position {pos}: P(unpaired) = {prob:.3f}  ({mrna[pos-1]})")

Step 7: Command-Line RNAfold and Output Parsing

Use the RNAfold CLI for batch folding via subprocess, then parse the output.

# Fold a single sequence from stdin
echo "GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA" | RNAfold

# Output:
# GCGGAUUUAGCUCAGUUGGGAGAGCGCCAGACUGAAGAUCUGGAGGUCCUGUGUUCGAUCCACAGAAUUCGCACCA
# (((((((..((((........)))).(((((.......))))).....(((((.......)))))))))))).... (-31.30)

# Batch fold from FASTA file
RNAfold < sequences.fasta > structures.txt

# Generate base pair probability dot plot (PostScript)
RNAfold --noPS < sequences.fasta   # suppress PostScript output
RNAfold -p < sequences.fasta       # save dot plot as rna.ps

# Python: run RNAfold via subprocess and parse output
import subprocess
import re

def run_rnafold(sequence: str) -> tuple:
    """Run RNAfold CLI and return (structure, mfe) tuple."""
    result = subprocess.run(
        ["RNAfold", "--noPS"],
        input=sequence,
        capture_output=True, text=True, timeout=30
    )
    if result.returncode != 0:
        raise RuntimeError(f"RNAfold failed: {result.stderr}")

    lines = r