 Research
 Open Access
 Published:
Evolution of genes neighborhood within reconciled phylogenies: an ensemble approach
BMC Bioinformatics volume 16, Article number: S6 (2015)
Abstract
Context
The reconstruction of evolutionary scenarios for whole genomes in terms of genome rearrangements is a fundamental problem in evolutionary and comparative genomics. The DeCo algorithm, recently introduced by Bérard et al., computes parsimonious evolutionary scenarios for gene adjacencies, from pairs of reconciled gene trees. However, as for many combinatorial optimization algorithms, there can exist many cooptimal, or slightly suboptimal, evolutionary scenarios that deserve to be considered.
Contribution
We extend the DeCo algorithm to sample evolutionary scenarios from the whole solution space under the Boltzmann distribution, and also to compute Boltzmann probabilities for specific ancestral adjacencies.
Results
We apply our algorithms to a dataset of mammalian gene trees and adjacencies, and observe a significant reduction of the number of syntenic conflicts observed in the resulting ancestral gene adjacencies.
Background
The reconstruction of the evolutionary history of genomic characters along a given species tree is a longstanding problem in computational biology. This problem has been well studied for several types of genomic characters, for which efficient algorithms exist to compute parsimonious evolutionary scenarios; classical examples include genes and genomes sequences [1], gene content [2], and gene family evolution [3, 4]. Recently, Bérard et al. [5] extended the corpus of such results to syntenic characters. They introduced the notion of adjacency forest, that models the evolution of gene adjacencies within a phylogeny, motivated by the reconstruction of the architecture of ancestral genomes, and described an efficient dynamic programming (DP) algorithm, called DeCo, to compute parsimonious adjacency evolutionary histories. So far, DeCo is the only existing tractable model that considers the evolution of gene adjacencies within a general phylogenetic framework: other tractable models of genome rearrangements accounting for a given species phylogeny are either limited to singlecopy genes and ignore genespecific events [6], assume restrictions on the gene duplication events, such as considering only wholegenome duplication (see [7] and references there), or require a dated species phylogeny [8].
From a methodological point of view, most existing algorithms to reconstruct evolutionary scenarios along a species tree in a parsimony framework rely on dynamicprogramming along this tree, whose introduction can be traced back to Sankoff in the 1970s (see [9] for a recent retrospective on this topic). Recently, several works considered more general approaches for such parsimony problems that either explore a wider range of values for combinatorial parameters of parsimonious models [10] or consider several alternate histories for a given instance, chosen for example from the set of all possible cooptimal scenarios or from the whole solution space, including suboptimal solutions (see [11–13] for examples of this approach for the gene tree/species tree reconciliation problem).
The present work follows the later approach and extends the DeCo DP scheme toward an exploration of the whole solution space of adjacency histories, under the Boltzmann probability distribution, that assigns a probability to each solution defined in terms of its parsimony score. This principle of exploring the solution space of a combinatorial optimization problem under the Boltzmann probability distribution is sometimes known as the "Boltzmann ensemble approach". It was initially introduced in the context of RNA folding, where the probability of any given conformation at the thermodynamic equilibrium follows a Boltzmann distribution, i.e. a conformation s is observed for a given RNA w with probability ${e}^{{E}_{w,s}/kT}/{\mathcal{Z}}_{w}$, where E_{ w,s } is the freeenergy of conformation s over w, k is the Boltzmann constant, T is the temperature, and ${\mathcal{Z}}_{w}$ is the partition function of w. This latter quantity can be seen as a renormalization factor, and is key in the study of RNA thermodynamics, but its computation involves summing over an exponential number of conformations compatible with the RNA sequence. A major paradigm shift occurred in RNA research when McCaskill [14] showed in 1990 how an efficient algorithm for the partition function could be adapted from a DP scheme for energy minimization through a simple change of algebra. This seminal work also introduced a variant of the insideoutside algorithm [15] for computing basepairing probabilities.
While this Boltzmann ensemble approach has been used for a long time in RNA structure analysis, to the best of our knowledge it is not the case in comparative genomics, where exact probabilistic models have been favoured recently [16, 17]. However, probabilistic models still pose computational challenges for large datasets, and so far a probabilistic model does not exist for gene adjacencies, which motivates our work. In the specific case of the DeCo model, the ability to explore alternative cooptimal or slightly suboptimal solutions is crucial. Indeed, as DeCo models gene adjacencies, each ancestral gene can only be adjacent to at most two other genes, which is not considered in DeCo. However, the initial experiments using DeCo on mammalian gene trees resulted in hundreds of ancestral genes were involved in more than two ancestral gene adjacencies [5]. This raises the question of filtering inferred ancestral adjacencies to reduce the level of syntenic conflict, which can be done on the basis of their Boltzmann probabilities. We reason that some of the erroneouslypredicted adjacencies may result from combinatorial optimization artifacts and that features of a gene adjacency parsimonious evolutionary scenario that are not robust to considering alternative equivalent, or slightly worse, solutions should be considered as dubious.
Methods
Models
A phylogeny is a rooted tree which represents the evolutionary relationships of a set of elements represented by its nodes: internal nodes are ancestors, leaves are extant elements, and edges represent direct descents between parents and children. We consider here three kinds of phylogenies (illustrated in Figure 1): species trees, reconciled gene trees and adjacencies trees/forests. Trees we consider are always rooted. For a tree T and a node x of T , we denote by T (x) the subtree rooted at x. If x is an internal node, we assume it has either one child, denoted by a_{ x }, or two children, denoted by a_{ x } and b_{ x }.
Species trees
A species tree S is a binary tree that describes the evolution of a set of related species, from a common ancestor (the root of the tree), through the mechanism of speciation. For our purpose, species are identified with genomes, and genes are arranged linearly or circularly along chromosomes.
Reconciled gene trees
A reconciled gene tree is also a binary tree that describes the evolution of a set of genes, called a gene family, through the evolutionary mechanisms of speciation, gene duplication and gene loss, within the given species tree S. Therefore, each node of a gene tree G represents a gene loss, an extant gene or an ancestral gene. Ancestral genes are represented by the internal nodes of G, while gene losses and extant genes are represented by the leaves of G.
We denote by s(g) ∈ S the species of a gene g ∈ G, and by e(g) the evolutionary event that leads to the creation of the two children a_{ g } and b_{ g }. If g is an internal node of G, then e(g) is a speciation (denoted by Spec) if the species pair {s(a_{ g }), s(b_{ g })} equals the species pair {a_{s(g)}, b_{s(g)}}, or a gene duplication (GDup) if s(a_{ g }) = s(b_{ g }) = s(g). Finally, if g is a leaf, then e(g) indicates either a gene loss (GLoss) or an extant gene (Extant), in which case e(g) is not an evolutionary event.
Adjacency trees and forests
A gene adjacency is a pair of genes that appears consecutively along a chromosome. An adjacency tree represents the evolution of an ancestral adjacency through the evolutionary events of speciation, gene duplication, gene loss (these events, as described above, occur at the gene level and are modelled in the reconciled gene trees), and adjacency duplication (ADup), adjacency loss (ALoss) and adjacency break (ABreak), that are adjacencyspecific events.

The duplication of an adjacency {g_{1}, g_{2}} follows from the simultaneous duplication of both its genes g_{1} and g_{2} (with s(g_{1}) = s(g_{2}) and e(g_{1}) = e(g_{2}) = GDup), resulting in the creation of two distinct adjacencies each belonging to {a_{g1}, b_{g1}} × {a_{g2}, b_{g2}}.

An adjacency may disappear due to several events, such as the loss of exactly one (gene loss) or both (adjacency loss) of its genes, or a genome rearrangement that breaks the contiguity between the two genes (adjacency break).
Finally, to model the complement of an adjacency break, i.e. the creation of adjacencies through a genome rearrangement, adjacency gain (AGain) events are also considered, and result in the creation of a new adjacency tree. It follows that the evolution of the adjacency between two genes can be described by a forest of adjacency trees, called an adjacency forest. In this forest, each node v belongs to a species denoted by s(v), and is associated to an evolutionary event e(v) ∈ {Spec, GDup, ADup} if g is an internal node, or {Extant, GLoss, ALoss, ABreak} if v is a leaf. Finally, adjacency gain events are associated to the roots of the trees of the adjacency forest. So in the same way that a gene tree G evolves within the species S, an adjacency forest F describing the evolution of the adjacency between two gene families G_{1} and G_{2} evolves within S, G_{1} and G_{2}. We refer the reader to Figure 1 for an illustration.
Parsimony scores and the Boltzmann distribution
When considered in a parsimonious framework, the score of an adjacency forest F is the number of adjacency gains and breaks; other events are not considered as they are the byproducts of evolutionary events already accounted for in the score of the reconciled gene trees G_{1} and G_{2}. We denote by s_{ a }(F) the parsimony score of an adjacency forest F. Let $\mathcal{F}\left({G}_{1},\phantom{\rule{2.36043pt}{0ex}}{G}_{2}\right)$ be the set of all adjacency forests for G_{1} and G_{2}, including both optimal and suboptimal ones, where we assume that at least one extant adjacency is composed of extant genes from G_{1} and G_{2}.
We define the Boltzmann factor of an adjacency forest F as
The partition function associated to two trees G_{1} and G_{2} is obtained as
where kT is an arbitrary constant. The partition function implicitly defines a Boltzmann probability distribution over $\mathcal{F}\left({G}_{1},\phantom{\rule{2.36043pt}{0ex}}{G}_{2}\right)$, where the probability of an adjacency forest F is defined by:
By exponentially favouring adjacency forests with lower parsimony scores, the Boltzmann distribution provides an alternative way to probe the search space, which is heavily influenced by the choice of kT. Indeed, decreasing kT values will skew the Boltzmann distribution towards more parsimonious adjacency forests. Its limiting distributions are uniform over the whole search space (kT → +∞) or over the set of cooptimal forests (kT → 0) (see Figure 2 for an illustration).
A Boltzmann probability distribution on the set of all adjacency forests for a given instance also implies a well defined notion of probability for features of adjacency forests. For example, one can associate a probability to a specific potential ancestral adjacency (i.e. adjacency between two genes from a given ancestral species) as the ratio of the sum of the probabilities of the adjacency forests that contain this adjacency with the partition function.
Algorithms
DeCo, the algorithm described in [5] to compute a parsimonious adjacency forest, is a DP scheme constrained by S, G_{1} and G_{2}. We first present this algorithm, then describe how to extend it into an Boltzmann ensemble algorithm.
The DeCo DP scheme
Let G_{1} and G_{2} be two reconciled gene trees and g_{1} and g_{2} be two nodes, respectively of G_{1} and G_{2}, such that s(g_{1}) = s(g_{2}). The DeCo algorithm computes, for every such pair of nodes g_{1} and g_{2}, two quantities denoted by c_{1}(g_{1}, g_{2}) and c_{0}(g_{1}, g_{2}), that correspond respectively to the most parsimonious score of a parsimonious adjacency forest for the pairs of subtrees G(g_{1}) and G(g_{2}), under the hypothesis of a presence (c_{1}) or absence (c_{0}) of an ancestral adjacency between g_{1} and g_{2}. As usual in DP along a species tree, the score of a parsimonious adjacency forest for G_{1} and G_{2} is given by min(c_{1}(r_{1}, r_{2}), c_{0}(r_{1}, r_{2})) where r_{1} is the root of G_{1} and r_{2} the root of G_{2}.
So, c_{1}(g_{1}, g_{2}) and c_{0}(g_{1}, g_{2}) can be computed as the minimum of a sum of the scores of adjacency gains or breaks and, more importantly, of terms of the form c_{1}(x, y) and c_{0}(x, y) with (x, y) ∈ {g_{1}, a_{g1}, b_{g1}} × {g_{2}, a_{g2}, b_{g2}} − (g_{1}, g_{2}), using the two combinatorial operator min and +.
(Un)ambiguity of the DeCo DP scheme
As defined in [18], the ambiguity of a DP algorithm can be defined as follows: a DP explores a combinatorial solution space (here for DeCo, the space of all possible adjacency forests, including possible suboptimal solutions), that can be explicitly generated by replacing in the equations min by $\u22d3$ (the settheoretic union operator) and + by the Cartesian product × between combinatorial sets. A DP algorithm is then unambiguous if the unions are disjoint, i.e. the sets provided as its arguments do not overlap.
We claim that the DeCo dynamic programming scheme is unambiguous. Indeed, computing c_{1}(g_{1}, g_{2}) and c_{0}(g_{1}, g_{2}) branches on disjoint subcases that each involve a different set of terms c_{1}(x, y) and c_{0}(x, y). The only case that deserves a closer attention is the case where e(g_{1}) = e(g_{2}) = GDup, as a simultaneous duplication can be obtained by two successive duplications. But in this case, the number of AGain events is different (see Figure 3), which ensures the pairwise difference of solutions.
Stochastic backtrack algorithm through algebraic substitutions
As mentioned in [18], any unambiguous dynamic programming scheme can be adapted through algebraic changes to exhaustively generate the set of all adjacency forests, and also compute the corresponding partition function. To that purpose one simply needs to replace the arithmetic operators (min, +) with (Σ, ×), and to exponentiate any atomic cost $C\in \mathbb{R}$ into a (partial) Boltzmann factor e^{C/kT}(see Figure 3).
This precomputation allows us to sample adjacency forests under the Boltzmann distribution, by changing the deterministic backtrack used for maximum parsimony into a stochastic operation. Indeed, assume that the partition function version of the DeCo equation computes c_{1}(g_{1}, g_{2}) (resp. c_{0}(g_{1}, g_{2})) as ${\sum}_{i\in \left[1,{k}_{1}\right]}{t}_{i}$, where the t_{ i } denote the contribution to the partition function of one of the local alternatives within the DP scheme. The latter are typically computed recursively as combinations of atomic adjacency gain/break costs, and recursive terms of the form c_{1}(x, y) and c_{0}(x, y) with (x, y) ∈ {g_{1}, a_{g1}, b_{g1}} × {g_{2}, a_{g2}, b_{g2}} − {(g_{1}, g_{2})}.
Then a (possibly nonparsimonious) random solution can be generated recursively for c_{1}(g_{1}, g_{2}) (resp. c_{0}(g_{1}, g_{2})), by branching on some t_{ i } with probability t_{ i }/c_{1}(g_{1}, g_{2}) (resp. t_{ i }/c_{0}(g_{1}, g_{2})), and proceed recursively on each occurrence of a recursive term within the alternative t_{ i }. The correctness of the algorithm, i.e. the fact that the random process generates each adjacency forests with Boltzmann probability, follows immediately from general considerations on unambiguous DP schemes [18].
The stochastic nature of the backtrack does not affect its worstcase complexity. This Boltzmann sampling algorithm, for an instance composed of two gene trees G_{1} and G_{2} of respective sizes (number of leaves) n_{1} and n_{2}, has time complexity of $\mathcal{O}\left({n}_{1}\times {n}_{2}\right)$ for each backtrack.
Rescaling to avoid numerical overflows
The partition function values $\mathcal{Z}\left({g}_{1},\phantom{\rule{2.36043pt}{0ex}}{g}_{2}\right)$, handled during the computation, typically grow exponentially in the total number of nodes in G_{1} and G_{2}, and may end up overflowing the floating point data type used within the DP tables. Following practice in RNA folding prediction [19], we address this issue by iteratively applying an homogeneous rescaling of these values during the computation, to keep the values found in the DP table asymptotically close to 1, while still allowing for analysis of the Boltzmann distribution.
To that purpose, one introduces a rescaling factor α which is applied, as a multiplicative term, to some of the DP rules. A rescaling is homogeneous for a pair of (sub)trees (G_{1}(g_{1}), G_{2}(g_{2})) (abridged into (g_{1}, g_{2}) from now) when the number of occurrences of α, encountered during the generation of a given solution F, only depends on (g_{1}, g_{2}) and not on specific features of F. Let us denote by κ_{g1, g2 }the number of occurrences of α for (g_{1}, g_{2}), then the rescaled contribution of a given solution F is now ${e}^{\frac{{s}_{a}\left(F\right)}{kT}}\times {\alpha}^{{k}_{g1,\phantom{\rule{2.36043pt}{0ex}}g2}}$, while the rescaled partition function, computed by the modified DP scheme, is given by
A direct execution of the stochastic backtrack algorithm then returns each forest F with probability
In other words, the introduction of the rescaling does not induce any bias in the stochastic sampling, i.e. the sampling still follows a Boltzmann distribution.
On the other hand, α can be used to constrain the values ${\mathcal{Z}}_{\alpha}\left({g}_{1},\phantom{\rule{2.36043pt}{0ex}}{g}_{2}\right)$ to avoid numerical overflows. For instance, setting ${\alpha}^{*}=\mathcal{Z}{\left({G}_{1},\phantom{\rule{2.36043pt}{0ex}}{G}_{2}\right)}^{1/{k}_{g1,\phantom{\rule{2.36043pt}{0ex}}g2}}$ yields ${\mathcal{Z}}_{{\alpha}^{*}}\left({G}_{1},\phantom{\rule{2.36043pt}{0ex}}{G}_{2}\right)=1$. Furthermore, if the rescaling terms are regularly distributed during the execution of the DP scheme, then the intermediate values c_{01}(g_{1}, g_{2}) also typically remain close to 1, thereby avoiding numerical over/underflows. In practice, $\mathcal{Z}\left({g}_{1},\phantom{\rule{2.36043pt}{0ex}}{g}_{2}\right)$ is the end product of the computation, and thus cannot be used to determine a suitable value for α. However, any value that avoids numerical over/underflow can be used, so DeClone accepts as input a prescribed value for α. Note also that α can also be typically inferred from a partial computation, based on the first occurrence of an under/overflow in the DP matrices. To apply these concepts in the context of the DeCo DP scheme, we are left to find an homogeneous rescaling.
Fortunately, we observe that the number of recursive calls ${c}_{01}\left({g}_{1}^{\prime},\phantom{\rule{2.36043pt}{0ex}}{g}_{2}^{\prime}\right)$, where $e\left({g}_{1}^{\prime}\right)$ ≠ GDup and $e\left({g}_{2}^{\prime}\right)$ ≠ GDup, is provably constant within the solutions generated from any call c_{01}(g_{1}, g_{2}). For the sake of simplicity, we assume here that calls of the form c_{01}(g_{1}, g_{2}), where e(g_{1}) = GLoss (resp. e(g_{2}) = GLoss), are expanded into calls c_{01}(g_{1}, a_{g2}) and c_{01}(g_{1}, b_{g2}) (resp. c_{01}(a_{g1}, g_{2}) and c_{01}(b_{g1}, g_{2})), unless g_{2} (resp. g_{1}) is also a leaf. From this observation that can be tediously verified by induction, we adapt the DP scheme as illustrated by Figure 3.
InsideOutside algorithm
While the sampling algorithm described above provides a flexible, easy to implement, approach to analyze the Boltzmann distribution, it only allows for the computation of estimates for properties of interest (for example the occurrence of a specific ancestral adjacency in evolutionary scenarios), whose accuracy may critically depend on the number of samples, the  a priori unknown  variance of the underlying distribution, or other factors. However, whenever the property of interest, in conjunction with the DP scheme, fulfills certain technical conditions [18], it is possible to compute its expectation exactly in polynomial time, by transforming the DP scheme using a variant of the insideoutside algorithm.
More precisely, our objective is to compute the probabilities associated with each of the $\mathcal{O}\left({n}_{1}\times {n}_{2}\right)$ lefthandside (LHS) to righthandside (RHS) transitions in the DP recurrence. Let us denote by l → r an LHS/RHS transition, such that
and by F_{ l→r } the set of forests whose production borrows the l → r transition. The Boltzmann probability of (l → r) is then defined as
Since $\mathcal{Z}\left({g}_{1},\phantom{\rule{2.36043pt}{0ex}}{g}_{2}\right)$ is known, it is sufficient to compute the numerator of the above fraction, i.e. the total Boltzmann factor of the forests ${\mathcal{F}}_{l\to r}$ that feature (l → r). On the other hand, the number of forests in ${\mathcal{F}}_{l\to r}$ typically grows exponentially on n_{1} + n_{2}, so one must find an efficient strategy for computing this summation.
The principle of the insideoutside algorithm [15] is to decompose each of the executions, associated with a forest in ${\mathcal{F}}_{l\to r}$, into: a) an inside part, generated from the recursive calls in the RHS r; and b) an outside part, which denotes the context in which the LHS l appears, i.e. an execution of the DP scheme which features a recursive call to l, and is truncated at that point. Let us remark that the inside and outside parts are independent, i. e. any inside part can be combined with any outside part to form a valid execution of the DP scheme, and the score of the associated forest is simply obtained by summing the scores of its two parts. Thus, the total Boltzmann factor of the forests ${\mathcal{F}}_{l\to r},\phantom{\rule{2.36043pt}{0ex}}l:=\left({x}_{l},{g}_{l},{g}_{l}^{\prime}\right)$, can be decomposed as
where C_{ r } denotes the constant score increment in the RHS, and ${d}_{{x}_{l}}\left({g}_{l},{g}_{l}^{\prime}\right)$ is the outside partition function, i.e. the total Boltzmann factor of all outside parts that are truncated at l. This term can be computed in $\mathcal{O}\left({n}_{1}\times {n}_{2}\right)$ by inverting the DP scheme of Figure 3 in a purely generic, yet quite technical, fashion [18]. To limit the risk of mistakes in the derivation/implementation of DP equations for d_{01}(g_{1}, g_{2}), we implemented an ad hoc parser, based on the inversion principle described by Ponty and Saule [18].
Once the probabilities P (l → r) are known, it is possible to determine the probability of an (ancestral) adjacency (g_{1}, g_{2}) by simply summing over the probabilities of transitions that infer such an adjacency, i. e. that feature a recursive call of the form c1(g_{1}, g_{2}) within their RHS. Iterating this over all (g_{1}, g_{2}) pairs, one obtains an adjacency matrix, as shown in Figure 2.
Results and discussion
Data
We reanalyzed a dataset described in [5] composed of 5, 039 reconciled gene trees and 50, 389 extant gene adjacencies, forming 6, 074 DeCo instances, with genes taken from 36 mammalian genomes from the Ensembl database in 2012. In [5], these data were analyzed using the DeCo algorithm that computed a single parsimonious adjacency forest per instance. All together, these adjacency forests defined 112, 188 (resp. 96, 482) ancestral and extant genes (resp. adjacencies), where, by "ancestral adjacency", we mean adjacency that involves two genes g_{1} and g_{2} whose descendants in their respective gene trees satisfy that they do not belong to the same species s(g_{1}) (equal to s(g_{2})), i.e. g_{1} and g_{2} are prespeciation genes, that were not duplicated within their species (this choice is motivated by the fact that the reconstruction of ancestral genomes considers prespeciation genomes. More important, we can observe 5, 817 ancestral genes participating to three or more ancestral adjacencies, which represent a significant level of syntenic conflict (close to 5%), as a gene can only be adjacent to at most two neighboring genes along a chromosome.
DeCo scores, solution space
Unlike reconciled gene trees, whose mutation cost can be high, most adjacency forests have a relatively low cost, with only 32 instances leading to forest of score 5 or above, while the average number of parsimonious syntenic events (adjacency gain and break) is 1.25. This illustrates the fact that syntenic events, that are due to genome rearrangements, are rare evolutionary events, which suggests that parsimony is a relevant criterion for such characters, and that robustness of syntenic characters with respect to the whole solution space should be assessed in terms of optimal or slightly suboptimal evolutionary scenarios.
Boltzmann sampling and exact Boltzmann probabilities
For each instance, we sampled 1, 000 adjacency forests under the Boltzmann distribution, for three values of kT, 0.001, 0.1, 0.5, and recorded the frequency of all observed ancestral adjacencies. Then for the same values of kT, we computed the exact Boltzmann probability of all potential ancestral adjacencies using the insideoutside algorithm. The result observed were very similar whether sampling or exact probabilities were considered. However, the time required to compute exact Boltzmann probabilities is polynomial, so the exact Boltzmann approach based on the insideoutside algorithm should naturally be favoured in applications. In consequence, we discuss only the case of exact Boltzmann probabilities below.
The main difference between the three values of kT is that, with kT = 0.5, nonoptimal adjacency forests have a higher Boltzmann probability in the Boltzmann distribution, while kT = 0.1 skews the distribution toward optimal adjacency forests and slightly suboptimal ones, and kT = 0.01 ensures that the probability of suboptimal adjacency forests is extremely low and almost does not contribute to the partition function. We then looked at the numbers of ancestral adjacencies, genes and syntenic conflicts from ancestral adjacencies in terms of Boltzmann probability. Table 1 below summarizes the obtained results.
The difference observed between the results with different values of kT supports that parsimony is an appropriate criterion for looking at gene adjacency evolution. Indeed, in the results obtained with kT = 0.5, that gives a higher probability to nonoptimal adjacency forests, it appears that the number of conserved ancestral adjacencies drops sharply after probability 0.6, showing that very few ancestral adjacencies appear with high probability. However, with kT = 0.1 and kT = 0.01, by taking a high probability threshold (starting at a threshold of 0.6), we reduce significantly the number of syntenic conflicts while maintaining a relatively similar number of ancestral genes than the experiments described in [5]; this observation illustrates the potential of the ensemble approach compared to the classical dynamic approach that relies on a single arbitrary optimal solution. Next, the experiment with kT = 0.01 that considers only cooptimal scenarios (the probability of nonoptimal scenarios falls under the numerical precision) shows that, despite conserving only ancestral adjacencies with maximal support in terms of Boltzmann probability, a significant number of syntenic conflicts remains. We conjecture that this is due to errors in the considered reconciled gene trees, and it would be interesting to see if the information about highly supported conflicting adjacencies can be used to correct reconciled gene tree.
Conclusions
The main contribution of our work is an extension of the DeCo dynamic programming scheme to consider adjacency forests in a probabilistic framework, under the Boltzmann distribution. The application of our algorithms on a mammalian genes dataset, together with a simple thresholdbased, approach to filter ancestral adjacencies, proved to be effective to reduce significantly the number of syntenic conflicts, illustrating the interest of the ensemble approach. This preliminary work raises several questions and can be extended along several lines. Among them, we can cite two of immediate interest. First, given the Boltzmann probabilities of the adjacency gains and breaks associated to ancestral adjacencies, we could use them to compute a Maximum Expected Accuracy adjacency forest, which is a parsimonious adjacency forest in a scoring model where each event is weighted by Boltzmann probability (see [20] for an example of this approach for RNA secondary structures). This would provide a unique evolutionary scenario per instance. Next, we considered here an evolutionary model based on speciation, duplication and loss. A natural extension would be to include the event of lateral gene transfer in the model. Efficient reconciliation algorithms exist for several variants of this model [3, 4], together with an extension of DeCo, called DeCoLT [21]. DeCoLT is also based on dynamic programming, and it is likely that the techniques we developed in the present work also apply to this algorithm.
References
 1.
Liberles DA: Ancestral Sequence Reconstruction. 2007, Oxford University Press, Oxford, UK
 2.
Csürös M: Ancestral reconstruction by asymmetric Wagner parsimony over continuous characters and squared parsimony over distributions. RECOMBCG. Lecture Notes in Computer Science. 2008, Springer, Berlin, Germany, 5267: 7286. doi:10.1007/9783540879893_6
 3.
Bansal MS, Alm EJ, Kellis M: Efficient algorithms for the reconciliation problem with gene duplication, horizontal transfer and loss. Bioinformatics. 2012, 28 (12): 283291. doi:10.1093/bioinformatics/bts225
 4.
Doyon JP, Scornavacca C, Gorbunov KY, Szöllosi GJ, Ranwez V, Berry V: An efficient algorithm for gene/species trees parsimonious reconciliation with losses, duplications and transfers. RECOMBCG. Lecture Notes in Computer Science. 2010, Springer, Berlin, Germany, 6398: 93108. doi:10.1007/9783642161810_9
 5.
Bérard S, Gallien C, Boussau B, Szöllosi GJ, Daubin V, Tannier E: Evolution of gene neighborhoods within reconciled phylogenies. Bioinformatics. 2012, 28 (18): 382388. doi:10.1093/bioinformatics/bts374
 6.
Biller P, Feijão P, Meidanis J: Rearrangementbased phylogeny using the singlecutorjoin operation. IEEE/ACM Trans. Comput. Biology Bioinform. 2013, 10 (1): 122134. doi:10.1109/TCBB.2012.168
 7.
Gagnon Y, Blanchette M, ElMabrouk N: A flexible ancestral genome reconstruction method based on gapped adjacencies. BMC Bioinformatics. 2012, 13 (S19): 4doi:10.1186/1471210513S19S4
 8.
Ma J, Ratan A, Raney BJ, Suh BB, Zhang L, Miller W, Haussler D: DUPCAR: reconstructing contiguous ancestral regions with duplications. Journal of Computational Biology. 2008, 15 (8): 10071027. doi:10.1089/cmb.2008.0069
 9.
Csürös M: How to infer ancestral genome features by parsimony: Dynamic programming over an evolutionary tree. Models and Algorithms for Genome Evolution. 2013, Springer, Berlin, Germany, 2945. doi:10.1007/9781447152989_3
 10.
LibeskindHadas R, Wu YC, Bansal MS, Kellis M: Paretooptimal phylogenetic tree reconciliation. Bioinformatics. 2014, 30 (12): 8795. doi:10.1093/bioinformatics/btu289
 11.
Bansal MS, Alm EJ, Kellis M: Reconciliation revisited: Handling multiple optima when reconciling with duplication, transfer, and loss. Journal of Computational Biology. 2013, 20 (10): 738754. doi:10.1089/cmb.2013.0073
 12.
Scornavacca C, Paprotny W, Berry V, Ranwez V: Representing a set of reconciliations in a compact way. J. Bioinformatics and Computational Biology. 2013, 11 (2): doi:10.1142/S0219720012500254
 13.
Doyon JP, Hamel S, Chauve C: An efficient method for exploring the space of gene tree/species tree reconciliations in a probabilistic framework. IEEE/ACM Trans. Comput. Biology Bioinform. 2012, 9 (1): 2639. doi:10.1109/TCBB.2011.64
 14.
McCaskill JS: The equilibrium partition function and base pair binding probabilities for RNA secondary structure. Biopolymers. 1990, 29 (67): 11051119. doi:10.1002/bip.360290621
 15.
Baker JK: Trainable grammars for speech recognition. The Journal of the Acoustical Society of America. 1979, 65 (S1): 132132. doi:10.1121/1.2017061
 16.
Arvestad L, Lagergren J, Sennblad B: The gene evolution model and computing its associated probabilities. J. ACM. 2009, 56 (2): doi:10.1145/1502793.1502796
 17.
Mahmudi O, Sjöstrand J, Sennblad B, Lagergren J: Genomewide probabilistic reconciliation analysis across vertebrates. BMC Bioinformatics. 2013, 14 (S15): 10doi:10.1186/1471210514S15S10
 18.
Ponty Y, Saule C: A combinatorial framework for designing (pseudoknotted) RNA algorithms. Algorithms in Bioinformatics (Proceedings of WABI'11). Lecture Notes in Computer Science. Edited by: Przytycka, T., Sagot, M.F. 2011, Springer, Berlin Heidelberg, Germany, 6833: 250269. doi:10.1007/9783642230387_22
 19.
Hofacker IL, Fontana W, Stadler PF, Bonhoeffer LS, Tacker M, Schuster P: Fast folding and comparison of RNA secondary structures. Monatshefte für Chemie/Chemical Monthly. 1994, 125 (2): 167188. doi:10.1007/BF00818163
 20.
Clote P, Lou F, Lorenz WA: Maximum expected accuracy structural neighbors of an RNA secondary structure. BMC Bioinformatics. 2012, 13 (S5): 6doi:10.1186/1471210513S5S6
 21.
Patterson M, Szöllosi GJ, Daubin V, Tannier E: Lateral gene transfer, rearrangement, reconciliation. BMC Bioinformatics. 2013, 14 (S15): 4doi:10.1186/1471210514S15S4
Acknowledgements
J.P.P.Z. visit to Simon Fraser University was funded by the São Paulo Research Foundation (FAPESP).
Declarations
The publication charges this article were funded by the Simon Fraser University (SFU) Open Access fund.
This article has been published as part of BMC Bioinformatics Volume 16 Supplement 19, 2015: Brazilian Symposium on Bioinformatics 2014. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/16/S19
Author information
Affiliations
Corresponding author
Correspondence to Cedric Chauve.
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
All authors participated in all aspects of the study.
Cedric Chauve, Yann Ponty and João Paulo Pereira Zanetti contributed equally to this work.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Chauve, C., Ponty, Y. & Zanetti, J.P.P. Evolution of genes neighborhood within reconciled phylogenies: an ensemble approach. BMC Bioinformatics 16, S6 (2015). https://doi.org/10.1186/1471210516S19S6
Published:
Keywords
 Synteny
 Parsimony
 Gene adjacency
 Dynamic programming
 Ensemble approach