PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of amolbioBioMed CentralBiomed Central Web Sitesearchsubmit a manuscriptregisterthis articleAlgorithms for Molecular Biology : AMB
 
Algorithms Mol Biol. 2010; 5: 27.
Published online 2010 June 23. doi:  10.1186/1748-7188-5-27
PMCID: PMC2909241

Inverse folding of RNA pseudoknot structures

Abstract

Background

RNA exhibits a variety of structural configurations. Here we consider a structure to be tantamount to the noncrossing Watson-Crick and G-U-base pairings (secondary structure) and additional cross-serial base pairs. These interactions are called pseudoknots and are observed across the whole spectrum of RNA functionalities. In the context of studying natural RNA structures, searching for new ribozymes and designing artificial RNA, it is of interest to find RNA sequences folding into a specific structure and to analyze their induced neutral networks. Since the established inverse folding algorithms, RNAinverse, RNA-SSD as well as INFO-RNA are limited to RNA secondary structures, we present in this paper the inverse folding algorithm Inv which can deal with 3-noncrossing, canonical pseudoknot structures.

Results

In this paper we present the inverse folding algorithm Inv. We give a detailed analysis of Inv, including pseudocodes. We show that Inv allows to design in particular 3-noncrossing nonplanar RNA pseudoknot 3-noncrossing RNA structures-a class which is difficult to construct via dynamic programming routines. Inv is freely available at http://www.combinatorics.cn/cbpc/inv.html.

Conclusions

The algorithm Inv extends inverse folding capabilities to RNA pseudoknot structures. In comparison with RNAinverse it uses new ideas, for instance by considering sets of competing structures. As a result, Inv is not only able to find novel sequences even for RNA secondary structures, it does so in the context of competing structures that potentially exhibit cross-serial interactions.

1 Introduction

Pseudoknots are structural elements of central importance in RNA structures [1], see Figure Figure1.1. They represent cross-serial base pairing interactions between RNA nucleotides that are functionally important in tRNAs, RNaseP [2], telomerase RNA [3], and ribosomal RNAs [4]. Pseudoknot structures are being observed in the mimicry of tRNA structures in plant virus RNAs as well as the binding to the HIV-1 reverse transcriptase in in vitro selection experiments [5]. Furthermore basic mechanisms, like ribosomal frame shifting, involve pseudoknots [6].

Figure 1
Representations of RNA structures. The pseudoknot structure of the glmS ribozyme pseudoknot P1.1 [40] as a diagram (top) and as a planar graph (bottom).

Despite them playing a key role in a variety of contexts, pseudoknots are excluded from large-scale computational studies. Although the problem has attracted considerable attention in the last decade, pseudoknots are considered a somewhat "exotic" structural concept. For all we know [7], the ab initio prediction of general RNA pseudoknot structures is NP-complete and algorithmic difficulties of pseudoknot folding are confounded by the fact that the thermodynamics of pseudoknots is far from being well understood.

As for the folding of RNA secondary structures, Waterman et al [8,9], Zuker et al [10] and Nussinov [11] established the dynamic programming (DP) folding routines. The first mfe-folding algorithm for RNA secondary structures, however, dates back to the 60's [12-14]. For restricted classes of pseudoknots, several algorithms have been designed: Rivas and Eddy [15], Dirks and Pierce [16], Reeder and Giegerich [17] and Ren et al [18]. Recently, a novel ab initio folding algorithm Cross has been introduced [19]. Cross generates minimum free energy (mfe), 3-noncrossing, 3-canonical RNA structures, i.e. structures that do not contain three or more mutually crossing arcs and in which each stack, i.e. sequence of parallel arcs, see eq. (1), has size greater or equal than three. In particular, in a 3-canonical structure there are no isolated arcs, see Figure Figure22.

Figure 2
σ-canonical RNA structure. Each stack of "parallel" arcs has to have minimum size σ. Here we display a 3-canonical structure.

The notion of mfe-structure is based on a specific concept of pseudoknot loops and respective loop-based energy parameters. This thermodynamic model was conceived by Tinoco and refined by Freier, Turner, Ninio, and others [13,20-24].

1.1 k-noncrossing, σ-canonical RNA pseudoknot structures

Let us turn back the clock: three decades ago Waterman et al. [25], Nussinov et al. [11] and Kleitman et al. in [26] analyzed RNA secondary structures. Secondary structures are coarse grained RNA contact structures, see Figure Figure33.

Figure 3
The phenylalanine tRNA structure. The phenylalanine tRNA secondary structure represented as 2-noncrossing diagram (top) and as planar graph (bottom).

RNA secondary structures as well as RNA pseudoknot structures can be represented as diagrams, i.e. labeled graphs over the vertex set [n] = {1, ..., n} with vertex degrees ≤ 1, represented by drawing its vertices on a horizontal line and its arcs (i, j) (i < j), in the upper half-plane, see Figure Figure44 and Figure Figure1.1. Given an arc (i, j) we refer to (j - i) as its arc-length.

Figure 4
Secondary structure. Secondary structures are particular k-noncrossing diagrams, 2-noncrossing diagrams exhibit no crossings at all, therefore RNA secondary structures coincide with 2-noncrossing diagrams having minimum arc-length two.

Here, vertices and arcs correspond to the nucleotides A, G, U, C and Watson-Crick (A-U, G-C) and (U-G) base pairs, respectively.

In a diagram, two arcs (i1, j1) and (i2, j2) are called crossing if i1 <i2 <j1 <j2 holds. Accordingly, a k-crossing is a sequence of arcs (i1, j1), ..., (ik, jk) such that i1 <i2 < ... <ik <j1 <j2 < ... <jk. We call diagrams containing at most (k - 1)-crossings, k-noncrossing diagrams, see Figure Figure55.

Figure 5
k-noncrossing diagrams. We display a 4-noncrossing diagram containing the three mutually crossing arcs (1, 7), (4, 9), (5, 11) (drawn in red).

RNA secondary structures exhibit no crossings in their diagram representation, see Figure Figure33 and Figure Figure4,4, and are therefore 2-noncrossing diagrams satisfying some minimum arc-length condition. An RNA pseudoknot structure is therefore a k-noncrossing diagram for some k satisfying some minimum arc-length condition.

A structure in which any stack has at least size σ is called σ-canonical, where a stack of size σ is a sequence of "parallel" arcs of the form

equation image
(1)

A sequence of consecutive stacks, separated by unpaired nucleotides, An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i2.gif i.e. where

equation image

is called a stem of length r, see Figure Figure66.

Figure 6
Stems. A stem composed by a sequence of three nested stacks. Note that respective stacks only have to be separated by isolated nucleotides on either the left hand side or the right hand side but not necessarily both.

As a natural generalization of RNA secondary structures k-noncrossing RNA structures [27-29] were introduced. A k-noncrossing RNA structure of length n is k-noncrossing diagram over [n] without arcs of the form (i, i + 1). In the following we assume k = 3, i.e. in the diagram representation there are at most two mutually crossing arcs, a minimum arc-length of four and a minimum stack-size of three base pairs. The notion k-noncrossing stipulates that the complexity of a pseudoknot is related to the maximal number of mutually crossing bonds. Indeed, most natural RNA pseudoknots are 3-noncrossing [30].

1.2 Neutral networks

Before considering an inverse folding algorithm into specific RNA structures one has to have at least some rationale as to why there exists one sequence realizing a given target as mfe-configuration. In fact this is, on the level of entire folding maps, guaranteed by the combinatorics of the target structures alone. It has been shown in [31], that the numbers of 3-noncrossing RNA pseudoknot structures, satisfying the biophysical constraints grows asymptotically as c3n-52.03n, where c3 >0 is some explicitly known constant. In view of the central limit theorems of [32], this fact implies the existence of extended (exponentially large) sets of sequences that all fold into one 3-noncrossing RNA pseudoknot structure, S. In other words, the combinatorics of 3-noncrossing RNA structures alone implies that there are many sequences mapping (folding) into a single structure. The set of all such sequences is called the neutral network of the structure S [33,34], see Figure Figure7.7. The term "neutral network" as opposed to "neutral set" stems from giant component results of random induced subgraphs of n-cubes. That is, neutral networks are typically connected in sequence space.

Figure 7
Neutral network in sequence space. We display sequence space (left) and structure space (right) as grids. We depict a set of sequences that all fold into a particular structure. Any two of these sequences are connected by a red path. The neutral network ...

By construction, all the sequences contained in such a neutral network are all compatible with S. That is, at any two positions paired in S, we find two bases capable of forming a bond (A-U, U-A, G-C, C-G, G-U and U-G), see Figure Figure8.8. Let s' be a sequence derived via a point-mutation of s. If s' is again compatible with S, we call this mutation "compatible".

Figure 8
A structure and a particular compatible sequence. A structure and a particular compatible sequence organized in the segments of unpaired and paired bases.

Let C[S] denote the set of S-compatible sequences. The structure S motivates to consider a new adjacency relation within C[S]. Indeed, we may reorganize a sequence (s1, ..., sn) into the pair

equation image
(2)

where the uh denotes the unpaired nucleotides and the ph = (si, sj) denotes base pairs, respectively, see Figure Figure8.8. We can then view An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i5.gif and An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i6.gif as elements of the formal cubes An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i7.gif and An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i8.gif implying the new adjacency relation for elements of C[S].

Accordingly, there are two types of compatible neighbors in the sequence space u- and p-neighbors: a u-neighbor has Hamming distance one and differs exactly by a point mutation at an unpaired position. Analogously a p-neighbor differs by a compensatory base pair-mutation, see Figure Figure99.

Figure 9
Diagram representation of an RNA structure and its compatible neighbors. Diagram representation of an RNA structure (top) and its induced compatible neighbors in sequence space (bottom). Here the neighbors on the inner circle have Hamming distance one ...

Note, however, that a p-neighbor has either Hamming distance one (G-C [mapsto] G-U) or Hamming distance two (G-C [mapsto] C-G). We call a u- or a p-neighbor, y, a compatible neighbor. In light of the adjacency notion for the set of compatible sequences we call the set of all sequences folding into S the neutral network of S. By construction, the neutral network of S is contained in C[S]. If y is contained in the neutral network we refer to y as a neutral neighbor. This gives rise to consider the compatible and neutral distance of the two sequences, denoted by C(s, s') and N(s, s'). These are the minimum length of a C[S]-path and path in the neutral network between s and s', respectively. Note that since each neutral path is in particular a compatible path, the compatible distance is always smaller or equal than the neutral distance.

In this paper we study the inverse folding problem for RNA pseudoknot structures: for a given 3-noncrossing target structure S, we search for sequences from C[S], that have S as mfe configuration.

2 Background

For RNA secondary structures, there are three different strategies for inverse folding, RNAinverse, RNA-SSD and INFO-RNA[35-37].

They all generate via a local search routine iteratively sequences, whose structures have smaller and smaller distances to a given target. Here the distance between two structures is obtained by aligning them as diagrams and counting "0", if a given position is either unpaired or incident to an arc contained in both structures and "1", otherwise, see Figure Figure1010.

Figure 10
Distance of two structures. Positions paired differently in S1 and S2 are assigned a "1". There are two types of positions: I. p is contained in different arcs, see position 4, (4, 20) [set membership] S1 and (4, 17) [set membership] S2. II. p is unpaired in one structure ...

One common assumption in these inverse folding algorithms is, that the energies of specific substructures contribute additively to the energy of the entire structure. Let us proceed by analyzing the algorithms.

RNAinverse is the first inverse-folding algorithm that derives sequences that realize given RNA secondary structures as mfe-configuration. In its initialization step, a random compatible sequence s for the target T is generated. Then RNAinverse proceeds by updating the sequence s to s', s'' ... step by step, minimizing the structure distance between the mfe structure of s' and the target structure T. Based on the observation, that the energy of a substructure contributes additively to the mfe of the molecule, RNAinverse optimizes "small" substructures first, eventually extending these to the entire structure. While optimizing substructures, RNAinverse does an adaptive walk in order to decrease the structure distance. In fact, this walk is based entirely on random compatible mutations.

RNA-SSD inverse folds RNA secondary structures by initializing sequences using three specific subroutines. In the first a particular compatible sequence is generated, where non-complementary nucleotides to bases adjacent to helical regions are assigned. In the second nucleotides located in unpaired positions as well as helical regions are assigned at random, using specific (non-uniform) probabilities. The third routine constitutes a mechanism for minimizing the occurrence of undesired but favourable interactions between specific sequence segments. Following these subroutines, RNA-SSD derives a hierarchical decomposition of the target structure. It recursively splits the structure and thereby derives a binary decomposition tree rooted in T and whose leaves correspond to T-substructures. Each non-leaf node of this tree represents a substructure obtained by merging the two substructures of its respective children. Given this tree, RNA-SSD performs a stochastic local search, starting at the leaves, subsequently working its way up to the root.

INFO-RNA constructs sequences folding into a given secondary structure by employing a dynamic programming method for finding a well suited initial sequence. This sequence has a lowest energy with respect to the T. Since the latter does not necessarily fold into T, (due to potentially existing competing configurations) INFO-RNA then utilizes an improved (relative to the local search routine used in RNAinverse) stochastic local search in order to find a sequence in the neutral network of T. In contrast to RNAinverse, INFO-RNA allows for increasing the distance to the target structure. At the same time, only positions that do not pair correctly and positions adjacent to these are examined.

2.1 Cross

Cross is an ab initio folding algorithm that maps RNA sequences into 3-noncrossing RNA structures. It is guaranteed to search all 3-noncrossing, σ-canonical structures and derives some (not necessarily unique), loop-based mfe-configuration. In the following we always assume σ ≥ 3. The input of Cross is an arbitrary RNA sequence s and an integer N. Its output is a list of N 3-noncrossing, σ-canonical structures, the first of which being the mfe-structure for s. This list of N structures (C0, C1, ..., CN-1) is ordered by the free energy and the first list-element, the mfe-structure, is denoted by Cross(s). If no N is specified, Cross assumes N = 1 as default.

Cross generates a mfe-structure based on specific loop-types of 3-noncrossing RNA structures. For a given structure S, let α be an arc contained in S (S-arc) and denote the set of S-arcs that cross α by An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i9.gif. For two arcs α = (i, j) and α' = (i', j'), we next specify the partial order "[precedes]" over the set of arcs:

equation image

All notions of minimal or maximal elements are understood to be with respect to [precedes]. An arc α [set membership] An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i11.gif is called a minimal, β-crossing if there exists no α' [set membership] An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i11.gif such that α' [precedes] α. Note that α [set membership] An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i11.gif can be minimal β-crossing, while β is not minimal α-crossing. 3-noncrossing diagrams exhibit the following four basic loop-types:

(1) A hairpin-loop is a pair

equation image

where (i, j) is an arc and [i, j] is an interval, i.e. a sequence of consecutive, isolated vertices (i, i + 1, ..., j - 1, j).

(2) An interior-loop, is a sequence

equation image

where (i2, j2) is nested in (i1, j1). That is we have i1 < i2 < j2 < j1.

(3) A multi-loop, see Figure Figure1111[19], is the closed structure formed by

equation image
(3)
Figure 11
The standard loop-types. The standard loop-types: hairpin-loop (top), interior-loop (middle) and multi-loop (bottom). These represent all loop-types that occur in RNA secondary structures.

where An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i15.gif denotes the substructure over the interval [ωh, τh], subject to the condition that if all these substructures are simply stems, then there are at least two of them, see Figure Figure66.

A pseudoknot, see Figure Figure1212[19], consists of the following data:

Figure 12
Pseudoknots. Pseudoknot loops, formed by all blue vertices and arcs.

(P1) A set of arcs

equation image

where i1 = min{ih} and jt = max{jh}, such that

(i) the diagram induced by the arc-set P is irreducible, i.e. the dependency-graph of P (i.e. the graph having P as vertex set and in which α and α' are adjacent if and only if they cross) is connected and

(ii) for each (ih, jh) [set membership] P there exists some arc β (not necessarily contained in P) such that (ih, jh) is minimal β-crossing.

(P2) Any i1 <x <jt, not contained in hairpin-, interior- or multi-loops.

Having discussed the basic loop-types, we are now in position to state

Theorem 1 Any 3-noncrossing RNA pseudoknot structure has a unique loop-decomposition [19].

Figure Figure1313 illustrates the loop decomposition of a 3-noncrossing structure.

Figure 13
Loop decomposition. Here a hairpin-loop (I), an interior-loop (II), a multi-loop (III) and a pseudoknot (IV).

In order to discuss the organization of Cross, we introduce the basic idea behind motifs and skeleta, combinatorial structures used in the folding algorithm.

A motif is a 3-noncrossing structure, having only [precedes]-maximal stacks of size exactly σ, i.e. no stacks nested in other stacks, see Figure Figure14.14. Despite that motifs can exhibit complicated crossings, they can be inductively generated. A skeleton, S is a k-noncrossing structure such that

Figure 14
Motif. A 3-noncrossing, 3-canonical motif.

• its core, c(S) has no noncrossing arcs and

• its L-graph, L(S) is connected.

Here the core of a structure, c(S), is obtained by collapsing its stacks into single arcs (thereby reducing its length) and the graph L(S) is obtained by mapping arcs into vertices and connecting any two if they cross in the diagram representation of S, see Figure Figure15.15. A skeleton reflects all cross-serial interactions of a structure.

Figure 15
Skeleton and its L-graph. We display a skeleton (left) and its L-graph (right).

Having introduced motifs and skeleta we can proceed by discussing the general idea of Cross. The algorithm generates 3-noncrossing RNA structure "from top to bottom" via the following three subroutines:

I (SHADOW): In this routine we generate all maximal stacks of the structure. Note that a stack is maximal with respect to [precedes] if it is not nested in some other stack. This is derived by "shadowing" the motifs, i.e. their σ-stacks are extended "from top to bottom".

II (SKELETONBRANCH): Given a shadow, the second step of Cross consists in generating, the skeleta-tree. The nodes of this tree are particular 3-noncrossing structures, obtained by successive insertions of stacks. Intuitively, a skeleton encapsulates all cross-serial arcs that cannot be recursively computed. Here the tree complexity is controlled via limiting the (total) number of pseudoknots.

III (SATURATION): In the third subroutine each skeleton is saturated via DP-routines. After the saturation the mfe-3-noncrossing structure is derived.

Figure Figure1616 provides an overview on how the three subroutines are combined.

Figure 16
An outline of Cross. For illustration purposes we assume here σ = 1. The routines SHADOW, SKELETONBRANCH and SATURATION are depicted. Due to space limitations we only represent a few select motifs and for the same reason only one of the motifs ...

3 The algorithm

The inverse folding algorithm Inv is based on the ab initio folding algorithm Cross. The input of Inv is the target structure, T. The latter is expressed as a character string of ":( )[ ]{ }", where ":" denotes unpaired base and "( )", "[ ]", "{ }" denote paired bases.

In Algorithm 7.1, we present the pseudocodes of algorithm Inv. After validation of the target structure (lines 2 to 5 in Algorithm 7.1), similar to INFO-RNA, Inv constructs an initial sequence and then proceeds by a stochastic local search based on the loop decomposition of the target. This sequence is derived via the routine ADJUST-SEQ. We then decompose the target structure into loops and endow these with a linear order. According to this order we use the routine LOCAL-SEARCH in order to find for each loop a "proper" local solution.

3.1 ADJUST-SEQ

In this section we describe Steps 2 and 3 of the pseudocodes presented in Algorithm 7.1. The routine MAKE-START, see line 8, generates a random sequence, start, which is compatible to the target, with uniform probability.

We then initialize the variable seqmin via the sequence start and set the variable d = + ∞, where d denotes the structure distance between Cross(seqmin) and T.

Given the sequence start, we construct a set of potential "competitors", C, i.e. a set of structures suited as folding targets for start. In Algorithm 7.2 we show how to adjust the start sequence using the routine ADJUST-SEQ. Lines 3 to 36 of Algorithm 7.2, contain a For-loop, executed at most An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i17.gif times. Here the loop-length An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i17.gif is heuristically determined.

For all computer experiments setting the Cross-parameter N = 50, the subroutine executed in the loop-body consists of the following three steps.

Step I. Generating C0(λi) via Cross. Suppose we are in the ith step of the For-loop and are given the sequence λi-1 where λ0 = start. We consider Cross(λi-1, N), i.e. the list of suboptimal structures with respect to λi-1,

equation image

If An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i19.gif, then Inv returns λi-1. Else, in case of An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i20.gif, we set

equation image

Otherwise we do not update seqmin and go directly to Step II.

Step II. The competitors. We introduce a specific procedure that "perturbs" arcs of a given RNA pseudoknot structure, S. Let a be an arc of S and let l(a), r(a) denote the start- and end-point of a. A perturbation of a is a procedure which generates a new arc a', such that

equation image

Clearly, there are nine perturbations of any given arc a (including a itself), see Figure Figure1717.

Figure 17
Perturbations. Nine perturbations of an arc (i, j). Original arcs are drawn dotted, and the arcs incident to red bases are the perturbations.

We proceed by keeping a, replacing the arc a by a nontrivial perturbation or remove a, arriving at a set of ten structures ν(S, a).

Now we use this method in order to generate the set C1(λi-1) by perturbing each arc of each structure An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i23.gif. If An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i24.gif has Ah arcs, An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i25.gif, then

equation image

This construction may result in duplicate, inconsistent or incompatible structures. Here, a structure is inconsistent if there exists at least one position paired with more than one base, and incompatible if there exists at least one arc not compatible with λi-1, see Figures Figures1818 and and19.19. Here compatibility is understood with respect to the Watson-Crick and G-U base pairing rules. Deleting inconsistent and incompatible structures, as well as those identical to the target, we arrive at the set of competitors,

equation image
Figure 18
Inconsistent structures. The dotted arc is perturbed by shifting its end-point. This perturbation leads to a nucleotide establishing two base pairs, which is impossible.
Figure 19
Incompatible structures. We display a perturbation of the dotted arc leading to a structure that is incompatible to the given sequence.

Step III. Mutation. Here we adjust λi-1 with respect to T as well as the set of competitors, C(λi-1) derived in the previous step. Suppose An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i28.gif. Let p(S, w) be the position paired to the position w in the RNA structure S [set membership] C(λi-1), or 0 if position w is unpaired. For instance, in Figure Figure20,20, we have p(T, 1) = 4, p(T, 2) = 0 and p(T, 4) = 1. For each position w of the target T, if there exists a structure Ch(λi-1) [set membership] C(λi-1) such that p(Ch(λi-1, w) ≠ p(T, w) (see positions 5, 6, 9, and 11 in Figure Figure20)20) we modify λi-1 as follows:

Figure 20
Sequence mutation. Suppose the top and middle structures represent the set of competitors and the bottom structure is target. We display λi - 1(top sequence) and its mutation, λi (bottom sequence). Two nucleotides of base pairs not contained ...

1. unpaired position: If p(T, w) = 0, we update An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i29.gif randomly into the nucleotide An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i30.gif, such that for each Ch(λi-1) [set membership] C(λi-1), either p(Ch(λi-1), w) = 0 or An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i31.gif is not compatible with An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i32.gif where v = p(Ch(λi-1), w) < 0, See position 6 in Figure Figure2020.

2. start-point: If p(T, w) <w, set v = p(T, w), We randomly choose a compatible base pair (An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i33.gif) different from (An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i29.gif, An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i32.gif) such that for each Ch(λi-1) [set membership] C(λi-1), either p(Ch(λi-1), w) = 0 or An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i31.gif is not compatible with An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i34.gif, where u = p(Ch(λi-1), w) > 0 is the end-point paired with An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i29.gif in Ch(λi-1) (Figure (Figure20:20: (5, 9). The pair G-C retains the compatibility to (5, 9), but is incompatible to (5, 10)). By Figure Figure2121 we show feasibility of this step.

Figure 21
Sequence adjust. Mutations are always possible, suppose p is paired with q in T and p is paired with q1 in one competitor and q2 in another one. For a fixed nucleotide at p there are at most two scenarios, since a base can pair with at most two different ...

3. end-point: If 0 <p(T, w) <w, then by construction the nucleotide has already been considered in the previous step.

Therefore, updating all the nucleotides of λi-1, we arrive at the new sequence An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i35.gif.

Note that the above mutation steps heuristically decrease the structure distance. However, the resulting sequence is not necessarily incompatible to all competitors. For instance, consider a competitor Ch whose arcs are all contained T. Since λi is compatible with T, λi is compatible with Ch. Since competitors are obtained from suboptimal folds such a scenario may arise.

In practice, this situation represents not a problem, since these type of competitors are likely to be ruled out by virtue of the fact that they have a mfe larger than that of the target structure.

Accordingly we have the following situation, competitors are eliminated due to two, equally important criteria: incompatibility as well as minimum free energy considerations.

If the distance of Cross(λi) to T is less than or equal to dmin + 5, we return to Step I (with λi). Otherwise, we repeat Step III (for at most 5 times) thereby generating An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i36.gif and set An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i37.gif where d(Cross(An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i38.gif), T) is minimal.

The procedure ADJUST-SEQ employs the negative paradigm [16] in order to exclude energetically close conformations. It returns the sequence seqmiddle which is tailored to realize the target structure as mfe-fold.

3.2 DECOMPOSE and LOCAL-SEARCH

In this section we introduce two the routines, DECOMPOSE and LOCAL-SEARCH. The routine DECOMPOSE partitions T into linearly ordered energy independent components, see Figure Figure1313 and Section 2.1. LOCAL-SEARCH constructs iteratively an optimal sequence for T via local solutions, that are optimal to certain substructures of T.

DECOMPOSE: Suppose T is decomposed as follows,

equation image

where the Tw are the loops together with all arcs in the associated stems of the target.

We define a linear order over B as follows: Tw <Th if either

1. Tw is nested in Th, or

2. the start-point of Tw precedes that of Th.

In Figure Figure2222 we display the linear order of the loops of the structure shown in Figure Figure1313.

Figure 22
Ordering. Linear ordering of loops: a1 = [11, 19], b1 = [10, 20], a2 = [7, 37], b2 = [5, 39], a3 = [21, 42], b3 = [20, 44], a4 = [25, 47], b4 = [24, 48], a5 = [7, 47], b5 = [5, 48], a6 = [49, 57], b6 = [48, 59], a7 = [1, 63], b7 = [1, 65].

Next we define the interval

equation image

projecting the loop Tw onto the interval [l(Tw), r(Tw)] and bw = [l', r'] ⊃ aw, being the maximal interval consisting of aw and its adjacent unpaired consecutive nucleotides, see Figure Figure13.13. Given two consecutive loops Tw <Tw + 1, we have two scenarios:

• either bw and bw+1 are adjacent, see b5 and b6 in Figure Figure2222,

• or bw [subset, dbl equals] bw + 1, see b1 and b2 in Figure Figure2222.

Let An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i41.gif, then we have the sequence of intervals a1, b1, c1, ..., am', bm', cm'. If there are no unpaired nucleotides adjacent to aw, then aw = bw and we simply delete all such bw. Thereby we derive the sequence of intervals I1, I2, ..., Im. In Figure Figure2323 we illustrate how to obtain this interval sequence: here the target decomposes into the loops T1, T2 and we have I1 = [3, 5], I2 = [3, 6], I3 = [2, 9], and I4 = [1, 10].

Figure 23
Example of the interval sequence. Loops and their induced sequence of intervals.

LOCAL-SEARCH: Given the sequence of intervals I1, I2, ..., Im. We proceed by performing a local stochastic search on the subsequences An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i42.gif (initialized via seq = seqmiddle and where s|[x, y] = sxsx + 1 ... sy). When we perform the local search on An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i43.gif, only positions that contribute to the distance to the target, see Figure Figure10,10, or positions adjacent to the latter, will be altered. We use the arrays U1, U2 to store the unpaired and paired positions of T. In this process, we allow for mutations that increase the structure distance by five with probability 0.1. The latter parameter is heuristically determined. We iterate this routine until the distance is either zero or some halting criterion is met.

4 Discussion

The main result of this paper is the presentation of the algorithm Inv, freely available at http://www.combinatorics.cn/cbpc/inv.html

Its input is a 3-noncrossing RNA structure T, given in terms of its base pairs (i1, i2) (where i1 <i2). The output of Inv is an RNA sequences s = (s1s2...sn), where sh [set membership] {A, C, G, G} with the property Cross(s) = T, see Figure Figure2424.

Figure 24
The UTR pseudoknot of bovine coronavirus. Its diagram representation and several sequences of its neutral network as constructed by Inv.

The core of Inv is a stochastic local search routine which is based on the fact that each 3-noncrossing RNA structure has a unique loop-decomposition, see Theorem 1 in Section 2.1. Inv generates "optimal" subsequences and eventually arrives at a global solution for T itself. Inv generalizes the existing inverse folding algorithm by considering arbitrary 3-noncrossing canonical pseudoknot structures. Conceptually, Inv differs from INFO-RNA in how the start sequence is being generated and the particulars of the local search itself.

As discussed in the introduction it has to be given an argument as to why the inverse folding of pseudoknot RNA structures works. While folding maps into RNA secondary structures are well understood, the generalization to 3-noncrossing RNA structures is nontrivial. However the combinatorics of RNA pseudoknot structures [27,28,38] implies the existence of large neutral networks, i.e. networks composed by sequences that all fold into a specific pseudoknot structure. Therefore, the fact that it is indeed possible to generate via Inv sequences contained in the neutral networks of targets against competing pseudoknot configurations, see Figure Figure2424 and Figure Figure2525 confirms the predictions of [31].

Figure 25
Pseudoknot PKI. The Pseudoknot PKI of the internal ribosomal entry site (IRES) region [41], its diagram representation and three sequences of its neutral network as constructed by Inv.

An interesting class are the 3-noncrossing nonplanar pseudoknot structures. A nonplanar pseudoknot structure is a 3-noncrossing structure which is not a bi-secondary structure in the sense of Stadler [30]. That is, it cannot be represented by noncrossing arcs using the upper and lower half planes. Since DP-folding paradigms of pseudoknots folding are based on gap-matrices [15], the minimal class of "missed" structures (given the implemented truncations) are exactly these, nonplanar, 3-noncrossing structures. In Figure Figure2626 we showcase a nonplanar RNA pseudoknot structure and 3 sequences of its neutral network, generated by Inv.

Figure 26
Example of nonplanar structure. A nonplanar 3-noncrossing RNA structure together with three sequences realizing them as mfe-structures.

As for the complexity of Inv, the determining factor is the subroutine LOCAL-SEARCH. Suppose that the target is decomposed into m intervals with the length [ell]1, ...., [ell]m. For each interval, we may assume that line 2 of LOCAL-SEARCH runs for fh times, and that line 14 is executed for gh times. Since LOCAL-SEARCH will stop (line 4) if Tstart = T (line 3), the remainder of LOCAL-SEARCH, i.e. lines 7 to 41 run for (fh - 1) times, each such execution having complexity O([ell]h). Therefore we arrive at the complexity

equation image

where c([ell]) denotes the complexity of the Cross. The multiplicities fh and gh depend on various factors, such as start, the random order of the elements of U1, U2 (see Algorithm 7.3) and the probability p. According to [32] the complexity of c([ell]h) is An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i45.gif and accordingly the complexity of Inv is given by

equation image

In Figure Figure2727 we present the average inverse folding time of several natural RNA structures taken from the PKdatabase [39]. These averages are computed via generating 200 sequences of the target's neutral networks. In addition we present in Table Table11 the total time for 100 executions of Inv for an additional set of RNA pseudoknot structures.

Figure 27
Fitting of mean inverse folding time (seconds) over sequence length via 2 cubic spines. For n = 35, ..., 75 we choose a natural pseudoknot structure from the PKdatabase and display the average inverse folding time based on sampling 200 sequences of the ...
Table 1
Inverse folding times for 100 executions of Inv for various RNA pseudoknot structures.

5 Competing interests

The authors declare that they have no competing interests.

6 Authors' contributions

All authors contributed equally to this paper. All authors read and approved the final manuscript.

7 Appendix

7.1 Algorithm 7.1 - INVERSE-FOLD

Input: k-noncrossing target structure T

Output: an RNA sequence seq

Require: k ≤ 3 and T is composed of ":( ) [ ] { }"

Ensure: Cross(seq) = T

1. ▻ Step 1: Validate structure

2. if false = CHECK-STRU(T) then

3.    print incorrect structure

4.    return NIL

5. end if

6.

7. ▻ Step 2: Generate the start sequence

8. start ← MAKE-START(T)

9.

10. ▻ Step 3: Adjust the start sequence

11. seqmiddle ← ADJUST-SEQ(start, T)

12.

13. ▻ Step 4: Decompose T and derive the ordered intervals.

14. Interval array I

15. m ← |I| ▻ I satisfies Im = T

16.

17. ▻ Step 5: Stochastic Local Search

18. seq seqmiddle

19. for all intervals in the array Iw do

20.    l ← start-point(Iw)

21.    r ← end-point(Iw)

22.    s' seq|[l, r] ▻ get sub-sequence

23.    seq|[l, r] LOCAL-SEARCH(s', Iw)

24. end for

25.

26. ▻ Step 6: output

27. if seqmin = Cross(seq) then

28.    return seq

29. else

30.    print Failed!

31.    return NIL

32. end if

7.2 Algorithm 7.2 - ADJUST-SEQ

Input: the original start sequence start

Input: the target structure T

Output: an initialized sequence seqmiddle

1. n ← length of T

2. dmin ← + ∞, seqmin start

3. for i = 1 to An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i47.gif do

4.    ▻ Step I: generate the set C0(λi - 1) via Cross

5.    C0(λi - 1) ← Cross(λi - 1, N)

6.    An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i48.gif

7. if d = 0 then

8.       return λi - 1

9.    else if d <dmin then

10.       dmin d, seqmin λi - 1

11.    end if

12.

13.    ▻ Step II: generate the competitor set C(λi - 1)

14.    C1(λi-1) ← ϕ

15.    for all An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i49.gif [set membership] C1(λi-1) do

16.       for all arc An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i50.gif of An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i49.gifdo

17.          An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i51.gif

18.       end for

19.    end for

20.    C(λi - 1) =

21.    {An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i52.gif is valid}

22.

23.    ▻ Step III: mutation

24.    seq λi - 1

25.    for w = 1 to n do

26.       if [exists] Ch(λi-1) [set membership] C(λi-1) s.t. p(Ch, w) ≠ p(T, w) then

27.          seq[w] ← random nucleotide or pair, s.t. seq [set membership] C[T] and seq [negated set membership] C[Ch(λi-1)]

28.       end if

29.    end for

30.    Tseq Cross(seq)

31.    if d(Tseq, T) <dmin + 5 then

32.       seqmiddle seq

33.    else if Step III run less than 5 times then

34.       goto Step III

35.    end if

36. end for ▻ loop to line 3

37.

38. return seqmiddle

7.3 Algorithm 7.3 - LOCAL-SEARCH

Input:seqmiddle

Input: the target T

Output: seq

Ensure: Cross(seq) = T

1. seq seqmiddle

2. if Cross(seq) = T then

3.    return seq

4. end if

5. decompose T and derive the ordered intervals

6. I ← [I1, I2, ..., Im]

7. for all Iw in I do

8.    ▻ Phase I: Identify positions

9.    An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i53.gif ▻ initialize dmin

10.

11.    derive U1 via An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i54.gif

12.    derive U2 via An external file that holds a picture, illustration, etc.
Object name is 1748-7188-5-27-i54.gif

13.

14.    ▻ Phase II: Test and Update

15.    for all p in U1 do

16.       random T compatible mutate seqp

17.    end for

18.    for all [p, q] in U2 do

19.       random T compatible mutate seqp

20.    end for

21.

22.    E ϕ

23.    for all p [set membership] U1, U2 do

24.       d d(T, Cross(seqp))

25.       if d <dmin then

26.          dmin d, seq seqp

27.          goto Phase I

28.       else if dmin <d <dmin + 5 then

29.          goto Phase I with the probability 0.1

30.       end if

31.       if d = dmin then

32.          E E [union or logical sum] {seq}

33.       end if

34.    end for

35.    seq e0 [set membership] E, where e0 has the lowest mfe in E

36.    if Phase I run less than 10 n times then

37.       goto Phase I

38.    end if

39. end for

40. return seq

8 Acknowledgements

We are grateful to Fenix W.D. Huang for discussions. Special thanks belongs to the two anonymous referee's whose thoughtful comments have greatly helped in deriving an improved version of the paper. This work was supported by the 973 Project, the PCSIRT of the Ministry of Education, the Ministry of Science and Technology, and the National Science Foundation of China.

References

  • Westhof E, Jaeger L. RNA pseudoknots. Curr Opin Struct Biol. 1992;2(3):327–333. doi: 10.1016/0959-440X(92)90221-R. [Cross Ref]
  • Loria A, Pan T. Domain structure of the ribozyme from eubacterial ribonuclease P. RNA. 1996;2:551–563. [PubMed]
  • Staple DW, Butcher SE. Pseudoknots: RNA structures with diverse functions. PLoS Biol. 2005;3(6):e213. doi: 10.1371/journal.pbio.0030213. [PMC free article] [PubMed] [Cross Ref]
  • Konings DA, Gutell RR. A comparison of thermodynamic foldings with comparatively derived structures of 16S and 16S-like rRNAs. RNA. 1995;1:559–574. [PubMed]
  • Tuerk C, MacDougal S, Gold L. RNA pseudoknots that inhibit human immunodeficiency virus type 1 reverse transcriptase. Proc Natl Acad Sci USA. 1992;89(15):6988–6992. doi: 10.1073/pnas.89.15.6988. [PubMed] [Cross Ref]
  • Chamorro A, Manko VS, Denisova TE. New exact solution for the exterior gravitational field of a charged spinning mass. Phys Rev D. 1991;44(10):3147–3151. doi: 10.1103/PhysRevD.44.3147. [PubMed] [Cross Ref]
  • Lyngsø RB, Pedersen CNS. RNA pseudoknot prediction in energy-based models. J Comput Biol. 2000;7(3-4):409–427. doi: 10.1089/106652700750050862. [PubMed] [Cross Ref]
  • Smith TF, Waterman MS. RNA secondary structure: A complete mathematical analysis. Math Biol. 1978;42:257–266.
  • Waterman MS, Smith TF. Rapid dynamic programming methods for RNA secondary structure. Adv Appl Math. 1986;7(4):455–464. doi: 10.1016/0196-8858(86)90025-4. [Cross Ref]
  • Zuker M, Stiegler P. Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information. Nucl Acids Res. 1981;9:133–148. doi: 10.1093/nar/9.1.133. [PMC free article] [PubMed] [Cross Ref]
  • Nussinov B, Jacobson AB. Fast algorithm for predicting the secondary structure of single-stranded RNA. Proc Natl Acad Sci USA. 1980;77(11):6309–6313. doi: 10.1073/pnas.77.11.6309. [PubMed] [Cross Ref]
  • Fresco JR, Alberts BM, Doty P. Some molecular details of the secondary structure of ribonucleic acid. Nature. 1960;188:98–101. doi: 10.1038/188098a0. [PubMed] [Cross Ref]
  • Jun IT, Uhlenbeck OC, Levine MD. Estimation of Secondary Structure in Ribonucleic Acids. Nature. 1971;230(5293):362–367. doi: 10.1038/230362a0. [PubMed] [Cross Ref]
  • Delisi C, Crothers DM. Prediction of RNA secondary structure. Proc Natl Acad Sci USA. 1971;68(11):2682–2685. doi: 10.1073/pnas.68.11.2682. [PubMed] [Cross Ref]
  • Rivas E, Eddy SR. A dynamic programming algorithm for RNA structure prediction including pseudoknots. J Mol Biol. 1999;285(5):2053–2068. doi: 10.1006/jmbi.1998.2436. [PubMed] [Cross Ref]
  • Dirks RM, Lin M, Winfree E, Pierce NA. Paradigms for computational nucleic acid design. Nucleic Acids Res. 2004;32(4):1392–1403. doi: 10.1093/nar/gkh291. [PMC free article] [PubMed] [Cross Ref]
  • Reeder J, Giegerich R. Design, implementation and evaluation of a practical pseudoknot folding algorithm based on thermodynamics. BMC Bioinformatics. 2004;5(104):2053–2068. [PMC free article] [PubMed]
  • Ren J, Rastegari B, Condon A, Hoos H. Hotkonts: Heuristic prediction of RNA secondary structures including pseudoknots. RNA. 2005;15:1494–1504. doi: 10.1261/rna.7284905. [PubMed] [Cross Ref]
  • Huang FWD, Peng WWJ, Reidys CM. Folding 3-noncrossing RNA pseudoknot structures. J Comp Biol. 2009;16(11):1549–75. doi: 10.1089/cmb.2008.0194. [PubMed] [Cross Ref]
  • Borer PN, Dengler B, Tinoco JI, Uhlenbeck OC. Stability of ribonucleic acid doublestranded helices. J Mol Biol. 1974;86(4):843–853. doi: 10.1016/0022-2836(74)90357-X. [PubMed] [Cross Ref]
  • Papanicolaou C, Gouy M, Ninio J. An energy model that predicts the correct folding of both the tRNA and the 5S RNA molecules. Nucleic Acids Res. 1984;12:31–44. doi: 10.1093/nar/12.1Part1.31. [PMC free article] [PubMed] [Cross Ref]
  • Turner DH, Sugimoto N, Freier SM. RNA structure prediction. Ann Rev Biophys Biophys Chem. 1988;17:167–192. doi: 10.1146/annurev.bb.17.060188.001123. [PubMed] [Cross Ref]
  • Walter AE, Turner DH, Kim J, Lyttle MH, Muller P, Mathews DH, Zuker M. Coaxial stacking of helixes enhances binding of oligoribonucleotides and improves predictions of RNA folding. Proc Natl Acad Sci USA. 1994;91(20):9218–9222. doi: 10.1073/pnas.91.20.9218. [PubMed] [Cross Ref]
  • Xia T, SantaLucia JJ, Burkard ME, Kierzek R, Schroeder SJ, Jiao X, Cox C, Turner DH. Thermodynamic parameters for an expanded nearest-neighbor model for formation of RNA duplexes with Watson-Crick base pairs. Biochemistry. 1998;37(42):14719–13735. doi: 10.1021/bi9809425. [PubMed] [Cross Ref]
  • Waterman MS. Combinatorics of RNA hairpins and cloverleaves. Stud Appl Math. 1979;60:91–96.
  • D Kleitman BR. The number of finite topologies. Proc Amer Math Soc. 1970;25:276–282.
  • Jin EY, Qin J, Reidys CM. Combinatorics of RNA structures with pseudoknots. Bull Math Biol. 2008;70:45–67. doi: 10.1007/s11538-007-9240-y. [PubMed] [Cross Ref]
  • Jin EY, Reidys CM. Combinatorial Design of Pseudoknot RNA. Adv Appl Math. 2009;42(2):135–151. doi: 10.1016/j.aam.2008.06.003. [Cross Ref]
  • Chen WYC, Han HSW, Reidys CM. Random k-noncrossing RNA Structures. Proc Natl Acad Sci USA. 2009;106(52):22061–22066. doi: 10.1073/pnas.0907269106. [PubMed] [Cross Ref]
  • Stadler PF. RNA Structures with Pseudo-Knots. Bull Math Biol. 1999;61:437–467. doi: 10.1006/bulm.1998.0085. [PubMed] [Cross Ref]
  • Ma G, Reidys CM. Canonical RNA Pseudoknot Structures. J Comput Biol. 2008;15(10):1257–1273. doi: 10.1089/cmb.2008.0121. [PubMed] [Cross Ref]
  • Huang FWD, Reidys CM. Statistics of canonical RNA pseudoknot structures. J Theor Biol. 2008;253(3):570–578. doi: 10.1016/j.jtbi.2008.04.002. [PubMed] [Cross Ref]
  • Reidys CM, Stadler PF, Schuster P. Generic properties of combinatory maps: neutral networks of RNA secondary structures. Bull Math Biol. 1997;59(2):339–397. doi: 10.1007/BF02462007. [PubMed] [Cross Ref]
  • Reidys CM. Local connectivity of neutral networks. Bull Math Biol. 2008;71(2):265–290. doi: 10.1007/s11538-008-9356-8. [PubMed] [Cross Ref]
  • Hofacker I, Fontana W, Stadler P, Bonhoeffer L, Tacker M, Schuster P. Fast folding and comparison of RNA secondary structures. Chem Month. 1994;125(2):167–188. doi: 10.1007/BF00818163. [Cross Ref]
  • Andronescu M, Fejes AP, Hutter F, Hoos HH, A C. A New Algorithm for RNA Secondary Structure Design. J Mol Biol. 2004;336(2):607–624. doi: 10.1016/j.jmb.2003.12.041. [PubMed] [Cross Ref]
  • Busch A, Backofen R. INFO-RNA--a fast approach to inverse RNA folding. Bioinformatics. 2006;22(15):1823–1831. doi: 10.1093/bioinformatics/btl194. [PubMed] [Cross Ref]
  • Jin EY, Reidys CM. Central and local limit theorems for RNA structures. J Theor Biol. 2008;253(3):547–559. doi: 10.1016/j.jtbi.2007.09.020. [PubMed] [Cross Ref]
  • PseudoBase. http://www.ekevanbatenburg.nl/PKBASE/PKBGETCLS.HTML
  • The pseudoknot structure of the glmS ribozyme pseudoknot P1.1. http://www.ekevanbatenburg.nl/PKBASE/PKB00276.HTML
  • Pseudoknot PKI of the internal ribosomal entry site (IRES) region. http://www.ekevanbatenburg.nl/PKBASE/PKB00221.HTML
  • The pseudoknot of SELEX-isolated inhibitor (ligand 70.28) of HIV-1 reverse transcriptase. http://www.ekevanbatenburg.nl/PKBASE/PKB00066.HTML
  • Pseudoknot PK2 of E.coli tmRNA. http://www.ekevanbatenburg.nl/PKBASE/PKB00050.HTML
  • Pineapple mealybug wilt associated virus - 2. http://www.ekevanbatenburg.nl/PKBASE/PKB00270.HTML

Articles from Algorithms for Molecular Biology : AMB are provided here courtesy of BioMed Central