Abstract
Background
Understanding the evolutionary relationships among species based on their genetic information is one of the primary objectives in phylogenetic analysis. Reconstructing phylogenies for large data sets is still a challenging task in Bioinformatics.
Results
We propose a new distancebased clustering method, the shortest triplet clustering algorithm (STC), to reconstruct phylogenies. The main idea is the introduction of a natural definition of socalled krepresentative sets. Based on krepresentative sets, shortest triplets are reconstructed and serve as building blocks for the STC algorithm to agglomerate sequences for tree reconstruction in O(n^{2}) time for n sequences.
Simulations show that STC gives better topological accuracy than other tested methods that also build a first starting tree. STC appears as a very good method to start the tree reconstruction. However, all tested methods give similar results if balanced nearest neighbor interchange (BNNI) is applied as a postprocessing step. BNNI leads to an improvement in all instances. The program is available at http://www.bi.uniduesseldorf.de/software/stc/ webcite.
Conclusion
The results demonstrate that the new approach efficiently reconstructs phylogenies for large data sets. We found that BNNI boosts the topological accuracy of all methods including STC, therefore, one should use BNNI as a postprocessing step to get better topological accuracy.
Background
Reconstructing the evolutionary relationships among species based on their genetic information is one of the primary objectives in phylogenetic analysis. In recent years, numerous heuristics to reconstruct phylogenies for large data sets have been proposed [111]. In addition, parallel treereconstruction programs have been implemented [1215].
To date, distancebased methods introduced by CavalliSforza and Edwards [16] and by Fitch and Margoliash [17] appear most appropriate to reconstruct phylogenies based on thousands of sequences. These methods are a compromise between computational speed and topological accuracy [1,3,57] and run typically in O(n^{3}) time for n sequences [1,3,5] or in O(n^{2}) for recently suggested approaches [6,7]. Clustering algorithms form a major class of distancebased methods [18]. They do not have an explicit objective function that needs to be optimized. They rather group sequences (or taxa) iteratively to reconstruct a distancebased phylogenetic tree. UPGMA is a popular method to infer phylogenies with the constraint that a molecular clock is imposed on the evolutionary process. Other clustering approaches have been proposed to relax the molecular clock assumption [1,3,5,1921].
An attempt to boost the accuracy and to reduce the computational burden is the introduction of krepresentative set concepts [10,11]. krepresentative sets consist of at most k elements but retain the most important information from whole sets. In this paper, we extend our original approach [10] by introducing a more natural krepresentative set concept. In a nutshell, representative sets are regarded as components to construct shortest triplets, each of which comprises three closely related sequences from three krepresentative sets. The collection of shortest triplets serves as building block for a new distancebased clustering method called shortest triplet clustering algorithm (STC).
Results
Simulations were run on a PC cluster with 16 nodes. Each node has two 1.8 GHz processors and 2 GB RAM. SeqGen [22] was used to evolve sequences along trees using the Kimura twoparameter model [23] with a transition/transversion ratio of 2.0. We generated 100 simulated data sets of 500 sequences each with sequence lengths 500, 1000 and 2000 nucleotides (nt), respectively. As one model tree, we used the rbcl gene tree with diameter 0.36 substitutions per site as inferred from an alignment of 500 rbclgenes [10]. We call this the rbclsimulation.
In a second experiment, the socalled large simulation, tree topologies were drawn from the YuleHarding distribution [24], and edge lengths were drawn from an exponential distribution and subsequently rescaled such that the mean diameter of the tree was either 0.1, 0.5, 1.0, or 1.5. For each value of the diameter we generated 100 trees with 1000 sequences and 100 trees with 5000 sequences. Thus, a total of 800 trees were used.
Finally, we tested the accuracy and runtime of the STC and compared it with six other commonly used distancebased methods. More specifically, we investigate the performance of the NeighborJoining method (NJ) [1] implemented in PAUP* 4.0 [25], BIONJ [3], Weighbor 1.2 [5], Harmony Greedy Triplet and Four Point Condition (HGT/FP) [7] as well as Greedy Minimum Evolution (GME) and Balanced Minimum Evolution (BME) [6]. Unfortunately, no distancebased program is available for the disccovering method [4]. All methods were combined with DNADIST version 3.5 [26] and pairwise distances were corrected for multiple hits according to the model used in the simulation. Moreover, we examined the performance of all methods when the balanced nearest neighbor interchange (BNNI) [6] is used as a postprocessing step.
Further, to illustrate the performance of STC we reanalyzed the 96taxon alignments of sequence length 500 nt, that were analyzed in [6] and available at http://www.lirmm.fr/~guindon/simul/ webcite. The 6000 trees were split into three groups called "slow" (0.2 substitutions per site), "moderate" (0.4 substitutions per site) and "fast" (1.0 substitutions per site). We call this the reanalyzed simulation.
The accuracy of a tree reconstruction method for a simulated data set is measured by the Robinson and Foulds (RF) distance [27] between the inferred tree and the model tree used to generate the data set. The RF distance between two trees is the number of bipartitions present in one of the two trees but not the other, divided by the number of possible bipartitions. Thus, the smaller the RF distance between two trees the closer are their topologies. In other words, the smaller the RF distance is between the inferred tree and the model tree the higher is the topological accuracy of the tree reconstruction method.
In the following we discuss the results of the rbclsimulation, and the large simulation and the reanalyzed simulation.
rbclsimulation
Table 1 shows that the STC outperforms all other methods analyzed in terms of topologlcal accuracy. For instance, the RF distance between the STCtree and the model tree is on average 0.177 (with respect to the sequence length of 500 nt) and better than NJ (0.190), slightly better than the second best method BME (0.184) and much better than HGT/FP (0.512). Table 1 also demonstrates that all tested methods including STC give higher topologlcal accuracy when the sequence length is increased. However, Table 2 shows that other methods in combination with BNNI outperform STC without BNNI. The combination of STC and BNNI shows similar performance as the combinations of NJ (BIONJ, Weighbor) and BNNI and, slightly better results than the combination of GME (HGT/FP) and BNNI.
Table 1. The average Robinson and Foulds distance of 100 simulated data sets of 500 sequences each with sequence lengths 500, 1000 and 2000 nt (rbcl simulation). Methods are used without BNNI.
Table 2. The average Robinson and Foulds distance of 100 simulated data sets of 500 sequences each with sequence lengths 500, 1000 and 2000 nt (rbcl simulation). Methods are used with BNNI.
Large simulation
Due to the increase in runtime, Weighbor could not be tested. Table 3 and 4 show that STC gives better results than the other methods independent of the diameter. All methods display a decrease in accuracy when the number of sequences changes from 1000 to 5000. As shown in Table 5 and 6, BNNI boosts the accuracy of all methods including STC. All methods give similar results when being used together with BNNI.
Table 3. The average Robinson and Foulds distance of 100 simulated data sets of 1000 taxa for each tree diameter 0.1, 0.5, 1.0 and 1.5 and with sequence length 1000 nt (large simulation). Methods are used without BNNI.
Table 4. The average Robinson and Foulds distance of 100 data sets of 5000 taxa for each tree diameter 0.1, 0.5, 1.0 and 1.5 and with sequence length 1000 nt (large simulation). Methods are used without BNNI.
Table 5. The average Robinson and Foulds distance of 100 simulated data sets of 1000 taxa for each tree diameter 0.1, 0.5, 1.0 and 1.5 and with sequence length 1000 nt (large simulation). Methods are used with BNNI.
Table 6. The average Robinson and Foulds distance of 100 data sets of 5000 taxa for each tree diameter 0.1, 0.5, 1.0 and 1.5 and with sequence length 1000 nt (large simulation). Methods are used with BNNI.
Reanalyzed simulation
Except for STC, the accuracies for the other methods displayed in Table 7 and 8 were taken from [6]. Table 7 shows that STC outperforms the other methods in terms of topological accuracy with the exception that Weighbor is slightly better than STC with respect to the slow simulation group. If BNNI is applied, all methods exhibit an almost identical performance (see Table 8).
Table 7. The average RF distance of the 96taxon alignments of sequence length 500 nt, that were analyzed in [6]. The 6000 trees were split into three groups called "slow" (0.2 substitutions per site), "moderate" (0.4 substitutions per site) and "fast" (1.0 substitutions per site). Except for STC, the accuracies for the other methods were taken from [6]. Methods are used without BNNI.
Table 8. The average RF distance of the 96taxon alignments of sequence length 500 nt, that were analyzed in [6]. The 6000 trees were split into three groups called "slow" (0.2 substitutions per site), "moderate" (0.4 substitutions per site) and "fast" (1.0 substitutions per site). Except for STC, the accuracies for the other methods were taken from [6]. Methods are used with BNNI.
Another look at the performance
Instead of looking at the average RF distance, we suggest to take a closer look at the simulated data. For each simulated data set, that is subjected to the STC and six other tree reconstruction methods mentioned above, we compute the RF distance between the reconstructed tree and the model tree for all methods. Figure 1 illustrates the results for the large simulation when comparing STC with NJ (left column) and STC with the second best method BME (right column). In each diagram specified by the number of taxa and reconstruction methods, 400 points are displayed, that resulted from 100 simulations for each of the treediameters (0.1, 0.5, 1.0 and 1.5). Although four treediameters were studied only two clouds of points are discernible, where the cloud in the northeast corner of each diagram represents the simulations with the treediameter 0.1. The remaining 300 points gather in the southwest cloud because the RFdistances from trees with diameter 0.5, 1.0, 1.5 are not substantially different from each other (see Table 3 and 4). More precisely, the horizontal and vertical axes indicate the RF distances of STC and NJ (or BME), respectively. Each point in the graph presents the RF distance for a simulated data set. Points above the dotted line are examples where the RF distance of the STCtree is less than the RF distance of the NJtree or BMEtree. Thus, the STC gives higher topological accuracy than NJ or BME with respect to the simulated data set. For example, Figure 1a illustrates the comparison between STC and NJ with respect to 1000 taxa data sets. 379 out of 400 points are above the diagonal, thus, STC gives better results than NJ in about 95% of the simulations. For the remaining 21 alignments (points), two methods showed the same RF distance. Finally, we found 19 points below the diagonal in which case NJ outperforms STC. For the large simulation (5000 taxa), NJ is worse than STC in all cases. However, the second best method BME is better than STC in 11% and 5% of the cases with respect to 1000 and 5000 sequence data sets.
Figure 1. The comparisons of topological accuracy between STC, NJ and BME for the large simulation. Each point in the graph presents the Robinson and Foulds (RF) distance for a simulated data set. Points above the dotted line are examples where the RF distance of the STCtree is less than the RF distance of the NJtree or BMEtree. Thus, the STC gives higher topological accuracy than NJ or BME with respect to the simulated data set.
Figure 2 shows the same analysis for the rbcl simulation. It shows that with increasing sequence length the cloud of points moves towards zero. From Figure 2 we learn that in some instances NJ (or BME) performs better (with regard to the RF distance) than STC, i.e. 20%, 12%, 8% (or 34%, 17%, 14%) of the simulations for sequence lengths 500, 1000 and 2000 nt, respectively.
Figure 2. The comparisons of topological accuracy between STC, NJ and BME for the rbcl simulation. Each point in the graph presents the Robinson and Foulds (RF) distance for a simulated data set. Points above the dotted line are examples where the RF distance of the STCtree is less than the RF distance of the NJtree or BMEtree. Thus, the STC gives higher topological accuracy than NJ or BME with respect to the simulated data set.
Similar results hold for the other methods. These results are summarized in Table 9 where we show the percentage of simulations in which STC is at least as good as the other methods.
Table 9. The percentage of cases where STC is at least as good as other tested methods in terms of RF distance. The number in parentheses is the percentage of cases where STC is equally good as other tested methods. Methods are used without BNNI.
Again, if BNNI is applied we observe that no substantial difference among the various approaches. The accuracy of the methods is mostly determined by BNNI (see Table 10).
Table 10. The percentage of cases where STC is better than other tested methods in terms of RF distance. The number in parentheses is the percentage of cases where STC is worse than other tested methods. Methods are used with BNNI.
Conclusion
We are presenting krepresentative sets which allow us to design a fast and accurate method to reconstruct phylogenies from large data sets with 1000 or more taxa. Simulations show that STC gives better results than other tested methods in terms of topological accuracy. However, if BNNI is introduced as a subsequent optimization step, the differences in the performance disappear. All methods show more or less the same accuracy. Thus, one should apply BNNI to improve the topological accuracy.
The time to reconstruct a tree of up to 1000 sequences is not really an issue for all tested distancebased methods, with the exception of Weighbor. Weighbor needed about 19 minutes to reconstruct a tree with 500 sequences, thus it is only applicable to data sets with up to some hundred sequences. For data sets with up to 1000 sequences, the remaining methods needed less than one minute to output a tree, thus the difference between methods in terms of runtime is not significant. For data sets with 5000 sequences, STC (GME, HGT/FP or BME) with BNNI took about 2.0 (2.5, 3.0 or 3.5) minutes to reconstruct a tree. NJ (BIONJ) with BNNI were slower and consumed approximately six minutes to output a tree. In short, the combination of STC and BNNI efficiently reconstruct trees for large data sets in both terms of topological accuracy and runtime.
Finally, we did not systematically evaluate the impact of the number of representatives k. We present some preliminary results for k = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90 and 100. Figure 3 shows that the RF distance of STC decreases when k grows from 1 to 5. This proves our intuition that a too small number of triplets leads to an inaccurate estimate of path lengths and edge lengths. When k ranges from 5 to 10, the RF distance remains more or less unchanged. For k ≥ 10, the RF distance increases steadily indicating a loss of accuracy. The decrease in accuracy is explained by the inclusion of triplets with large distances which include noise and disturb the reconstruction. Thus, we chose k = 5 as a good compromise between the accuracy and computational complexity for all data sets. That is, the practical complexity of the STC algorithm is only O(n^{2}).
Figure 3. The impact of the number of representatives k. The RF distance of STC decreases when k grows from 1 to 5. When k ranges from 5 to 10, the RF distance remains more or less unchanged. For k ≥ 10, the RF distance increases steadily indicating a loss of accuracy.
Methods
In this section we introduce a new clustering algorithm to reconstruct phylogenies based on distance matrices.
Additive distances
Let S = {s_{1}, s_{2},..., s_{n}} be a set of n objects (typically contemporary sequences/taxa), let D = D(uv) be a distance matrix where D(uv) is the distance between two objects u and v.
Definition 1
The distance matrix D is additive if and only if it satisfies the fourpoint condition [28]: for any quartet {u, v, w, x},
D(uv) + D(wx) ≤ max{D(uw) + D(vx), D(ux) + D(vw)}.
In this case, the objects s ∈ S are related by a tree T = (V, E) where V is the set of vertices such that S ⊂ V and E = {{v_{1}, v_{2}}v_{1}, v_{2 }∈ V} is the set of edges. A vertex with one adjacent edge is called a leaf, all other vertices are called internal nodes. We let L ⊂ V be the leaf set of the tree T. Note that we typically require L ⊆ S in the phylogenetic setting.
If D is additive, then there exists a map and a length function _{} such that
for all u, v ∈ S where p(φ (u), φ(v)) is the unique path connecting φ(u) and φ(v) in T and _{} denotes the distances between vertices in T (cf. [29]). ℓ(e) is called edge length of the edge e. To avoid unnecessary complication, we consider only onetoone maps from S on the leaf set L of T. If D is additive, the reconstruction of tree T and ℓ is trivial. If D is not additive, methods are available that try to fit a tree T to D with respect to an objective function (cf. [30]). Thus, in the following we consider arbitrary distance matrices and we want to reconstruct a tree together with a length function .
Estimating edge lengths using triplets
We consider a subset X of S, then induces a map on a subtree of T such that the relationships of objects in X are displayed by the subtree with leaf set φ(X). The complement S_{0}(X) = S  X we will call the unclassified object set, because the relationships of objects in S_{0}(X) to X is not known from the subtree. Note that we will use S_{0 }instead of S_{0}(X) if X is clear from the context.
Let denote T_{r }= (V_{r}, E_{r}) a rooted tree with root r and leaf set L_{r}, and let S_{r }be a subset of S such that φ(S_{r}) = L_{r}. For convenience, we use S_{r }and L_{r }interchangeably.
Now, we consider the most simple edge length estimation problem. That is, we would like to estimate the edge lengths for a triplet tree {a, b, c} with distance matrix D (see Figure 4a). Edge lengths are estimated as follows
Figure 4. On the left, estimation of edge lengths ℓ(arabc), ℓ(brabc) and ℓ(crabc) of the triplet tree {a, b, c}. On the right, estimation of path length ℓ(s_{0}rs_{0}s_{1}s_{2}) and edge lengths ℓ(r_{1}rs_{0}s_{1}s_{2}), ℓ(r_{2}rs_{0}s_{1}s_{2}) based on the triplet tree {s_{0}, s_{1}, s_{2}}.
Now consider a rooted T_{r }with the inferred treelike metric . The rooted tree T_{r }consists of two rooted subtrees and (see Figure 4b). For convenience, we will use T_{i }instead of if r_{i }is clear from the context. The leaf set S_{r }= {S_{1 }∪ S_{2}} where S_{r }⊂ S and S_{0 }= S  S_{r }is not represented in T_{r}. Then we can compute
for each triplet (s_{0}, s_{1}, s_{2}) ∈ (S_{0 }× S_{1 }× S_{2}).
With (s_{1}r_{1}) and (s_{2}r_{2}) we denote the known distances of s_{1 }and s_{2 }to their roots r_{1 }and r_{2}. Thus, we can compute for each triplet {s_{0}, s_{1}, s_{2}} the lengths ℓ(r_{1}r) and ℓ(r_{2}r) as
Note that, if D is additive and T_{1}, T_{2 }are isometric subtrees of T, the lengths ℓ(r_{1}r) and ℓ(r_{2}r) do not depend on the choice of the triplet {s_{0}, s_{1}, s_{2}}.
Regardless of additivity considerations, we may define the average length for a fixed s_{0 }∈ S_{0 }as
We can estimate edge lengths ℓ(r_{1}r) and ℓ(r_{2}r) by using all possible triplets as
Recovering a tree from a distance matrix
The largest path length criterion
We want to reconstruct a tree T = (V, E) with respect to a distance matrix D such that D_{T }represents D. To this end, we use triplets and the notation of a rooted tree T_{r }together with Equations 4 and 5.
Our algorithm starts with the observation that if we take an arbitrarily rooted tree T_{m }with m ∈ S and length function , then there must be a pair of leaves (neighboring leaves) that share an immediate most recent common ancestor mrca which is farthest away from the root m with respect to . In Figure 5, the pair (3, 4) satisfies this condition, we say this pair fulfills the largest path length criterion. The largest path length criterion easily generalizes to arbitrarily rooted subtrees T_{i }and T_{j }of T_{m}, where all descendants from the roots of T_{i }and T_{j }are in the vertex sets V_{i }or V_{j}, respectively.
Figure 5. The tree is rooted at leaf 5. In the tree, leaves 3 and 4 with the largest path length from their most recent common ancestor to the root 5 are neighbors.
Let be the set of rooted subtrees of T_{m }(each leaf l ∈ L_{m }is considered as a rooted subtree T_{l}). Now consider two disjoint rooted subtrees T_{i }and T_{j }of T_{m }where i, j ∈ V_{m}. Then the distance ℓ(m, mrcamS_{i}S_{j}) from the mrca of T_{i }and T_{j }to m is computed according to Equation 4, where S_{i }and S_{j }are the leaf sets of T_{i }and T_{j}, respectively. Then we pick
(, ) = argmax{ℓ(m, mrcamS_{i}S_{j})  T_{i}, T_{j }∈ } (6)
as a pair of neighbors (if we detect more than one pair, we randomly select one). By construction, (, ) fulfills the largest path length criterion.
If D is additive, ℓ(m, mrcamS_{i}S_{j}) is exactly the path length from the mrca of (T_{i}, T_{j}) to m. In other words, the path length from the mrca of (, ) to m is large stand (, ) is a true neighboring pair. However, in real applications D is rarely additive, therefore the root m is selected so as to avoid noise from stochastic errors involved with large distance estimates [17]. To this end, m is selected such that the distance from the farthest object to root m is minimal,
med = argmin_{m'∈S}{max{D(m'x)x = 1,..., n}} (7)
med is called a median object.
Moreover, to reduce the computational complexity of finding a pair of neighbors (, ) using Equation 6, we store for each T_{i }∈ its potential neighbor T_{i' }∈ such that
T_{i' }= argmax{ℓ(med, mrcamed, S_{i}, S_{j})T_{j }∈ }. (8)
Now the neighboring pair (, ) fulfilling the largest path length criterion is determined as follows
(, ) = argmax{ℓ(med, mrcamed, S_{i}, S_{i'})T_{i }∈ }. (9)
In the following, we present a natural clustering algorithm to reconstruct trees based on distance matrices and the largest path length criterion
Clustering Algorithm
• Initial step: Find the median object med using Equation 7. Set = {T_{1},..., T_{n}}  {T_{med}}. Find for each T_{i }∈ its potential neighbor T_{i' }∈ using Equation 8.
• Selection step (largest path length criterion): Find the neighboring pair (, ) using Equation 9.
• Agglomeration step: Combine and into a new rooted tree with root i_{0}j_{0}, and estimate new edge lengths of using Equation 5. Delete and and add to . Find the potential neighbor for the new rooted tree . using Equation 8, and replace T_{i' }for each T_{i }∈ by if is its potential neighbor.
• Stopping step: If  > 1 goto the Selection step, otherwise output the tree.
This algorithm is similar to approaches described elsewhere [1921], however, an essential difference is that we estimate path lengths and edge lengths by using triplets.
Local rearrangement
The heart of the clustering algorithm is the largest path length criterion, at which the path length from the mrca of (T_{i}, T_{j}) to med is estimated by ℓ(med, mrcamed, S_{i}, S_{j}) using Equation 4. Thus, as path length we take the average of the lengths obtained from at most O(n^{2 }triplets {med, s_{i}, s_{j}} ∈ med × S_{i }× S_{j}. This average may not be the representative estimate of the true path length. Moreover the root med may be too far way from the mrca and this leads to an inaccurate estimate of the path length.
To take these problems into account, we extend the clustering algorithm. To this end, imagine the algorithm has clustered T_{i }and T_{j }with corresponding disjoint leaf sets S_{i}, S_{j }⊂ S (we have finished the agglomeration step). Thus, we have created the newly rooted tree T_{{ij} }with leaf set S_{ij }= {S_{i }∪ S_{j}} and the set of unclassified objects S_{0}(S_{ij}) = S  S_{ij}. In the following, we describe the nearest neighbor interchange operation around the root of T_{i }upon condition that T_{i }consists of two rooted subtrees T_{x}, T_{y }(Figure 6a). First, we estimate average path lengths from the unclassified object set S_{0}(S_{ij}) to the mrca of (T_{x}, T_{y}), (T_{x}, T_{j}) and (T_{y}, T_{j}) as
Figure 6. Reconstruction of new rooted tree T_{{ij} }using the the preorder traversal procedure based on the largest average path length criterion. If (T_{x}, T_{y}) is the neighboring pair, we stick to the suggested grouping of T_{i }and T_{j }(see Figure 6a). Otherwise, if (T_{x}, T_{j}) or (T_{y}, T_{j}) is the neighboring pair, we switch to the trees displayed in Figure 6b or 6c, respectively.
For convenience, we will use ℓ(S_{0}(S_{ij})S_{x}S_{y}) instead of ℓ(S_{0}(S_{ij})S_{x}S_{y}S_{x}S_{y}). We now use the average path lengths from Equation 10 to decide which pair of subtrees among (T_{x}, T_{y}), (T_{x}, T_{j}) and (T_{y}, T_{j}) is preferred. More specifically, if
ℓ(S_{0}(S_{ij})S_{x}S_{y}) ≥ max{ℓ(S_{0}(S_{ij})S_{x}S_{j}), ℓ(S_{0}(S_{ij})S_{y}S_{j})}
we stick to the suggested grouping of T_{x }and T_{y }(see Figure 6a). Otherwise, if ℓ(S_{0}(S_{ij})S_{x}S_{j}) or ℓ(S_{0}(S_{ij})S_{y}S_{j}) is larger than the remaining average path lengths, we swap T_{y }and T_{j }or T_{x }and T_{j }as displayed in Figure 6b or 6c, respectively. Note that, this decision can be considered as a correction of the largest path length criterion by taking all possible triplets into account. We call the correction the largest average path length criterion.
We now explain the preorder traversal procedure [31] to reconstruct the rooted tree T_{i }using the nearest neighbor interchange operation based on the largest average path length criterion (T_{i }is a subtree of T_{{ij} }= (T_{i}, T_{j})):
Preorder traversal procedure (T_{i})
• Step 1: If T_{i }is a single leaf, return.
• Step 2: Otherwise, T_{i }consists of two subtrees T_{x }and T_{y}. Do the nearest neighbor interchange operation around the root of T_{i }based on the largest average path length criterion (Equation 10). If T_{x }and T_{j }(or T_{y }and T_{j}) were exchanged, estimate new edge lengths using Equation 5.
• Step 3: Apply the preorder traversal procedure to two rooted subtrees of T_{i}.
Representative sets and shortest triplets
For a set S of sequences (or taxa), the (genetic) distance matrix D is typically not additive due to stochastic errors [17]. Larger distances between two sequences are less accurately estimated. This leads to a low performance of both the clustering algorithm and the preorder traversal procedure for divergent data sets.
In earlier work [10,11], we have presented simple representative concepts to reduce stochastic error involved in large distances. Here, we extend our work by introducing the socalled krepresentative set concept. We use now genetic distances instead of topological distances (all edges have length 1). Our motivation is to reduce the computational complexity and to exclude objects far away from the root under consideration. In the clustering algorithm, the path length from the mrca of (T_{i}, T_{j}) to med (see Figure 7) can be estimated by two approaches. The first method picks randomly one pair (s_{i}, s_{j}) ∈ S_{i }× S_{j }then computes
Figure 7. we select only min(k, S_{i}) and min(k, S_{j}) closest leaves to the root of T_{i }and T_{j }with respect to the path length, respectively, i.e. for k = 3 we pick {1, 2} from and {4, 5, 6} from . The leaf set {1, 2} (or {4, 5, 6}) is the 3representative leaf set of the rooted subtree . (or ).
The second approach takes the average distance
Both approaches suffer from noise. Estimating the path length using Equation 11 may be inaccurate because it randomly picks a pair (s_{i}, s_{j}) which may not be really representative. Equation 12 may be problematic, especially since it might be susceptible to noise, due to the possibility of including long distances with large stochastic errors.
To overcome these problems, we select only min(k, S_{i}) and min(k, S_{j}) closest leaves to the root of T_{i }and T_{j }with respect to the path length, respectively. To illustrate, for k = 3 we pick {1, 2} from T_{i }and {4, 5, 6} from T_{j }in Figure 7.
We now define as the set of min(k, S_{i}) closest leaves to the root of T_{i}. is called the krepresentative leaf set. Hereafter, we estimate similar to Equation 4 the path length from the mrca of (T_{i}, T_{j}) to med as
which is only based on the krepresentative leaf sets. Now we can perform the clustering algorithm with reduced complexity. However, we also want to improve the preorder traversal procedure. The average path length from the unclassified object set S_{0}(S_{ij}) to the mrca of (T_{i}, T_{j}) is estimated by Equation 10 which also suffers from noise. To overcome this problem, we select only min(k, S_{0}(S_{ij})) unclassified objects closest to the root of tree T_{{ij} }with respect to distances where s_{0 }∈ S_{0}(S_{ij}). We call the subset, denoted (S_{ij}), krepresentative unclassified object set.
We now define a shortest triplet, which contains three representatives of the three krepresentative sets. By construction, are close to the root of T_{{ij} }and close to each other. Therefore, the edge length estimates based on shortest triplet {} are less susceptible to estimation errors.
We now rewrite Equation 10 to estimate the average path length from the representative unclassified object set (S_{ij}) to the mrca of (T_{i}, T_{j}) using only shortest triplets as
In short, the preorder traversal procedure uses only shortest triplets to estimate path lengths as well as edge lengths.
Shortest triplet clustering algorithm (STC)
We introduce now the shortest triplet clustering algorithm by combining the clustering algorithm, the local rearrangement, the krepresentative sets, and the shortest triplets approach.
Shortest triplet clustering algorithm (STC)
• Initial step:
 (i): Find the median object med using Equation 7.
 (ii): Set = {T_{1},..., T_{n}}  {T_{med}} and for each T_{i }∈ its representative leaf set = {i}.
 (iii): Find for each T_{i }∈ its potential neighbor T_{i' }∈ using Equation 8.
• Selection step (largest path length criterion): Find the neighboring pair (, ) using Equation 9.
• Agglomeration step:
 (i): Combine and into a new rooted tree with root i_{0}j_{0}, and estimate new edge lengths of using Equation 5 based on shortest triplets.
 (ii): Compute the krepresentative leaf set of . based on krepresentative leaf sets and of and , respectively.
 (iii): Compute the krepresentative unclassified object set of .
 (iv): Delete and and add to .
 (v): Find the potential neighbor for the new rooted tree using Equation 8 based on representative sets, and replace T_{i' }for each T_{i }∈ by if is its potential neighbor.
• Local rearrangement step: Apply the preorder traversal procedure to the rooted subtrees and of the new rooted tree based on only shortest triplets.
• Stopping step: If  > 1, goto Selection step, otherwise output the tree.
The complexity of STC
Now we briefly describe the complexity of the STC. At the initial step, (i), (ii), and (iii) are done in O(n^{2}), O(n) and O(n^{2}) time, respectively. Thus, the complexity of the initial step is O(n^{2}). The selection step is done in O(n). At the agglomeration step, (i), (ii), (iii), (iv), and (v) are done in O(k^{3}), O(k), O(nk^{2}), O(1), and O(nk^{2}) time, respectively. Thus, the complexity of the agglomeration step is O(nk^{2 }+ k^{3}). Finally, we are estimating the complexity of the preorder traversal procedure based on only shortest triplets. Step 1 is done in constant time. Step 2, the nearest neighbor interchange operation around the root of T_{i }costs O(k^{3}). Estimating new edge lengths is done in O(k^{3}) time. Recomputing the krepresentative leaf set of T_{i }based on krepresentative leaf sets of its rooted subtrees T_{x }and T_{y }costs O(k) time. Finally, recomputing the krepresentative unclassified object set (S_{i}) of T_{i }based on the krepresentative leaf set of T_{j }and the krepresentative unclassified object set (S_{ij}) of T_{{ij} }is done in O(k) time. Thus, the complexity of step 2 is O(k^{3}). Step 3 is done in constant time. Step 1, step 2, and step 3 are repeated O(n) times so the complexity of the preorder traversal procedure is O(nk^{3}).
In the STC algorithm, the selection step, the agglomeration step and the local rearrangement step are repeated (n  2) times so the overall complexity of the STC algorithm is O(n^{2}k^{3}). Practically, we chose k = 5 as a good compromise between the accuracy and computational complexity for all data sets. That is, the practical complexity of the STC algorithm is only O(n^{2}).
Authors' contributions
Both authors participated in the design of the study and performed the analysis. LSV implemented the algorithms. Both authors wrote and approved the final manuscript.
Acknowledgements
We would like to express special thanks to Heiko Schmidt for his technical support. We thank Gunter Weiss, Ingo Ebersberger, Tanja Gesell and Jutta Buschbom for carefully reading the manuscript. We acknowledge the use of supercomputing resources of the ZAM/NIC at the Research Center Jiilich. We thank three anonymous referees for helpful comments.
References

Saitou N, Nei M: The Neighbor – joining Method: A New Method for Reconstructing Phylogenetic Trees.
Mol Biol Evol 1987, 4:406425. PubMed Abstract  Publisher Full Text

Strimmer K, von Haeseler A: Quartet Puzzling: A Quartet Maximum – Likelihood Method for Reconstructing Tree Topologies.

Gascuel O: BIONJ: An Improved Version of the NJ Algorithm Based on a Simple Model of Sequence Data.
Mol Biol Evol 1997, 14:685695. PubMed Abstract  Publisher Full Text

Huson DH, Nettles SM, Warnow TJ: DiskCovering, a FastConverging Method for Phylogenetic Reconstruction.
J Comput Biol 1999, 6:369386. PubMed Abstract  Publisher Full Text

Bruno WJ, Socci ND, Halpern AL: Weighted Neighbor Joining: A Likelihood BasedApproach to DistanceBased Phylogeny Reconstruction.
Mol Biol Evol 2000, 17:189197. PubMed Abstract  Publisher Full Text

Desper R, Gascuel O: Fast and Accurate Phylogeny Reconstruction Algorithms Based on the MinimumEvolution Principle.
J Comput Biol 2002, 9:687706. PubMed Abstract  Publisher Full Text

Csürös M: Fast Recovery of Evolutionary Trees with Thousands of Nodes.
J Comput Biol 2002, 9:277297. PubMed Abstract  Publisher Full Text

Guindon S, Gascuel O: A Simple, Fast and Accurate Algorithm to Estimate Large Phylogenies by Maximum Likelihood.
Syst Biol 2003, 52:696704. PubMed Abstract  Publisher Full Text

Stamatakis A, Ludwig T, Meier H: RAxMLIII: A fast program for maximum likelihoodbased inference of large phylogenetic trees.
Bioinformatics 2005, 21:456463. PubMed Abstract  Publisher Full Text

Vinh LS, von Haeseler A: IQPNNI: Moving fast through tree space and stopping in time.
Mol Biol Evol 2004, 21:15651571. PubMed Abstract  Publisher Full Text

Vinh LS, Schmidt HA, von Haeseler A: PhyNav: A Novel Approach to Reconstruct Large Phylogenies. In Proceedings of the 28th Annual German Classification Society Conference (GfKl 2004). Dortmund, Germany; 2004:in press.

Charleston MA: HitchHiking: A Parallel Heuristic Search Strategy, Applied to the Phylogeny Problem.
J Comput Biol 2001, 8:7991. PubMed Abstract  Publisher Full Text

Brauer MJ, Holder MT, Dries LA, Zwickl DJ, Lewis PO, Hillis DM: Genetic Algorithms and Parallel Processing in MaximumLikelihood Phylogeny Inference.
Mol Biol Evol 2002, 19:17171726. PubMed Abstract  Publisher Full Text

Schmidt HA, Strimmer K, Vingron M, von Haeseler A: TREEPUZZLE: Maximum likelihood phylogenetic analysis using quartets and parallel computing.
Bioinformatics 2002, 18:502504. PubMed Abstract  Publisher Full Text

Schmidt HA, von Haeseler A: Maximum Likelihood Analysis Using TREEPUZZLE. In Current Protocols in Bioinformatics. Edited by Baxevanis AD, Davison DB, Page RDM, Stormo G, Stein L. New York, USA: Wiley and Sons; 2003:6.6.16.6.25.

CavalliSforza L, Edwards AWF: Phylogenetic analysis: Models and estimation procedures.
Am J Hum Genet 1967, 19:233257. PubMed Abstract

Fitch W, Margoliash E: Construction of Phylogenetic trees.
Science 1967, 155:279284. PubMed Abstract

Hartigan AJ: Clustering Algorithms. John Wiley and Sons, Inc; 1975.

Farris J: On the phenetic approach to vertebrate classification.
1977, 17:823850.

Klotz NKRB LC, Mitchell RM: Calculation of evolutionary trees from sequence data.
Proc Natl Acad Sci USA 1979, 76:45164520. PubMed Abstract  PubMed Central Full Text

Li WH: Simple method for constructing phylogenetic trees from distance matrices.
Proc Natl Acad Sci USA 1981, 78:10851089. PubMed Abstract  PubMed Central Full Text

Rambaut A, Crassly NC: SeqGen: An application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees.
Comput Appl Biosci 1997, 13:235238. PubMed Abstract

Kimura M: A Simple Method for Estimating Evolutionary Rates of Base Substitutions through Comparative Studies of Nucleotide Sequences.
J Mol Evol 1980, 16:111120. PubMed Abstract

Harding EF: The probabilities of rooted treeshapes generated by random bifurcation.

Swofford DL, Olsen GJ, Waddell PJ, Hillis DM: Phylogeny Reconstruction. In Molecular Systematics. 2nd edition. Edited by Hillis DM, Moritz C, Mable BK. Sunderland, Massachusetts: Sinauer Associates; 1996:407514.

Felsenstein J: PHYLIP (Phylogeny Inference Package) version 3.5c. [http://evolution.genetics.washington.edu/phylip.html] webcite
Department of Genetics, University of Washington, Seattle 1993.

Robinson DR, Foulds LR: Comparison of phylogenetic trees.
Mathematical Biosciences 1981, 53:131147. Publisher Full Text

Buneman P: The recovery of trees from measures of dissimilarity. In Mathematics in the archaeological and historical sciences. Edited by Hodson, Lendall, Tautu. Edinburgh: Edinburgh university press; 1971.

Semple C, Steel M: Phylogenetics. OXFORD univerity press; 2003.

Felsenstein J: Inferring Phylogenies. Sunderland, Massachusetts: Sinauer Associates; 2004.

Aho AV, Hopcroft JE, Ullman JD: The Design and Analysis of Computer Algorithms. AddisonWesley Publishing Company; 1974.