RNA secondary structure prediction

The ViennaRNA package is a collection of programs for the prediction and analysis of RNA secondary structures. Each of them is interfaced by a separate class in the biotite.application_v2.viennarna subpackage. The most prominent one is RNAfold, that predicts the secondary structure with the minimum free energy for a given sequence. Its interface is the FoldApp. As example we take the famous tRNA-Phe from yeast, which forms the cloverleaf structure found in tRNAs.

import biotite.sequence as seq
import biotite.application_v2.viennarna as viennarna

# Note that RNA sequences use 'T' instead of 'U' in Biotite
trna = seq.NucleotideSequence(
    "GCGGATTTAGCTCAGTTGGGAGAGCGCCAGACTGAAGATCTGG"
    "AGGTCCTGTGTTCGATCCACAGAATTCGCACCA"
)
app = viennarna.FoldApp()
folded = app.run(trna).result()
print(folded.dot_bracket)
print(folded.free_energy)

(((((((..((((........)))).(((((.......))))).....(((((.......))))))))))))....
-22.4

The result is a FoldedRNA, that describes the structure in dot-bracket notation: Each pair of matching brackets is a base pair and each dot is an unpaired base. The four stems of the cloverleaf are clearly visible. For further processing, the structure is also available as an array of base pairs, giving the indices of the paired bases in the sequence.

base_pairs = folded.base_pairs()
print(base_pairs[:5])

[[ 0 71]
 [ 1 70]
 [ 2 69]
 [ 3 68]
 [ 4 67]]

The energy parameters of the prediction depend on the temperature. Hence, by increasing the temperature the structure gets less stable, until it eventually melts.

for temperature in [37, 60, 80]:
    folded_at_temp = app.run(trna, temperature=temperature).result()
    print(f"{temperature} °C: {folded_at_temp.dot_bracket}")
    print(f"       {folded_at_temp.free_energy:6.2f} kcal/mol")

37 °C: (((((((..((((........)))).(((((.......))))).....(((((.......))))))))))))....
       -22.40 kcal/mol
60 °C: ((((((((........((((.((((((..((((...........))))..))))))..))))..))))))))....
        -9.48 kcal/mol
80 °C: ((((.......((((((((.......)))).)))).............(((((.......)))))...))))....
        -1.20 kcal/mol

If some base pairs or unpaired bases are already known, e.g. from experiments, they can be given as constraints via the pairs, paired and unpaired parameters to guide the prediction.

Plotting the structure

The dot-bracket notation is not very intuitive to read. Fortunately, ViennaRNA also provides RNAplot, that computes 2D coordinates for each base, so that the structure can be drawn as a graph. The corresponding PlotApp returns these coordinates as array, that can be plotted with Matplotlib.

import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

coord = viennarna.PlotApp().run(folded.dot_bracket).result()

fig, ax = plt.subplots(figsize=(6.0, 6.0))
# Connect consecutive bases along the backbone
ax.plot(*coord.T, color="black", linewidth=1, zorder=1)
# Connect paired bases
ax.add_collection(
    LineCollection(
        [(coord[i], coord[j]) for i, j in base_pairs],
        color="silver",
        linewidth=1,
        zorder=1,
    )
)
# Draw each base as circle with its symbol, using 'U' for RNA
ax.scatter(*coord.T, s=140, color="white", edgecolor="black", zorder=2)
for (x, y), symbol in zip(coord, str(trna).replace("T", "U")):
    ax.text(x, y, symbol, ha="center", va="center", fontsize=8, zorder=3)
ax.set_aspect("equal")
ax.axis("off")
fig.tight_layout()

Beyond the minimum free energy structure

The minimum free energy structure is only the most stable one among many possible structures. RNAsubopt, interfaced by the SuboptApp, enumerates all structures within a given energy range above the minimum.

suboptimal = viennarna.SuboptApp().run(trna, energy_range=0.5).result()
for folded_alternative in sorted(
    suboptimal, key=lambda folded: folded.free_energy
):
    print(folded_alternative.dot_bracket, folded_alternative.free_energy)

(((((((..((((........)))).(((((.......))))).....(((((.......)))))))))))).... -22.4
((((((.....((((((.((......)).))))))...((((..((((........))))..)))))))))).... -22.1
(((((((..((((........)))).((((((....).))))).....(((((.......)))))))))))).... -22.1
(((((((....((((((.((......)).)))))).((((....))))(((((.......)))))))))))).... -22.1

The other classes of the subpackage cover further programs from ViennaRNA, for example CofoldApp for the structure of two hybridized strands, AlifoldApp for the consensus structure of an alignment or PKplexApp for structures containing pseudoknots.