BMC Evolutionary Biology BioMed Central Methodology article

Background Phylogenies of rapidly evolving pathogens can be difficult to resolve because of the small number of substitutions that accumulate in the short times since divergence. To improve resolution of such phylogenies we propose using insertion and deletion (indel) information in addition to substitution information. We accomplish this through joint estimation of alignment and phylogeny in a Bayesian framework, drawing inference using Markov chain Monte Carlo. Joint estimation of alignment and phylogeny sidesteps biases that stem from conditioning on a single alignment by taking into account the ensemble of near-optimal alignments. Results We introduce a novel Markov chain transition kernel that improves computational efficiency by proposing non-local topology rearrangements and by block sampling alignment and topology parameters. In addition, we extend our previous indel model to increase biological realism by placing indels preferentially on longer branches. We demonstrate the ability of indel information to increase phylogenetic resolution in examples drawn from within-host viral sequence samples. We also demonstrate the importance of taking alignment uncertainty into account when using such information. Finally, we show that codon-based substitution models can significantly affect alignment quality and phylogenetic inference by unrealistically forcing indels to begin and end between codons. Conclusion These results indicate that indel information can improve phylogenetic resolution of recently diverged pathogens and that alignment uncertainty should be considered in such analyses.


Background
The Minimum Evolution (ME) principle is a commonly used principle to estimate phylogenetic trees of a set Γ of n species (taxa) given an n × n symmetric matrix D = {d ij } of evolutionary distances.First introduced by Kidd and Sgaramella-Zonta [1] and subsequently reinterpreted by Rzhetsky and Nei [2,3], the ME principle aims at finding a phylogeny characterized by minimal sum of branch lengths, under the auxiliary criteria that branches have a positive length and the pair-wise distances on the tree are not smaller than the directly observed pair-wise differences.Its biological justification is based on the fact that, when unbiased estimates of the true distances are available, the correct phylogenetic tree has an expected length shorter than any other possible tree [2,3]  where || • || 1 is the -vector norm; v is a vector of the 2n -3 edge lengths; X is a n(n -1)/2 × (2n -3) topological matrix [4] encoding a phylogenetic tree as an unrooted binary tree with the n taxa in Γ as terminal vertices (leaves); is the set of all the topological matrices; finally, f(• , • , •) defines the level of compatibility among the distances in D and the distances induced by the phylogenetic tree edges.Any optimal solution (X*, v*) of problem (1) defines a phylogenetic tree satisfying the minimum evolution principle.A topological matrix X is an Edge-Path incidence matrix of a Tree (EPT) (see [5], and additional files 1 and 2) that encodes a tree as follows: any generic entry x ij,k is set to 1 if the edge k belongs to the path from the leaf i to the leaf j, 0 otherwise.In the rest of the paper we refer to problem (1) as the ME problem.
The distance matrix D of problem (1) is estimated from the dataset, e.g., accordingly to any method described in [6][7][8][9][10][11][12].Condition f(D, X, v) = 0 typically imposes that, for any given EPT matrix X, v minimizes the (weighted) sum of the square values of the differences between the distances in D and the corresponding distances induced by the phylogenetic tree edges [6,13].In particular, under the unweighted least-square (also called Ordinary Least-Squares (OLS)) [2]: where X † is the Moore-Penrose pseudoinverse of X, and D Δ is a vector whose components are obtained by taking row per row the entries of the strictly upper triangular matrix of D.
Others [14] and [15] have suggested the use of a Weighted Least-Squares (WLS) function: where W is a strictly positive definite diagonal matrix whose entries w ij represent weights associated to leaves i and j.Finally, Hasegawa et al. [16] introduced a Generalized Least-Squares (GLS) function in which v is computed using: where C is a strictly positive definite symmetric matrix representing the covariance matrix of D. To avoid the occurrence of negative branch lengths [14,17], problem (1) can be modified as follows: Problem 2. Minimum Evolution under Linear Programming (LP) Unfortunately, both problems (1) and ( 2) are -hard [18].In this context, let us observe that, given Γ, the cardinality of is: where n!! is the double factorial of n.Hence, the number of topological matrices grows exponentially with the number of leaves ( [6], p. 25, and see additional files 1 and 2).
Problem (1) has received great attention from the scientific community such that exact and approximate algorithms to solve it have been developed.Exact algorithms for solving problem (1) are typically based on an exhaustive approach (i.e., enumerating all possible trees X).As an example, PAUP* 4.0 [19] allows exhaustive search for datasets containing up to 12 taxa.A number of heuristics were also developed in the last 20 years.E.g., Rzhetsky and Nei [2,3] (i) start from a Neighbor-Joining (NJ) tree [20,21], (ii) apply a local search generating topologies within a given topological distance (see [2]) from the NJ tree, and (iii) return the best topology found.Kumar [22] further improved the approach as follows: starting from a topology, a leaf l is selects at each step and all possible assignments of l on the topology are tested.Despite that the neighborhood size in Kumar's approach is larger than in Rzhetsky and Nei's algorithm, it requires examining a number of topologies that is, at most, an exponential function of the number of leaves n: (n -1)!/2, and generates solution in a shorter computing time.Finally, Bryant and Waddell [4] implemented programming optimisation and Desper and Gascuel [23] introduced a greedy search that both improved speed and accuracy of the search.
Here, we introduce the Ant Colony Optimization (ACO) algorithm for estimating phylogenies under the minimum evolution principle, and show that ACO has the potential to compete with other widely-used methods.
ACO (see [24,25] for an introduction, and [26,27] for recent reviews) is a widely-used metaheuristic approach for solving hard combinatorial optimization problems.ACO is inspired from the pheromone trail laying and following behavior of real ants.ACO implements indirect communication among simple agents, called (artificial) ants.Communication is mediated by (artificial) pheromone trails implemented as a probabilistic model to which the ants adapt during the algorithm's execution to reflect their search experience.ACO has proven a successful technique for numerous -hard combinatorial optimization problems (see [28]), although no application to the ME phylogeny problem is currently known.Our specific implementation of the ACO algorithm exploits a stochastic version of the Neighbor-Joining (NJ) algorithm [20,21] to explore tree space.

Iterative addition
Given a set Γ of taxa, let us define a partial tree as a m-leaf tree whose leaves are taxa of a subset Γ' ⊂ Γ, with m = |Γ'|.Moreover, given a partial tree, with node set V and edge set E, let us say that we add/insert a leaf i (not yet in Γ') on the edge (r, s) ∈ E (i.e., the edge joining the nodes r, s ∈ V), and generate a new partial tree with node set = V ∪ {i, t} and edge set = E ∪ {(r, t), (t, s), (t, i)}\{(r, s)}.In other words, we add a leaf i on an edge, divide that edge with a new node t, and join the leaf i to t.All algorithms described here build complete phylogenetic trees by iteratively adding one leaf at a time on the edges of a partial tree.

Primal bound
To generate a first upper bound [5] of the ME problem, we adapted the Sequential Addition (SA) greedy algorithm [6].The Sequential Addition algorithm is less prone, than NJ, to generate a systematic local optimum at the end of the search (i.e., starting from "too good" a primal bound may lead to inefficient results [29]).
The pseudo-code of our version of the Sequential Addition algorithm is presented in Figure 1.In the initialization step, we arbitrarily chose a subset Γ' ⊆ Γ of m ≤ n leaves, and we generate as initial m-leaf partial tree, i.e., an optimal solution of the problem (1) when only m leaves are considered.A each iteration, we join the leaf i to all possible leaves already present in Γ', and choose the solution that minimize tree length (we break possible ties randomly), hence, generating a new partial tree and new set Γ' = Γ' ∪ {i}.We iterate the procedure until a tree with n leaves is obtained.Finally, fixing the topology matrix , we determine the optimal edge weights by imposing f (D, , v = 0, and return the length of the tree, i.e., the upper bound on the optimal solution of the ME problem. Unfortunately, the computation complexity of our heuristic is O((2m -5)!! + n(n -m) 2 ).AT each iteration, given a partial tree, i.e., a k-leaf phylogenetic tree of the leaves in Γ' ⊂ Γ with k = |Γ'|, and a leaf i not in Γ', the procedure generates all the different (k + 1)-leaf partial trees that can be obtained by adding the leaf i in each edge of the current partial tree.

The ant colony optimization algorithm
The specific ACO algorithm for the minimum evolution problem (hereafter ACO-ME) that we introduce here (cf.pseudo-code in Figure 2), is a hybrid between the Max-Min Ant System (MMAS) [30,31] and the Approximate Nondeterministic Tree Search (ANTS) [32].Both methods are modifications of the original Ant System approach [33].
The core of the ACO-ME algorithm is the iteration phase, where the ants generate a set of trees.Then, starting from the trees in , a local search is performed until a locally optimal tree is found and compared with the current-best tree.If stopping conditions are met the procedure ends, otherwise the iteration phase is repeated.
Each ant builds a phylogenetic tree by iteratively adding a leaf at a time to a partial tree.Following a relation-learning model [34], the choices performed by an ant about (i) which leaf to insert, and (ii) where to add it on the partial tree are based of a set of parameters {τ ij } called pheromone trails.The values of the pheromone trail parameters {τ ij } represent a stochastic desirability that a leaf i shares a direct common ancestor with a vertex j on a partial tree.The ants generate a new set of trees and the pheromone trail parameters are updated at the end of the main iteration phase.
Let us now consider the algorithm in more details.It uses two identical data structures: s* and s k .The former stores the current-best complete reconstruction (solution) known, whereas the latter stores the best complete reconstruction obtained by the ants during iteration k.The algorithm also uses a variable n a , i.e., the number of artificial ants.How to set the value of n a is discussed in the Parameter settings section.In the initialization phase, s* is first set to the reconstruction obtained by the Sequential Addi High-level pseudo-code for the Sequential Addition heuristic Figure 2 High-level pseudo-code for the Sequential Addition heuristic.tion algorithm and s k is set to null, then the pheromone trail parameters are updated.We implemented the MMAS [30,31] method of pheromone update, where τ min ≤ τ ij ≤ τ max .Here, we set τ min and τ max to 0.0001 and 0.9999, respectively [35].In the initialization phase, the pheromone trail parameters {τ ij } are set to 0.5, i.e., all positions for leaf insertion have the same desirability.
Before describing the iteration phase, let us introduce some definitions.Let be a partial tree with k leaves,

V(
) the set of vertices of , and the set of leaves of .Let us also use the recursive distance definition of [23,36]: if A and B are two non-intersecting subtrees from a tree , then the average distance between A and B is: In the iteration phase, each artificial ant r generates a complete phylogenetic tree using the ConstructCompleteReconstruction(r) procedure, as illustrated in Figure 3: ant r randomly selects four leaves from the set Γ, and builds a partial tree , k = 4, then, ant r (i) chooses, among the leaves not yet inserted in the partial topology, the leaf i defining the smallest distance d ij , j ∈ , and (ii) computes the probability that i has a common ancestor with the vertex j ∈ V ( ) using the formula suggested by ANTS [32]: High-level pseudo-code for the ACO algorithm Figure 3 High-level pseudo-code for the ACO algorithm.
where η ij represents the "heuristic desirability" that leaf i shares a common ancestor with a vertex j of V( ) (whereas τ ij represents the corresponding "stochastic desirability").Finally, α ∈ [01] allows the relative weighting of heuristic and stochastic desirabilities.The heuristic desirability η ij is computed as: where , i.e., the sum of the distances from i to the leaves not yet inserted in the partial tree divided the number of leaves inserted in the partial tree.
Note that η ij , Δ ij , and u i correspond to the quantities used in the Neighbor-Joining algorithm [20,21](see also [23]).Hence, computation of the vector p i = {p ij }, for all i ∈ Γ, can be interpreted as the stochastic application of the Neighbor-Joining algorithm.A possible problem (not observed yet in practice in our analyses) is that η ij can take negative values.Finally, ant r randomly chooses a vertex j on the basis of the probabilities p i , and the leaf i is added to the tree.
At the end of the construction phase, a set of trees is obtained and a 2-OPT local search (with best-improvement and without candidate list [37,38]) is iteratively performed on each tree: two randomly-chosen leaves are swapped and the tree length is evaluated.Swap i is performed on the new tree if swap i-1 generated an improvement, otherwise it is performed on the old tree.To reduce the 2-OPT computational overhead, we perform no more than 10 swappings on each tree in .If the best tree generated by the 2-OPT local search is shorter than the tree in s*, both s* and s k are updated, otherwise only s k is updated.
The pheromone update completes the iteration phase: each entry τ ij is updated following: where where κ ∈ ‫ޒ‬ and ρ, the pheromone evaporation rate, are two tuning constants, s best is one of the tree s* or s k (see below), and the length of s best .When applying equation ( 8), if τ ij is greater than τ max or smaller than τ min , then its value is set to τ max or τ min , respectively.We set to ρ 0.

Parameter settings
We evaluated the performances of the ACO-ME algorithm under different values of the parameter κ(0.1, 0.5, and 1), and different numbers of ants (1 to 10).For each of the 30 possible combinations of these parameters values, we run ACO-ME for 1000 iterations.As suggested elsewhere (see [29]), we do not consider colony sizes larger than 10.
Relative performances are measured using a normalized index as in [39][40][41]: where is the best solution found under parameter value k using dataset j, whereas and are respectively the best and worst solutions found on the instance j using the parameter value k.By definition, performance index values are in the interval [0, 1].The optimal parameter value exhibits the smallest relative performance index (see box-and-whisker plot histograms in Figure (4,5,6).Figures 4, 5, and 6 indicate that, for small, medium, and large datasets, the optimal combinations of number of ants/κ are 7/1, 10/0.5, and 8/0.5, respectively.However, differences of performances are not spectacular among different combinations of parameter values (except that , . x j k ( ) x j min x j max performances are generally very low when a single at is used).

Experimental evaluation
We first used a set of distance matrices generated from real datasets: the dataset "551314.nex" that includes 55 RBCL sequences of 1314 nucleotides each, and the dataset "Zilla500.nex"that includes 500 RBCL sequences of 1428 nucleotides each.These datasets are available at [42].Note that sequences in these datasets were aligned using Clus-talX [43] and columns including gaps were excluded before computing pairwise distances.Second, we generated (i) 10 artificial instances of 20 taxa (also called small instances); (ii) 10 artificial instances of 50 taxa (also called medium instances); and (iii) 10 artificial instances of 100 taxa (also called large instances).Each artificial instance was generated by random sampling of taxa and partial character reshuffling of the Zilla500.nexdata set.More explicitly, after random selection of the 20 or 50 or 100 taxa, we randomly reshuffled characters among taxa, for 50 percents of the aligned columns.As the reshuffling makes the dataset prone to yield undefined pairwise distances [44], we simply used the absolute number of differences between sequence pairs for generating the distance matrix.Edge lengths were computed using the standard OLS because WLS and GLS can potentially lead to inconsistent results et al. [45].All numerical experiments were performed on a workstation Apple 64-bit Power Mac G5 dual processor dual core, with 8 Gb of RAM, and OS X.The ACO-ME source code is written in C/C++ and compiled using IBM XL C/C++ compiler version 6.0.We compared the quality (total length) of trees generated by the ACO-ME algorithm to those obtained using a classical hill-climbing algorithm (implemented in PAUP* 4.0 [19]) after a fixed run time of 1 minute.The starting tree was generated using the Neighbor-Joining algorithm [20,21], and the TBR branch-swapping operator [6] was used for exploring the solution space.PAUP* 4.0 was used with and without the "Steepest Descent" (SD) option.When SD is activated, all possible TBR are tried, and the rearrangement producing the largest decrease in tree length is selected, inducing a computational overhead similar to that of the 2-OPT local search implemented in our ACO-ME algorithm.Each algorithm was run 30 times on each of the two real datasets.Figure 7a  Colony size Normalized Ranking thermore, Figure 7a and 7b indicate that the relative performances of ACO-ME, in comparison to hill climbing, increase with larger datasets.Note that, contrary to our simple implementation of ACO-ME, the implementation of ME in PAUP* 4.0 [19] incorporates procedures [4,23] that greatly speed-up the OLS (reaching a complexity O(n 2 )).We trust that implementation of these procedures in combination with further tuning of the ACO parameters (number of ants, relative weights of the heuristic information and stochastic pheromone parameters, etc) would lead to better performances of the ACO-ME algorithm.Figure 8a and 8b indicate that the relative performances described above are relatively stable trough time, especially for large data sets (at any time during the run, ACO-ME has similar performances than "hill-climbing without SD" and better performances than "hill-climbing with SD").

Conclusion
We introduce here an Ant Colony Optimization algorithm (ACO) for the phylogeny estimation problem under the minimum evolution principle and demonstrate the feasibility of this approach.Although much improve-ment in performances can probably be obtained through (i) modification of the local search phase, (ii) tuning of the ACO parameters (number of ants, relative weights of the heuristic information and stochastic pheromone parameters, etc), and (iii) implementation of speed-up procedures and optimization of the code, the current implementation of our algorithm already demonstrates that the ant colony metaphor can efficiently solve instances of the phylogeny inference problem.Comparison of score vs. running time for hill-climbing with steepest descent (line labeled "1"), hill-climbing without steepest decent (line labeled "2"), and ACO-ME (line labeled "3") on datasets of 55 (a) and 500 (b) taxa.

Problem 1 .
compatible with the distances in D.More formally, the ME principle can be expressed in terms of the following optimization problem: Minimum Evolution under Least Square (LS) Published: 15 November 2007 BMC Evolutionary Biology 2007, 7:228 doi:10.1186/1471-2148-7-228 min

 1
Principle of the ACO-ME algorithmFigure Principle ofthe ACO-ME algorithm.The core iteration includes three main steps: (i) the pheromone update phase during which artificial ants walk on a graph with all possible connections among the n taxa and (n -2) internal nodes, and lay a trail of volatile pheromone on the branches of the starting tree; (ii) the stochastic construction phase during which new trees are built using both the heuristic information of the pairwise distances and the stochastic process guided by the newly-updated pheromone trail matrix (ants follow a given edge with a probability which is a function of the amount of pheromone on that edge); and (iii) the 2-OPT local search phase that corresponds to a local search using taxon swapping.The curved arrow indicates the stochastic jump of an ant from one edge to another.See text for details.

4
and 7b show that ACO-ME performances are intermediate between hill-climbing with SD, and hill-climbing without SD.Fur-Normalized ranking of the ACO algorithm performances with small datasets (20 taxa) and κ = 0.1 (a), κ = 05 (b), and κ = 1 (c) versus colony size n a Figure Normalized ranking of the ACO algorithm performances with small datasets (20 taxa) and κ = 0.1 (a), κ = 05 (b), and κ = 1 (c) versus colony size n a .

6 7
Normalized ranking of the ACO algorithm performances with large datasets (100 taxa) and κ = 0.1 (a), κ = 0.5 (b), and κ = 1 (c) versus colony size n a Figure Normalized ranking of the ACO algorithm performances with large datasets (100 taxa) and κ = 0.1 (a), κ = 0.5 (b), and κ = 1 (c) versus colony size n a .between ACO-ME and hill-climbing (with and without Steepest Descent, SD) after a fixed run time of 1 minute on datasets of 55 (a) and 500 (b) taxa Figure Comparison of performances between ACO-ME and hill-climbing (with and without Steepest Descent, SD) after a fixed run time of 1 minute on datasets of 55 (a) and 500 (b) taxa.A paired Wilcoxon test indicates that ACO-ME performances are significantly better (p-value = 3.92e-2 for 55 taxa dataset, and p-value = 6.821e -4 for 500 taxa dataset) than those of hill-climbing with SD, but significantly worst (p-value = 4.71e -3 for 55 taxa dataset, and p-value = 4.53e -4 for 500 taxa dataset) than those of hill-climbing without SD.Comparison of score vs. running time for hill-climbing with steepest descent (line labeled "1"), hill-climbing without steepest decent (line labeled "2"), and ACO-ME (line labeled "3") on datasets of 55 (a) and 500 (b) taxa Figure8 1, κ to κρ ∈ [10 -2 , 10 -1 ], and α to 0.7.Fine-tuning of these parameters might have a significant impact on search efficiency but such a systematic analysis is out of the scope of a proof-of concept for the use of ACO-ME.Finally, if the objective function does not decrease after 30 iteration, ACO-ME chooses s k as s best instead of s* for the pheromone updating; if the objective function does not decrease after 30 additional iterations, then all {τ ij } are reset to 0.5 and s* is used for pheromone updating.