Minimum Free Energy (MFE)

General Remarks

According to the laws of Thermodynamics, the most probable conformation (here, secondary structure \(s\)) is the one with least free energy (\(E(s)\)), since:

\[p(s) \propto e^{-\beta E(s)} \text{ with } \beta = \frac{1}{RT}\]

where \(R\) is the gas constant and \(T\) the thermodynamic temperature. Hence, a reasonable RNA secondary structure representative is the structure with minimum free energy (MFE). The basic idea behind an efficient dynamic programming MFE prediction algorithm was first described by Zuker and Stiegler [1981]. Our library implements various variants of the original algorithmic idea that can be applied to:

• single sequences,

• multiple sequence alignments (MSA),

• multiple interacting sequences (RNA-RNA interaction), and

• circular sequences

All algorithms are usually accessible in two forms:

• a Simplified MFE Interface that requires minimal input and user interaction, but, at the same time is limited to default setting, and

• a Generic Interface that allows for a fine grained selection of algorithmic details, e.g. parameter settings or constraints

The latter provides functions that operate on vrna_fold_compound_t objects, which, in turn, can be fully customized to the specific needs for the use-case.

See also

The full list of available API symbols for global MFE prediction can be found in Global MFE Prediction.

Simplified MFE Interface

The most simple interface to these implementations might be the vrna_fold() function. It only requires an input RNA sequence in uppercase RNA alphabet (AUCG) to return the MFE in \(kcal \cdot mol^{-1}\) and a corresponding secondary structure in Dot-Bracket Notation.

Single Sequence MFE

Below, you’ll find an example that illustrates how to predict the MFE and a corresponding secondary structure for a single (linear) RNA sequence:

C

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <ViennaRNA/mfe/global.h>
#include <ViennaRNA/utils/basic.h>

int
main()
{
  /* The RNA sequence */
  char  *seq = "GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA";

  /* allocate memory for MFE structure (length + 1) */
  char  *structure = (char *)vrna_alloc(sizeof(char) * (strlen(seq) + 1));

  /* predict Minmum Free Energy and corresponding secondary structure */
  float mfe = vrna_fold(seq, structure);

  /* print sequence, structure and MFE */
  printf("%s\n%s [ %6.2f ]\n", seq, structure, mfe);

  /* cleanup memory */
  free(structure);

  return 0;
}

Python

import RNA

# The RNA sequence
seq = "GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA"

# compute minimum free energy (MFE) and corresponding structure
(ss, mfe) = RNA.fold(seq)

# print output
print("{}\n{} [ {:6.2f} ]".format(seq, ss, mfe))

Perl

use RNA;

# The RNA sequence
my $seq = "GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA";

# compute minimum free energy (MFE) and corresponding structure
my ($ss, $mfe) = RNA::fold($seq);

# print output
printf "%s\n%s [ %6.2f ]\n", $seq, $ss, $mfe;

The result of the above code should look like:

GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA
..(((((........)))))(((((((...)))))))..... [  -8.80 ]

Consensus MFE

In cases where you want to predict the consensus MFE structure for a set of aligned RNA sequences (MSA), consider the following example:

C

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <ViennaRNA/mfe/global.h>
#include <ViennaRNA/utils/basic.h>
#include <ViennaRNA/utils/alignments.h>

int
main()
{
  /* The RNA sequence alignment */
  const char  *sequences[] = {
    "CUGCCUCACAACGUUUGUGCCUCAGUUACCCGUAGAUGUAGUGAGGGU",
    "CUGCCUCACAACAUUUGUGCCUCAGUUACUCAUAGAUGUAGUGAGGGU",
    "---CUCGACACCACU---GCCUCGGUUACCCAUCGGUGCAGUGCGGGU",
    NULL /* indicates end of alignment */
  };

  /* compute the consensus sequence */
  char        *cons = consensus(sequences);

  /* allocate memory for MFE consensus structure (length + 1) */
  char        *structure = (char *)vrna_alloc(sizeof(char) * (strlen(sequences[0]) + 1));

  /* predict Minmum Free Energy and corresponding secondary structure */
  float       mfe = vrna_alifold(sequences, structure);

  /* print consensus sequence, structure and MFE */
  printf("%s\n%s [ %6.2f ]\n", cons, structure, mfe);

  /* cleanup memory */
  free(cons);
  free(structure);

  return 0;
}

Python

import RNA

# The RNA sequence alignment
sequences = [
    "CUGCCUCACAACGUUUGUGCCUCAGUUACCCGUAGAUGUAGUGAGGGU",
    "CUGCCUCACAACAUUUGUGCCUCAGUUACUCAUAGAUGUAGUGAGGGU",
    "---CUCGACACCACU---GCCUCGGUUACCCAUCGGUGCAGUGCGGGU"
]

# compute the consensus sequence
cons = RNA.consensus(sequences)

# predict Minmum Free Energy and corresponding secondary structure
(ss, mfe) = RNA.alifold(sequences);

# print output
print("{}\n{} [ {:6.2f} ]".format(cons, ss, mfe))

Perl

use RNA;

# The RNA sequence alignment
my @sequences = (
    "CUGCCUCACAACGUUUGUGCCUCAGUUACCCGUAGAUGUAGUGAGGGU",
    "CUGCCUCACAACAUUUGUGCCUCAGUUACUCAUAGAUGUAGUGAGGGU",
    "---CUCGACACCACU---GCCUCGGUUACCCAUCGGUGCAGUGCGGGU"
);

# compute the consensus sequence
my $cons = RNA::consensus(\@sequences);

# predict Minmum Free Energy and corresponding secondary structure
my ($ss, $mfe) = RNA::alifold(\@sequences);

# print output
printf "%s\n%s [ %6.2f ]\n", $cons, $ss, $mfe;

Note

Further simplified functions that compute the MFE for more special cases include:

• vrna_circfold() - Treats the RNA sequence as circular sequence (Hofacker and Stadler [2006])

• vrna_alifold() - Uses a Multiple Sequence Alignment as input to compute consensus MFE

• vrna_circalifold() - Treats the aligned sequences as circular sequences

Generic Interface

The generic MFE interface requires one to create a vrna_fold_compound_t object first before calling the vrna_mfe() function (or method). Adaptations to the algorithm parameters can be either introduced through the second or third vrna_fold_compound_t constructor parameter (see also the vrna_fold_compound() function), or later via the various vrna_fold_compound_t associated methods (or functions in C).

The following example demonstrates how to change some of the model parameters before actually computing the MFE (structure) under the non-default settings.

C

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <ViennaRNA/model.h>
#include <ViennaRNA/fold_compound.h>
#include <ViennaRNA/utils/basic.h>
#include <ViennaRNA/mfe/global.h>

int
main()
{
  char      *seq = "GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA";

  /* allocate memory for MFE structure (length + 1) */
  char      *structure = (char *)vrna_alloc(sizeof(char) * (strlen(seq) + 1));

  /* create a new model details structure to store the Model Settings */
  vrna_md_t md;

  /* ALWAYS set default model settings first! */
  vrna_md_set_default(&md);

  /* Mark sequence as being circular, change temperature to
   * non-default and change dangle model
   */
  md.circ         = 1;    /* Treat RNA input as circular sequence */
  md.temperature  = 25.0; /* 25 Deg Celcius */
  md.dangles      = 1;    /* Switch to dangles-model 1 */

  /* create a fold compound */
  vrna_fold_compound_t  *fc = vrna_fold_compound(seq, &md, VRNA_OPTION_DEFAULT);

  /* predict Minmum Free Energy and corresponding secondary structure */
  float                 mfe = vrna_mfe(fc, structure);

  /* print sequence, structure and MFE */
  printf("%s\n%s [ %6.2f ]\n", seq, structure, mfe);

  /* cleanup memory */
  free(structure);
  vrna_fold_compound_free(fc);

  return 0;
}

Python

import RNA

# The RNA sequence
seq = "GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA"

# create a new model details structure with
# non-default settings
md = RNA.md(
        circ        = 1,    # Circular RNA sequence
        temperature = 25.0, # 25 Deg Celcius
        dangles     = 1     # Dangle Model 1
      )

# create a fold compound
fc = RNA.fold_compound(seq, md)

# predict Minmum Free Energy and corresponding secondary structure
(ss, mfe) = fc.mfe()

# print sequence, structure and MFE
print("{}\n{} [ {:6.2f} ]".format(seq, ss, mfe))

Perl

use RNA;

# The RNA sequence
my $seq = "GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA";

# create a new model details structure
my $md = new RNA::md();

# change temperature and dangle model
$md->{circ}         = 1;    # Circular RNA sequence
$md->{temperature}  = 25.0; # 25 Deg Celcius
$md->{dangles}      = 1;    # Dangle Model 1

# create a fold compound
my $fc = new RNA::fold_compound($seq, $md);

# predict Minmum Free Energy and corresponding secondary structure
my ($ss, $mfe) = $fc->mfe();

# print sequence, structure and MFE
printf "%s\n%s [ %6.2f ]\n", $seq, $ss, $mfe;

This should produce the following output:

GAGUAGUGGAACCAGGCUAUGUUUGUGACUCGCAGACUAACA
..(((((........)))))(((((((...)))))))..... [  -5.64 ]

Backtracking

The implementations for mfe/backtracking are generally agnostic with respect to whether local or global structure prediction is in place.

See also

The full list of available API symbols for backtracking of MFE structures can be found in Backtracking MFE structures.

Local MFE Prediction

In addition to the global MFE prediction algorithm(s) above, we also provide implementations for the local (sliding window) MFE prediction, in particular for

• Single sequences,

• Multiple sequence alignments (MSA), and

Note, that our implementation scans an RNA sequence (or MSA) from the 3’ to the 5’ end, and reports back locally optimal (consensus) structures, the corresponding free energy, and the position of the sliding window in global coordinates.

For any particular RNA sequence (or MSA) multiple locally optimal (consensus) secondary structures may be predicted. Thus, we tried to implement an interface that allows for an effortless conversion of the corresponding hits into any target data structure. As a consequence, we provide two distinct ways to retrieve the corresponding predictions, either

• through directly writing to an open FILE stream on-the-fly, or

• through a callback function mechanism.

The latter allows one to store the results in any possible target data structure. Our implementations then pass the results through the user-implemented callback as soon as the prediction for a particular window is finished.

Return Valued Based Example

Here is some example code that shows how to use the simple local MFE prediction API:

C

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <ViennaRNA/mfe/local.h>
#include <ViennaRNA/utils/basic.h>

int
main()
{
  /* a long RNA sequence input */
  char  *seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCA" \
               "GGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUACGUUGGUAUC" \
               "AUGUCCUCCGUCCUAACAAUUAUACAUCGAGAGGCAAAAUUUCUAAUCCGGGGUC" \
               "AGUGAGCAUUGCCAUUUUAUAACUCGUGAUCUCUCGCUACUUAGGCGAUCCCUGC" \
               "CAAUGAGGGUCAAGGAGUUGAAUUAUCGGGCCACAUCGACGUGGCCUUUACGGCC" \
               "AGGUAAUUCAAAGGCCUCAAGUCCU";

  /* allocate memory for MFE structure (length + 1) */
  char  *structure = (char *)vrna_alloc(sizeof(char) * (strlen(seq) + 1));

  /* predict local MFE structures with maximum base pair distance of 50 nt */
  float e = vrna_Lfold(seq, 50, NULL);

  /* cleanup memory */
  free(structure);

  return 0;
}

Python

import RNA

# a long RNA sequence input
longseq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCA" \
          "GGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUACGUUGGUAUC" \
          "AUGUCCUCCGUCCUAACAAUUAUACAUCGAGAGGCAAAAUUUCUAAUCCGGGGUC" \
          "AGUGAGCAUUGCCAUUUUAUAACUCGUGAUCUCUCGCUACUUAGGCGAUCCCUGC" \
          "CAAUGAGGGUCAAGGAGUUGAAUUAUCGGGCCACAUCGACGUGGCCUUUACGGCC" \
          "AGGUAAUUCAAAGGCCUCAAGUCCU"

# predict local MFE structures with maximum base pair distance of 50 nt
e = RNA.Lfold(longseq, 50)

Perl

use RNA;

# a long RNA sequence input
my $longseq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCA" .
              "GGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUACGUUGGUAUC" .
              "AUGUCCUCCGUCCUAACAAUUAUACAUCGAGAGGCAAAAUUUCUAAUCCGGGGUC" .
              "AGUGAGCAUUGCCAUUUUAUAACUCGUGAUCUCUCGCUACUUAGGCGAUCCCUGC" .
              "CAAUGAGGGUCAAGGAGUUGAAUUAUCGGGCCACAUCGACGUGGCCUUUACGGCC" .
              "AGGUAAUUCAAAGGCCUCAAGUCCU";

# predict local MFE structures with maximum base pair distance of 50 nt
my $e = RNA::Lfold($longseq, 50);

In the example above, we didn’t provide an output file stream as third argument for vrna_Lfold(). Hence, the output was automatically written to stdout.

Callback Based Example

In contrast to the previous example, here is a variant that uses the callback based local MFE prediction API:

C

#include <stdlib.h>
#include <stdio.h>
#include <string.h>

#include <ViennaRNA/mfe/local.h>
#include <ViennaRNA/utils/basic.h>

/* a callback function that simply prints the resulting
 * local MFE prediction
 */
void
callback(unsigned int start,
         unsigned int end,
         const char *structure,
         float en,
         void *data)
{
  if (structure) {
    printf("%.*s\n%s (%6.2f) - [%u:%u]\n",
           end - start + 1,
           (char *)data) + start - 1,
           structure,
           en,
           start,
           end);
  }
}

int
main()
{
  /* a long RNA sequence input */
  char  *seq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCA" \
               "GGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUACGUUGGUAUC" \
               "AUGUCCUCCGUCCUAACAAUUAUACAUCGAGAGGCAAAAUUUCUAAUCCGGGGUC" \
               "AGUGAGCAUUGCCAUUUUAUAACUCGUGAUCUCUCGCUACUUAGGCGAUCCCUGC" \
               "CAAUGAGGGUCAAGGAGUUGAAUUAUCGGGCCACAUCGACGUGGCCUUUACGGCC" \
               "AGGUAAUUCAAAGGCCUCAAGUCCU";

  /* predict local MFE structures with maximum base pair distance of 50 nt */
  float e = vrna_Lfold_cb(seq, 50, callback, (void *)seq);

  return 0;
}

Python

import RNA

# a long RNA sequence input
longseq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCA" \
          "GGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUACGUUGGUAUC" \
          "AUGUCCUCCGUCCUAACAAUUAUACAUCGAGAGGCAAAAUUUCUAAUCCGGGGUC" \
          "AGUGAGCAUUGCCAUUUUAUAACUCGUGAUCUCUCGCUACUUAGGCGAUCCCUGC" \
          "CAAUGAGGGUCAAGGAGUUGAAUUAUCGGGCCACAUCGACGUGGCCUUUACGGCC" \
          "AGGUAAUUCAAAGGCCUCAAGUCCU"

def callback(start, end, structure, en, data = None):
    if structure:
        print(f"{data[start - 1:end]}\n{structure} ({en:6.2f}) - [{start}-{end}]")

# predict local MFE structures with maximum base pair distance of 50 nt
e = RNA.Lfold_cb(longseq, 50, callback, longseq)

Perl

use RNA;

# a long RNA sequence input
my $longseq = "AUUUCCACUAGAGAAGGUCUAGAGUGUUUGUCGUUUGUCAGAAGUCCCUAUUCCA" .
              "GGUACGAACACGGUGGAUAUGUUCGACGACAGGAUCGGCGCACUACGUUGGUAUC" .
              "AUGUCCUCCGUCCUAACAAUUAUACAUCGAGAGGCAAAAUUUCUAAUCCGGGGUC" .
              "AGUGAGCAUUGCCAUUUUAUAACUCGUGAUCUCUCGCUACUUAGGCGAUCCCUGC" .
              "CAAUGAGGGUCAAGGAGUUGAAUUAUCGGGCCACAUCGACGUGGCCUUUACGGCC" .
              "AGGUAAUUCAAAGGCCUCAAGUCCU";

sub callback {
  my ($start, $end, $structure, $en, $data) = @_;
  if ($structure) {
    printf(substr($data, start - 1, end - start + 1) .
           $structure .
           "(%6.2f) - [%u:%u]",
           en,
           start,
           end);
  }
}

# predict local MFE structures with maximum base pair distance of 50 nt
my $e = RNA::Lfold_cb($longseq, 50, \&callback, \&longseq);

See also

The full list of available API symbols for local (sliding-window) MFE prediction can be found in Local MFE Prediction.