Suboptimals and Representative Structures#
Suboptimal Structures sensu Zuker#
The algorithm to compute optimal secondary structures that contain a particular base pair has been published by Zuker [1989] and is based on ideas for predicting structures for circular RNAs, in particular viroids, as presented in Steger et al. [1984].
Suboptimal Structures within an Energy Band around the MFE#
Callback Based Suboptimals#
#include <stdlib.h>
#include <stdio.h>
#include <ViennaRNA/fold_compound.h>
#include <ViennaRNA/utils/basic.h>
#include <ViennaRNA/utils/strings.h>
#include <ViennaRNA/subopt/basic.h>
void
subopt_callback(const char *structure,
float energy,
void *data)
{
/* simply print the result and increase the counter variable by 1 */
if (structure)
printf("%04d\t%s\t%6.2f\n", (*((int *)data))++, structure, energy);
}
int
main()
{
char *sequence = "GGGGAAAACCCC"
/* Create a fold compound for the sequence */
vrna_fold_compound_t *fc = vrna_fold_compound(seq, NULL, VRNA_OPTION_DEFAULT);
int counter = 0;
/*
* call subopt to enumerate all secondary structures in an energy band of
* 5 kcal/mol of the MFE and pass it the address of the callback and counter
* variable
*/
vrna_subopt_cb(fc, 500, &subopt_callback, (void *)&counter);
/* cleanup memory */
free(seq);
vrna_fold_compound_free(fc);
return 0;
}
import RNA
sequence = "GGGGAAAACCCC"
counter = 0
# Print a subopt result and increase counter by 1
def print_subopt_result(structure, energy, data):
if structure != None:
data += 1
print(f"{data:04d}\t{structure}\t{energy:6.2f}")
# Get model details with unique ML decomposition
md = RNA.md(uniq_ML = 1)
# Create a 'fold_compound' for our sequence
fc = RNA.fold_compound(sequence, md)
# Enumerate all structures 500 dacal/mol = 5 kcal/mol arround
# the MFE and print each structure using the function above
fc.subopt_cb(500, print_subopt_result, counter);
Boltzmann Sampling#
Below follow a few Python examples that demonstrate the usage of the
Boltzmann sampling implementations. To prepare the Boltzmann sampling
implementation, we first need to compute the partition function and thereby
fill the dynamic programming matrices required for the stochastic backtracking:
Prerequisites#
import RNA
sequence = "UGGGAAUAGUCUCUUCCGAGUCUCGCGGGCGACGGGCGAUCUUCGAAAGU \
GGAAUCCGUACUUAUACCGCCUGUGCGGACUACUAUCCUGACCACAUAGU"
"""
First we prepare a fold_compound object and run the usual
steps to compute the partition function
"""
# create model details
md = RNA.md()
# activate unique multibranch loop decomposition
md.uniq_ML = 1
# create fold compound object
fc = RNA.fold_compound(sequence, md)
# compute MFE
(ss, mfe) = fc.mfe()
# rescale Boltzmann factors according to MFE
fc.exp_params_rescale(mfe)
# compute partition function to fill DP matrices
fc.pf()
"""
Now we are ready to perform Boltzmann sampling
"""
Let us now step by step go through the different ways to create samples of (sub)structures from the ensemble.
Sampling Substructures#
We begin with sub-structures of a particular length that start at the 5’ end, i.e. pos 1, of the input sequence:
# 1. backtrace a single sub-structure of length 10
print(fc.pbacktrack5(10))
# 2. backtrace a single sub-structure of length 50
print(fc.pbacktrack5(50))
# 3. backtrace multiple sub-structures of length 10 at once
for s in fc.pbacktrack5(20, 10):
print(s)
# 4. backtrace multiple sub-structures of length 50 at once
for s in fc.pbacktrack5(100, 50):
print(s)
def store_structure(structure, data):
"""
A simple callback function that stores
a structure sample into a list
"""
if structure:
data.append(structure)
# Draw a sample of 100 sub-structures of length 50
ss = []
i = fc.pbacktrack5(100, 50, store_structure, ss)
for s in ss:
print(s)
Of course one can sample substructures from anywhere within the input sequence as well:
# 1. backtrace a single substructure from the subsequence [10:50]
print(fc.pbacktrack_sub(10, 50))
# 2. backtrace multiple substructures from the subsequence [10:50]
for s in fc.pbacktrack_sub(100, 10, 50):
print(s)
# 3. backtrace multiple substructures from the subsequence [10:50] non-redundantly
for s in fc.pbacktrack_sub(100,
10,
50,
RNA.PBACKTRACK_NON_REDUNDANT):
print(s)
def store_structure(structure, data):
"""
A simple callback function that stores
a structure sample into a list
"""
if structure:
data.append(structure)
# list container to store our structures
ss = []
# 1. backtrace multiple substructures from the subsequence [10:50]
i = fc.pbacktrack_sub(100, 10, 50, store_structure, ss)
# 2. backtrace multiple substructures from the subsequence [10:50] non-redundantly
i = fc.pbacktrack_sub(100,
10,
50,
store_structure,
ss,
RNA.PBACKTRACK_NON_REDUNDANT)
for s in ss:
print(s)
Sampling Full-length Structures#
Finally, let us draw full-length structures:
# 1. backtrace a single structure (full length)
print(fc.pbacktrack())
# 2. backtrace multiple structures at once
for s in fc.pbacktrack(100):
print(s)
# 3. backtrace multiple structures non-redundantly
for s in fc.pbacktrack(100, RNA.PBACKTRACK_NON_REDUNDANT):
print(s)
def store_structure(structure, data):
"""
A simple callback function that stores
a structure sample into a list
"""
if structure:
data.append(structure)
# list container to store our structures
ss = []
# 1. backtrace multiple full-length structures in callback mode
i = fc.pbacktrack(100, store_structure, ss)
# 2. non-redundantly backtrace multiple full-length structures in callback mode
i = fc.pbacktrack(100,
store_structure,
ss,
RNA.PBACKTRACK_NON_REDUNDANT)
for s in ss:
print(s)
Non-redundant Sampling#
The non-redundant mode for Boltzmann sampling (activated by
RNA.PBACKTRACK_NON_REDUNDANT) removes each sample that has been drawn
from the ensemble. Thus, no structure will be sampled more than once.
Whenever you drew a few structures in non-redundant mode and want to continue sampling from the remaining sub-ensemble that is left, you can use the resume capability for this sampling process:
# perform 15 iterations and draw 500 non-redundant structures each time
# (with resume option)
num_samples = 500
iterations = 15
d = None # pbacktrack memory object
s_list = []
for i in range(0, iterations):
d, ss = fc.pbacktrack(num_samples, d, RNA.PBACKTRACK_NON_REDUNDANT)
s_list.extend(ss)
for s in s_list:
print(s)
def store_structure(structure, data):
"""
A simple callback function that stores
a structure sample into a list
"""
if structure:
data.append(structure)
# Non-redundantly backtrace multiple full length structures
# in callback mode with resume option
num_samples = 500
iterations = 15
d = None # pbacktrack memory object
s_list = []
for i in range(0, iterations):
d, i = fc.pbacktrack(num_samples,
store_structure,
s_list,
d,
RNA.PBACKTRACK_NON_REDUNDANT)
for s in ss:
print(s)
Maximum Expected Accuracy (MEA) Structure#
Another useful structure representative one can determine from base pairing probabilities \(p_{ij}\) is the structure that exhibits the maximum expected accuracy (MEA). By assuming the base pair probability is a good measure of correctnes of a pair \((i,j)\), the expected accuracy of a structure is
\[EA(s) = \sum_{(i,j) \in s} 2 \gamma p_{ij} + \sum_{i \notin s} p^u_i\]
with unpaired probabilities \(p^u_i = 1 - \sum_{j \neq i} p_{ij}\) and weighting factor \(\gamma\) that allows us to weight paired against unpaired positions. The MEA algorithm attempts to find the structure that maximizes \(EA\). Higher values of \(\gamma\) result in more base pairs of lower probability and thus higher sensitivity. Low values of \(\gamma\) result in structures containing only highly likely pairs, i.e. high specificity.
Example MEA structure Prediction
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ViennaRNA/data_structures.h>
#include <ViennaRNA/params/basic.h>
#include <ViennaRNA/utils/basic.h>
#include <ViennaRNA/eval/structures.h>
#include <ViennaRNA/mfe/global.h>
#include <ViennaRNA/partfunc/basic.h>
#include <ViennaRNA/structures/mea.h>
int
main(int argc,
char *argv[])
{
char *seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA";
/* get a vrna_fold_compound with default settings */
vrna_fold_compound_t *vc = vrna_fold_compound(seq, NULL, VRNA_OPTION_DEFAULT);
/* call MFE function */
double mfe = (double)vrna_mfe(vc, NULL);
/* rescale parameters for Boltzmann factors */
vrna_exp_params_rescale(vc, &mfe);
/* call PF function */
FLT_OR_DBL en = vrna_pf(vc, NULL);
/* compute centroid structure */
float mea;
char *mea_structure = vrna_MEA(fc, 1.0, &mea);
/* print MEA structure, its free energy and maximum expected accuracy */
printf("%s\n%s (%6.2f MEA=%6.2f)\n",
seq,
mea_structure,
vrna_eval_structure(vc, mea_structure),
mea);
/* free MEA structure */
free(mea_structure);
/* free memory occupied by vrna_fold_compound */
vrna_fold_compound_free(vc);
return EXIT_SUCCESS;
}
#!/usr/bin/python
#
import RNA
seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA"
# create fold_compound data structure
# (required for all subsequently applied algorithms)
fc = RNA.fold_compound(seq)
# compute MFE and MFE structure
(mfe_struct, mfe) = fc.mfe()
# rescale Boltzmann factors for partition function computation
fc.exp_params_rescale(mfe)
# compute partition function
(pp, pf) = fc.pf()
# compute MEA structure
(mea_struct, mea) = fc.MEA(1.0)
# compute free energy of MEA structure
mea_en = fc.eval_structure(mea_struct)
# print only the MEA structure, its free energy
# in kcal/mol and the maximum exprected accuracy
print(f"{seq}\n{mea_struct} ({mea_en:6.2f} d={mea:.2f})")
#!/usr/bin/perl
#
use RNA;
my $seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA";
# create fold_compound data structure
# (required for all subsequently applied algorithms)
my $fc = new RNA::fold_compound($seq);
# compute MFE and MFE structure
my ($mfe_struct, $mfe) = $fc->mfe();
# rescale Boltzmann factors for partition function computation
$fc->exp_params_rescale($mfe);
# compute partition function
my ($pp, $pf) = $fc->pf();
# compute MEA structure
my ($mea_struct, $mea) = $fc->MEA(1.0);
# compute free energy of MEA structure
my $mea_en = $fc->eval_structure($mea_struct);
# print only the MEA structure, its free energy
# in kcal/mol and the maximum expected accuracy
printf("%s\n%s (%6.2f MEA = %.2f)\n", $seq, $mea_struct, $mea_en, $mea);
The examples above should all print the following:
AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA
.....((((..(((...)))..))))((((.(....).)))).((.(((.....))).)).. ( -8.30 MEA = 47.47)
The Centroid Structure#
The centroid is the structure with the minimal average distance to all other structures
\[<d(s)> = \sum_{(i,j) \in s} (1-p_{ij}) + \sum_{(i,j) \notin s} p_{ij}\]
Thus, the centroid is simply the structure containing all pairs with \(p_ij > 0.5\)
Example Centroid structure Prediction
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ViennaRNA/data_structures.h>
#include <ViennaRNA/params/basic.h>
#include <ViennaRNA/utils/basic.h>
#include <ViennaRNA/eval/structures.h>
#include <ViennaRNA/mfe/global.h>
#include <ViennaRNA/partfunc/basic.h>
#include <ViennaRNA/structures/centroid.h>
int
main(int argc,
char *argv[])
{
char *seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA";
/* get a vrna_fold_compound with default settings */
vrna_fold_compound_t *vc = vrna_fold_compound(seq, NULL, VRNA_OPTION_DEFAULT);
/* call MFE function */
double mfe = (double)vrna_mfe(vc, NULL);
/* rescale parameters for Boltzmann factors */
vrna_exp_params_rescale(vc, &mfe);
/* call PF function */
FLT_OR_DBL en = vrna_pf(vc, NULL);
/* compute centroid structure */
double dist;
char *cent = vrna_centroid(vc, &dist);
/* print centroid structure, its free energy and mean distance to the ensemble */
printf("%s\n%s (%6.2f d=%6.2f)\n",
seq,
cent,
vrna_eval_structure(vc, cent),
dist);
/* free centroid structure */
free(cent);
/* free memory occupied by vrna_fold_compound */
vrna_fold_compound_free(vc);
return EXIT_SUCCESS;
}
#!/usr/bin/python
#
import RNA
seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA"
# create fold_compound data structure
# (required for all subsequently applied algorithms)
fc = RNA.fold_compound(seq)
# compute MFE and MFE structure
(mfe_struct, mfe) = fc.mfe()
# rescale Boltzmann factors for partition function computation
fc.exp_params_rescale(mfe)
# compute partition function
(pp, pf) = fc.pf()
# compute centroid structure
(centroid_struct, dist) = fc.centroid()
# compute free energy of centroid structure
centroid_en = fc.eval_structure(centroid_struct)
# print only the centroid structure, its free energy
# in kcal/mol and the distance to the ensemble
print(f"{seq}\n{centroid_struct} ({centroid_en:6.2f} d={dist:.2f})")
#!/usr/bin/perl
#
use RNA;
my $seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA";
# create fold_compound data structure
# (required for all subsequently applied algorithms)
my $fc = new RNA::fold_compound($seq);
# compute MFE and MFE structure
my ($mfe_struct, $mfe) = $fc->mfe();
# rescale Boltzmann factors for partition function computation
$fc->exp_params_rescale($mfe);
# compute partition function
my ($pp, $pf) = $fc->pf();
# compute centroid structure
my ($centroid_struct, $dist) = $fc->centroid();
# compute free energy of centroid structure
my $centroid_en = $fc->eval_structure($centroid_struct);
# print only the centroid structure, its free energy
# in kcal/mol and the distance to the ensemble
printf("%s\n%s (%6.2f d = %.2f)\n", $seq, $centroid_struct, $centroid_en, $dist);
The examples above should all print the following:
AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCAGGUACGA
.....((((..((.....))..))))((((........)))).((.(((.....))).)).. ( -7.20 d = 8.96)