RNAlib-2.4.3
Examples

Example programs using the RNAlib C library

Using the 'old' API

The following program exercises most commonly used functions of the library. The program folds two sequences using both the mfe and partition function algorithms and calculates the tree edit and profile distance of the resulting structures and base pairing probabilities.

Note
This program uses the old API of RNAlib, which is in part already marked deprecated. Please consult the RNAlib API v3.0 page for details of what changes are necessary to port your implementation to the new API.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "utils.h"
#include "fold_vars.h"
#include "fold.h"
#include "part_func.h"
#include "inverse.h"
#include "RNAstruct.h"
#include "treedist.h"
#include "stringdist.h"
#include "profiledist.h"
void main()
{
char *seq1="CGCAGGGAUACCCGCG", *seq2="GCGCCCAUAGGGACGC",
*struct1,* struct2,* xstruc;
float e1, e2, tree_dist, string_dist, profile_dist, kT;
Tree *T1, *T2;
swString *S1, *S2;
float *pf1, *pf2;
FLT_OR_DBL *bppm;
/* fold at 30C instead of the default 37C */
temperature = 30.; /* must be set *before* initializing */
/* allocate memory for structure and fold */
struct1 = (char* ) space(sizeof(char)*(strlen(seq1)+1));
e1 = fold(seq1, struct1);
struct2 = (char* ) space(sizeof(char)*(strlen(seq2)+1));
e2 = fold(seq2, struct2);
free_arrays(); /* free arrays used in fold() */
/* produce tree and string representations for comparison */
xstruc = expand_Full(struct1);
T1 = make_tree(xstruc);
S1 = Make_swString(xstruc);
free(xstruc);
xstruc = expand_Full(struct2);
T2 = make_tree(xstruc);
S2 = Make_swString(xstruc);
free(xstruc);
/* calculate tree edit distance and aligned structures with gaps */
tree_dist = tree_edit_distance(T1, T2);
free_tree(T1); free_tree(T2);
printf("%s\n%s %3.2f\n", aligned_line[0], aligned_line[1], tree_dist);
/* same thing using string edit (alignment) distance */
string_dist = string_edit_distance(S1, S2);
free(S1); free(S2);
printf("%s mfe=%5.2f\n%s mfe=%5.2f dist=%3.2f\n",
aligned_line[0], e1, aligned_line[1], e2, string_dist);
/* for longer sequences one should also set a scaling factor for
partition function folding, e.g: */
kT = (temperature+273.15)*1.98717/1000.; /* kT in kcal/mol */
pf_scale = exp(-e1/kT/strlen(seq1));
/* calculate partition function and base pair probabilities */
e1 = pf_fold(seq1, struct1);
/* get the base pair probability matrix for the previous run of pf_fold() */
bppm = export_bppm();
pf1 = Make_bp_profile_bppm(bppm, strlen(seq1));
e2 = pf_fold(seq2, struct2);
/* get the base pair probability matrix for the previous run of pf_fold() */
bppm = export_bppm();
pf2 = Make_bp_profile_bppm(bppm, strlen(seq2));
free_pf_arrays(); /* free space allocated for pf_fold() */
profile_dist = profile_edit_distance(pf1, pf2);
printf("%s free energy=%5.2f\n%s free energy=%5.2f dist=%3.2f\n",
aligned_line[0], e1, aligned_line[1], e2, profile_dist);
}

In a typical Unix environment you would compile this program using:

cc ${OPENMP_CFLAGS} -c example.c -I${hpath}

and link using

cc ${OPENMP_CFLAGS} -o example -L${lpath} -lRNA -lm

where ${hpath} and ${lpath} point to the location of the header files and library, respectively.

Note
As default, the RNAlib is compiled with build-in OpenMP multithreading support. Thus, when linking your own object files to the library you have to pass the compiler specific ${OPENMP_CFLAGS} (e.g. '-fopenmp' for gcc) even if your code does not use openmp specific code. However, in that case the OpenMP flags may be ommited when compiling example.c

Using the 'new' v3.0 API

The following is a simple program that computes the MFE, partition function, and centroid structure by making use of the v3.0 API.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ViennaRNA/eval.h>
#include <ViennaRNA/fold.h>
int main(int argc, char *argv[]){
char *seq = "AGACGACAAGGUUGAAUCGCACCCACAGUCUAUGAGUCGGUGACAACAUUACGAAAGGCUGUAAAAUCAAUUAUUCACCACAGGGGGCCCCCGUGUCUAG";
char *mfe_structure = vrna_alloc(sizeof(char) * (strlen(seq) + 1));
char *prob_string = vrna_alloc(sizeof(char) * (strlen(seq) + 1));
/* get a vrna_fold_compound with default settings */
/* call MFE function */
double mfe = (double)vrna_mfe(vc, mfe_structure);
printf("%s\n%s (%6.2f)\n", seq, mfe_structure, mfe);
/* rescale parameters for Boltzmann factors */
/* call PF function */
FLT_OR_DBL en = vrna_pf(vc, prob_string);
/* print probability string and free energy of ensemble */
printf("%s (%6.2f)\n", prob_string, en);
/* 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 (%6.2f d=%6.2f)\n", cent, vrna_eval_structure(vc, cent), dist);
/* free centroid structure */
free(cent);
/* free pseudo dot-bracket probability string */
free(prob_string);
/* free mfe structure */
free(mfe_structure);
/* free memory occupied by vrna_fold_compound */
return EXIT_SUCCESS;
}

Perl Examples

Using the Flat Interface

Example 1: "Simple MFE prediction"

#!/usr/bin/perl
use warnings;
use strict;
use RNA;
my $seq1 = "CGCAGGGAUACCCGCG";
# compute minimum free energy (mfe) and corresponding structure
my ($ss, $mfe) = RNA::fold($seq1);
# print output
printf "%s [ %6.2f ]\n", $ss, $mfe;

Using the Object Oriented (OO) Interface

The 'fold_compound' class that serves as an object oriented interface for vrna_fold_compound_t

Example 1: "Simple MFE prediction"

#!/usr/bin/perl
use warnings;
use strict;
use RNA;
my $seq1 = "CGCAGGGAUACCCGCG";
# create new fold_compound object
my $fc = new RNA::fold_compound($seq1);
# compute minimum free energy (mfe) and corresponding structure
my ($ss, $mfe) = $fc->mfe();
# print output
printf "%s [ %6.2f ]\n", $ss, $mfe;

Python Examples

Using the Flat Interface

Example 1: "Simple MFE prediction"

import RNA
sequence = "CUCGUCGCCUUAAUCCAGUGCGGGCGCUAGACAUCUAGUUAUCGCCGCAA"
# compute minimum free energy (mfe) and corresponding structure
(structure, mfe) = RNA.fold(sequence)
# print output
print "%s\n%s [ %6.2f ]" % (structure, mfe)

Using the Object Oriented (OO) Interface

The 'fold_compound' class that serves as an object oriented interface for vrna_fold_compound_t

Example 1: "Simple MFE prediction"

import RNA;
sequence = "CGCAGGGAUACCCGCG"
# create new fold_compound object
fc = RNA.fold_compound(sequence)
# compute minimum free energy (mfe) and corresponding structure
(ss, mfe) = fc.mfe()
# print output
print "%s [ %6.2f ]" % (ss, mfe)

Example 2: "Use callback while generating suboptimal structures"

import RNA
sequence = "GGGGAAAACCCC"
# Set global switch for unique ML decomposition
RNA.cvar.uniq_ML = 1
subopt_data = { 'counter' : 1, 'sequence' : sequence }
# Print a subopt result as FASTA record
def print_subopt_result(structure, energy, data):
if not structure == None:
print ">subopt %d" % data['counter']
print "%s" % data['sequence']
print "%s [%6.2f]" % (structure, energy)
# increase structure counter
data['counter'] = data['counter'] + 1
# Create a 'fold_compound' for our sequence
a = RNA.fold_compound(sequence)
# Enumerate all structures 500 dacal/mol = 5 kcal/mol arround
# the MFE and print each structure using the function above
a.subopt_cb(500, print_subopt_result, subopt_data);

Example 3: "Revert MFE to Maximum Matching using soft constraint callbacks"

import RNA
seq1 = "CUCGUCGCCUUAAUCCAGUGCGGGCGCUAGACAUCUAGUUAUCGCCGCAA"
# Turn-off dangles globally
RNA.cvar.dangles = 0
# Data structure that will be passed to our MaximumMatching() callback with two components:
# 1. a 'dummy' fold_compound to evaluate loop energies w/o constraints, 2. a fresh set of energy parameters
mm_data = { 'dummy': RNA.fold_compound(seq1), 'params': RNA.param() }
# Nearest Neighbor Parameter reversal functions
revert_NN = {
RNA.DECOMP_PAIR_HP: lambda i, j, k, l, f, p: - f.eval_hp_loop(i, j) - 100,
RNA.DECOMP_PAIR_IL: lambda i, j, k, l, f, p: - f.eval_int_loop(i, j, k, l) - 100,
RNA.DECOMP_PAIR_ML: lambda i, j, k, l, f, p: - p.MLclosing - p.MLintern[0] - (j - i - k + l - 2) * p.MLbase - 100,
RNA.DECOMP_ML_ML_STEM: lambda i, j, k, l, f, p: - p.MLintern[0] - (l - k - 1) * p.MLbase,
RNA.DECOMP_ML_STEM: lambda i, j, k, l, f, p: - p.MLintern[0] - (j - i - k + l) * p.MLbase,
RNA.DECOMP_ML_ML: lambda i, j, k, l, f, p: - (j - i - k + l) * p.MLbase,
RNA.DECOMP_ML_UP: lambda i, j, k, l, f, p: - (j - i + 1) * p.MLbase,
RNA.DECOMP_EXT_STEM: lambda i, j, k, l, f, p: - f.E_ext_loop(k, l),
RNA.DECOMP_EXT_STEM_EXT: lambda i, j, k, l, f, p: - f.E_ext_loop(i, k),
RNA.DECOMP_EXT_EXT_STEM: lambda i, j, k, l, f, p: - f.E_ext_loop(l, j),
RNA.DECOMP_EXT_EXT_STEM1: lambda i, j, k, l, f, p: - f.E_ext_loop(l, j-1),
}
# Maximum Matching callback function (will be called by RNAlib in each decomposition step)
def MaximumMatching(i, j, k, l, d, data):
return revert_NN[d](i, j, k, l, data['dummy'], data['params'])
# Create a 'fold_compound' for our sequence
fc = RNA.fold_compound(seq1)
# Add maximum matching soft-constraints
fc.sc_add_f(MaximumMatching)
fc.sc_add_data(mm_data, None)
# Call MFE algorithm
(s, mm) = fc.mfe()
# print result
print "%s\n%s (MM: %d)\n" % (seq1, s, -mm)