FastTree 2 – Approximately Maximum-Likelihood Trees for Large Alignments — TEI-tagged view

Morgan N. PricepersName, Paramvir S. DehalpersName, Adam P. ArkinpersName (2010date). “FastTree 2 – Approximately Maximum-Likelihood Trees for Large Alignments.” PLoS ONEtitle·j.

Object of study — scalable approximately maximum-likelihood phylogenetic inference for very large DNAterm and protein sequenceterm alignments

Research question — whether restricted tree-search heuristicsterm can retain phylogenetic accuracy while making maximum-likelihood-scale inference practical for hundreds of thousands of sequences

Contribution — FastTree 2term, its minimum-evolutionterm SPR and maximum-likelihood NNI pipeline, CATterm rate approximation, local support procedure, and comparative accuracy/scalability evaluation

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

FastTree 2term – Approximately Maximum-Likelihood Trees for Large Alignments Abstract We recently described FastTree, a tool for inferring phylogeniesterm for alignments with up to hundreds of thousands of sequences. Here, we describe improvements to FastTree that improve its accuracy without sacrificing scalability.

Where FastTree 1 used nearest-neighbor interchangesterm (NNIs) and the minimum-evolutionterm criterion to improve the tree, FastTree 2term adds minimum-evolutionterm subtree-pruning-regraftingterm (SPRs) and maximum-likelihood NNIs. FastTree 2term uses heuristicsterm to restrict the search for better trees and estimates a rate of evolutionterm for each site (the “CATterm” approximation). Nevertheless, for both simulated and genuine alignments, FastTree 2term is slightly more accurate than a standard implementation of maximum-likelihood NNIs (PhyMLterm 3 with default settings). Although FastTree 2term is not quite as accurate as methods that use maximum-likelihood SPRs, most of the splits that disagree are poorly supported, and for large alignments, FastTree 2term is 100–1,000 times faster. FastTree 2term inferred a topologyterm and likelihood-based local support valuesterm for 237,882 distinct 16S ribosomal RNAsterm on a desktop computer in 22 hours and 5.8 gigabytes of memoryterm·var.

FastTree 2term allows the inference of maximum-likelihood phylogeniesterm for huge alignments. FastTree 2term is freely available at http://www.microbesonline.org/fasttree .

Introduction Inferring evolutionary relationships, or phylogeniesterm, from families of related DNAterm or protein sequencesterm is a central method in computational biology. Sequence-based phylogeniesterm are widely used to understand the evolutionary relationships of organisms and to analyze the functions of genes.

The largest gene familiesterm already contain tens to hundreds of thousands of representatives, and with the rapid improvements in DNAterm sequencing, we expect even larger data sets to arrive soon. Large families can be aligned with profile-based methods that scale linearly with the number of sequencesterm·var ( http://hmmer.janelia.org/ ; [1] ). However, most methods for inferring phylogeniesterm from these alignments scale as O( ) or worse, where is the number of sequencesterm·var and is the length (width) of the alignment. Thus, inferring phylogeniesterm has become computationally challenging.

We recently described a scalable method for inferring phylogeniesterm, FastTree 1.0term [2] . FastTree 1.0term is based on the “minimum-evolutionterm” principle – it tries to find a topologyterm that minimizes the amount of evolution, or the sum of the branch lengthsterm. FastTree 1.0term uses a heuristicterm variant of neighbor joiningterm [3] , [4] to quickly find a starting treeterm·var and uses nearest-neighbor interchangesterm (NNIs) to refine the topologyterm. (A nearest-neighbor interchangeterm·def swaps a node and its neighbor; for example, it might change ((A,B),C,D) to ((A,C),B,D) or ((A,D),B,C).) FastTree implements these operations in O( ) space, where is the number of characters in the alphabet, by storing profiles of subtreesterm instead of distances between them. This requires far less memoryterm·var than storing the pairwise distancesterm, which is necessary for neighbor joiningterm and related approaches. This also allows for heuristicsterm that reduce the theoretical running timeterm·var to O( ). (FastTree 1.0term also included some O( ) steps, but these have since been removed, see http://www.microbesonline.org/fasttree/ChangeLog .) In comparison, computing all pairwise distancesterm, which is required with most minimum-evolutionterm approaches, requires O( ) time. The main limitation of FastTree 1.0term, as compared to other minimum-evolutionterm methods, is that it does not correct distances for multiple substitutions during its initial neighbor joiningterm phase. However, this is more than made up for by the NNIs. In practice, FastTree 1.0term is more accurate than most other minimum-evolutionterm methods, but not as accurate as maximum-likelihood methods [2] .

In the maximum-likelihood (ML) approachterm, evolution is explicitly modeled with a transition rate matrixterm, and the tree that best explains the data – the tree with the highest likelihoodterm – is the best tree [5] . The ML criterion ranks the trees but does not specify how to find a good topologyterm. Because ML phylogenetic inference is NP complete [6] , no practical method can guarantee that it will find the optimal topologyterm for a large alignment. The most scalable ML methods, such as PhyMLterm and RAxMLterm, begin with a starting treeterm·var produced by a faster method, and try to increase the likelihoodterm by optimizing individual branch lengthsterm and performing local rearrangements [7] – [9] . By re-optimizing only a few branch lengthsterm at each move, the cost of considering or performing a move can be reduced to O( ) time, where is the size of the alphabet. However, in practice, the number of moves grows as roughly O( ), and the optimization steps are inherently slow because they require numerical solving and iteration. This explains why both PhyMLterm and RAxMLterm can take over a day for just 1,000 protein sequencesterm. Estimating the reliability of the tree with the bootstrapterm [10] generally increases the computational requirements another 100-fold (although this can be reduced by reusing computations across replicates [11] ).

Here, we describe FastTree 2term, a tool for inferring ML trees for large alignments. Besides constructing an initial tree with neighbor joiningterm and improving it with minimum-evolutionterm NNIs, FastTree 2term uses minimum-evolutionterm subtree-pruning-regraftingterm (SPRs) [8] , [12] and ML NNIs to further improve the tree. (In subtree-pruning-regraftingterm, a subtree is removed from the tree and reinserted elsewhere, e.g., pruning and regrafting C might change ((A,B),(C,D),E) to ((A,(B,C)),D,E).) FastTree 2term uses heuristicsterm to reduce the search space and hence to maintain the scalability of both stages. Another justification for reducing the search space is that intensive tree search often finds small improvements in the tree's length or likelihoodterm, but these changes may not be statistically or biologically significant (e.g., [13] ). Briefly, FastTree's key heuristicsterm are: It uses “linear SPRs” to consider just of the possible SPR moves. At each node, it examines the shortest SPRs first and then extends the most promising candidates.

It searches for SPR moves for every subtree just twice, instead of iterating until convergenceterm·def.

During the ML phase, it limits the ML NNIs to at most rounds; in practice, it converges before it reaches this limit, but the limit ensures a predictable running timeterm·var.

It limits the effort to optimize model parameters and branch lengthsterm.

It abandons optimization for NNI moves that seem, after partial optimization, to significantly lower the likelihoodterm.

It does not try to improve parts of the tree that did not improve in recent rounds.

To account for the variation in rates across sites, FastTree uses the “CATterm” approximation [14] rather than the standard discrete gamma model with four rates ( ) [15] . Some sites evolve much more slowly than others, and the ideal way to account for this is to integrate the likelihoodterm at each site over the (unknown) relative evolutionary rate of that site, using a prior distribution over the relative rates such as a gamma distributionterm. However, these integrals are analytically intractable and computationally prohibitive. The “ ” approach is to use four rate categoriesterm to approximate the continuous gamma distributionterm. However, still requires four times more CPU time and memoryterm·var than a model with no rate variation across sites. Furthermore, for large alignments, the data tightly constrains the rate at each site. Thus, it is much faster, and just as accurate, to use a good estimate of the rate at each site (CATterm) rather than to sum over four potential rates ( ) [14] . FastTree selects the most likely rate for each site from among 20 fixed possibilities.

Because of the heuristicsterm, FastTree 2term is not guaranteed to reach a locally optimal likelihoodterm in tree space. However, at each step it does guarantee that the likelihoodterm increases (under the CATterm approximation). Thus, FastTree 2term is an approximately-maximum-likelihood methodterm.

We will show that in practice, FastTree 2term is slightly more accurate than a standard implementation of maximum-likelihood NNIs, PhyMLterm 3 with default settings [16] , [17] . Specifically, in simulations, FastTree 2term recovers a higher proportion of true splits, and on genuine alignments, FastTree 2term's topologies tend to have higher likelihoodsterm. FastTree's minimum-evolutionterm SPR moves give it a better starting treeterm·var than PhyMLterm's starting treeterm·var, which is obtained with BIONJ (a weighted variant of neighbor joiningterm [18] ). This more than makes up for FastTree's heuristicsterm, which reduce the intensity of search for ML NNIs but have little effect on accuracy. We also confirm that using the CATterm approximation instead of the model (which is itself an approximation of the continuous gamma distributionterm) has little effect on the quality of the tree.

Although FastTree 2term is significantly less accurate than ML methods that use SPR moves, such as PhyMLterm with slower settings or RAxMLterm, most of the splits that disagree are poorly supported, and FastTree is much faster. FastTree 2term can analyze alignments with tens or hundreds of thousands of sequences in under a day on a desktop computer. For alignments with 500 sequences or more, FastTree 2term is at least 100 times faster than either PhyMLterm 3.0 or RAxMLterm 7.2.1. FastTree 2term is faster than RAxMLterm 7 mostly because of less intensive ML search (NNIs instead of SPRs) and because RAxMLterm 7 optimizes branch lengthsterm under the model. However, FastTree also has a faster starting treeterm·var, and it initially increases the likelihoodterm more quickly than RAxMLterm 7 does.

Because of its speed, FastTree 2term is suitable for bootstrapping. However, to provide a quicker estimate of the tree's reliability, FastTree 2term provides local support valuesterm based on the Shimodaira-Hasegawa (SH) test [16] , [17] , [19] . FastTree 2term should be useful for reconstructing the tree of life and for analyzing the millions of uncharacterized proteins that are being identified by genome sequencing.

Results We compared FastTree's speed and accuracy to those of PhyMLterm 3.0 and RAxMLterm 7, the most popular maximum-likelihood methods. To measure the quality of the resulting trees, we measured the topological accuracyterm on simulated alignments and the likelihoodterm on genuine biological alignments.

Topological Accuracyterm in Simulations We tested FastTree on simulated protein alignments with 250 to 5,000 sequences [2] . These simulations were derived from diverse gene familiesterm that arise in genome-scale studies (“Collections of Orthologous Groups” or COGs, [20] ). The simulations include varying evolutionary rates across sites and include realistic placement of gaps. The simulations are available from the FastTree web site ( http://microbesonline.org/fasttree/#Sims ).

We defined the topological accuracyterm as the proportion of the splits in the true trees that are recoveredterm by each method. This is the converse of the topological (“Robinson-Foulds”) distance, scaled to range from 0 to 1. As shown in Table 1 , FastTree 2term was slightly more accurate than PhyMLterm 3 with default settings (NNI search), and much more accurate than minimum-evolutionterm or parsimony methods, but not as accurate as ML methods that use SPR moves. The differences in accuracy between FastTree 2term and the other methods were statistically significant (all , paired tests).

Topological accuracyterm of trees inferred from simulated alignments.

To test the practical significance of the additional true splits that are found by using ML SPR moves, we examined the local support valuesterm reported by PhyMLterm 3. We defined “strongly supportedterm” as having both SH-like local supports and approximate likelihood ratio test (aLRT) supportsterm [21] of 95% or higher. Only 16% of the true splits that are found by PhyMLterm 3 with SPR moves but missed by FastTree 2term were strongly supportedterm. The full distribution of support values is shown in Figure 1 . Conversely, among the strongly supportedterm splits that were found by PhyMLterm 3 with SPRs but not FastTree, 20% were incorrect. Thus, few of the additional true splits have high support, and of the splits that disagree, even the ones that have high support have a significant probability of being incorrect.

Local support valuesterm for splits found by PhyMLterm with SPR moves and/or FastTree. We examined local support valuesterm for the splits inferred by PhyMLterm 3.0 with + SPRs on simulated alignments with 250 protein sequencesterm. We classified PhyMLterm's splits as correct and found by both PhyMLterm and FastTree, correct but missed by FastTree, or incorrect. We show the distribution of support values for each class. The right-most bin includes the strongly supportedterm splits (0.95 to 1.0), and the gray dashed line shows the uniform distribution. The support values are PhyMLterm's minimum of the approximate likelihoodterm ratio test [21] and SH-like [16] , [17] local supports.

To understand why FastTree 2term was outperforming PhyMLterm 3 with NNI search, we ran PhyMLterm 3 with FastTree's minimum-evolutionterm tree as its starting treeterm·var. For the protein simulations with 250 sequences, this improved PhyMLterm's accuracy to 86.8%, which is statistically indistinguishable from FastTree's accuracy of 86.9% ( , paired test, ). We also confirmed that FastTree's minimum-evolutionterm phase yields more accurate starting treesterm·var than either PhyMLterm 3's approach of using BIONJ with maximum-likelihood distances or RAxMLterm's implementation of parsimony ( Table 1 ).

CAT-Based Branch Lengthsterm and Local Support Valuesterm Because FastTree 2term does not exhaustively optimize the likelihoodterm, and because it reports branch lengthsterm and local support valuesterm that were estimated using the CATterm approximation, we compared its branch lengthsterm and local support valuesterm to -based lengths and supports. Specifically, for the protein simulations with 250 sequences, we re-optimized the branch lengthsterm and computed local SH-like support values for the FastTree topologies with PhyMLterm 3 and . (For both tools, we used the Jones-Taylor-Thorton (JTT) model of amino acid evolution.) PhyMLterm's internal branch lengthsterm were well correlated with those from FastTree ( = 0.90). For branch lengthsterm of 1.0 or less, the average difference was just 0.01, and for branch lengthsterm between 0.01 and 1.0, the average percent difference was 13%. For internal branch lengthsterm on correct splits, FastTree agreed slightly better with the true lengths (median absolute difference of 0.062 for PhyMLterm and 0.059 for FastTree). Thus, the CATterm approximation gave acceptable branch lengthsterm.

If accurate branch lengthsterm are essential, however, then neither the CATterm approximation nor the standard approximation is sufficient. The approximation was introduced for alignments with just 10 sequences, and four discrete rate categoriesterm may not suffice to give accurate likelihoodsterm on larger alignments [15] , [22] . For alignments of 16S ribosomal RNAsterm, branch lengthsterm can be a factor of two shorter than branch lengthsterm ( Figure S1 ). As explained in Figure S1 , correcting by the average posterior rate reduces this problem, and FastTree can compute a fast but accurate approximation to -based lengths.

The local SH-like support values also showed a good correlation between FastTree and PhyMLterm ( = 0.90). For splits with local support valuesterm of at least 0.9 from either FastTree or PhyMLterm, the average absolute difference was just 0.008, which is not much greater than the sampling error. (For example, with 1,000 bootstrapsterm and 95% support, the standard deviation of support values due to sampling is = 0.007.) FastTree was less effective than PhyMLterm in distinguishing correct from incorrect splits, but the difference was slight: the area under the receiver operating curve (AOC) was 0.880 instead of 0.887 ( , test of [23] ).

Effectiveness of Heuristicsterm We then examined how the topological accuracyterm of FastTree 2term is affected by its heuristicsterm. As shown in Table 1 , the minimum-evolutionterm phase of FastTree, which uses linear SPRs, is not as accurate as FastMEterm 2, a minimum-evolutionterm method that performs exhaustive SPR moves [8] , [12] . FastMEterm computes distances between internal nodes differently from the minimum-evolutionterm phase of FastTree: FastMEterm uses averages of distances between sequences, while FastTree uses distances between profiles, which are averages of sequences. Nevertheless, FastTree 1 with only NNI moves gave very similar results as FastMEterm with only NNI moves [2] . Thus, we attribute the modest difference in accuracy of the minimum-evolutionterm methods with SPRs to FastTree's heuristicsterm. To eliminate this effect, we ran FastTree with the FastMEterm starting treeterm·var. To eliminate the effect of FastTree's ML heuristicsterm, we ran it with exhaustive ML NNIs, and with more exhaustive optimization of branch lengthsterm within each NNI (4 rounds of optimizing branch lengthsterm for each quartet, instead of 1–2 rounds). In combination, FastTree 2term with FastMEterm+SPR starting treesterm·var and exhaustive NNIs improved the accuracy on simulated alignments with 5,000 protein sequencesterm from 84.3% to 85.0%. This modest effect illustrates that all of FastTree's heuristicsterm have little effect on accuracy, and that removing them would improve the topologyterm little relative to adding ML SPRs (e.g., RAxMLterm 7.2.1 was 88.4% accurate).

We also tested FastTree on simulations with over 78,000 nucleotide sequences. These simulations are derived from a 16S ribosomal RNAterm alignment (see Methods ). The large size of these simulated alignments makes them a stringent test of FastTree's heuristicsterm. In these simulations, FastTree gave much more accurate topologies than exact neighbor joiningterm or Clearcut [24] , a faster heuristicterm variant of neighbor joiningterm ( Table 1 ). (To analyze such large alignments with exact neighbor joiningterm, we used NINJA [25] .) To verify that the heuristicsterm in FastTree's neighbor joiningterm phase do not reduce accuracy, we also ran FastTree with the exact neighbor-joining tree as its starting treeterm·var, before doing minimum-evolutionterm NNIs and SPRs and ML NNIs. This gave the same accuracy as the regular FastTree or as FastTree with the fastest settings of its heuristicsterm for the neighbor joiningterm phase (-fastest). All three variants found 92.10% of splits correctly.

It may seem surprising that FastTree can reach accurate topologies when it does not compare all pairs of sequences to each other. However, minimum-evolutionterm NNIs and SPRs are “consistent” – they find correct trees, even if the distances contain some errors, as long as the errors are much smaller than the internal branch lengthsterm [26] , [27] . In practice, the errors are often larger than the internal branch lengthsterm, but this still probably explains why NNIs and SPRs suffice to find most of the splits correctly.

Quality of Topologies for Biological Alignments To confirm that FastTree finds good topologies for genuine alignments, and not just in simulations, we tested it on 16S ribosomal RNAsterm and on protein families from COGterm. Although these families are quite large (up to 300,000 or 19,000 members, respectively), we first tested random subsets of just 500 sequences, so that we could run PhyMLterm 3 with . To measure the quality of the topologies from FastTree 2term, PhyMLterm 3, and RAxMLterm 7, we re-optimized the branch lengthsterm with a model (using RAxMLterm) and compared the resulting likelihoodsterm. As expected from the simulations, FastTree found better topologies than PhyMLterm 3 with and NNI moves, but not as good as RAxMLterm 7 ( Table 2 ).

Average log-likelihood for genuine alignments with 500 sequences.

We then tested FastTree and RAxMLterm on larger alignments of 16S rRNAs and COGs. For alignments with thousands of sequences, RAxMLterm 7.0.4 is a bit slow, so we used RAxMLterm 7.2.1, which introduced a fast convergenceterm·def option as well as other optimizations. With fast convergenceterm·def, RAxMLterm terminates the search if less than 1% of splits change during a round of SPR moves. As shown in Table 1 , for 5,000 proteins, RAxMLterm with fast convergenceterm·def is nevertheless quite accurate.

On the larger alignments, RAxMLterm 7.2.1's likelihoodsterm were much higher than FastTree's, and all of the differences in likelihoodterm were statistically significant (all , SH testterm using CONSEL [28] ). However, FastTree did find most of the splits in the RAxMLterm topologyterm that had strong support ( Table 3 ). For example, FastTree found 96–98% of RAxMLterm's splits that had global bootstrapterm of 90% or higher.

Comparison of RAxMLterm and FastTree's log likelihoodsterm, and the agreement of FastTree with RAxMLterm's well-supported splits, for large genuine alignments.

Running Timeterm·var and Memoryterm·var Required Finally, we compared the computational performance of FastTree, RAxMLterm, and PhyMLterm, on genuine alignments. As shown in Table 4 , for alignments with 500 sequences, FastTree is about 100 times faster than RAxMLterm 7.0.4 when using the same model of evolution, and even faster relative to PhyMLterm 3. For alignments with thousands of sequences, FastTree was still 100–800 times faster than RAxMLterm 7.2.1 with fast convergenceterm·def of SPRs, while PhyMLterm 3 did not complete in a reasonable amount of time.

Running timeterm·var and memoryterm·var usage on genuine alignments.

For one of the largest alignments existing today, containing 237,882 16S ribosomal RNAsterm, FastTree took less than a day and 5.8 GB of memoryterm·var on a desktop computer. For comparison, given that RAxMLterm took over 2 days for just 15,011 sequences, and optimistically assuming O( ) scaling, RAxMLterm would take around half a year for the full 16S alignment. Analyzing such alignments with traditional minimum evolutionterm·def approaches based on a distance matrix would also be prohibitive – just computing and storing all pairwise distancesterm for these sequences, without computing a topologyterm, would require roughly a day and a half and 113 GB of storage.

All of the FastTree times include the computation of local SH-like support values, while the other tools were run without support values. The local support valuesterm do not affect FastTree's running timeterm·var much. For example, across seven COG alignments with 2,500 protein sequencesterm each, the average time for FastTree to infer a tree is 345 seconds, and the average time for it to compute SH-like supports is 51 seconds. For the full alignment of 237,882 16S rRNAs, the supports required just one hour.

Much of the time in RAxMLterm 7.2.1 is spent optimizing the branch lengthsterm under the model, even though the CATterm approximation is used to search for a good topologyterm. (RAxMLterm can also perform SPR moves under the model, but we ran it with SPR moves under the CATterm model only, followed by optimizing branch lengthsterm under , because RAxMLterm 7.2.1 does not report CAT-based branch lengthsterm.) If branch lengthsterm are not required, such as during bootstrapping, then RAxMLterm can be 2–3 times faster than shown in Table 4 . For example, for 15,011 16S rRNAs, if the phase is removed, then RAxMLterm 7.2.1 takes 30 hours instead of 64 hours, which is still about 45 times slower than FastTree. The phase of RAxMLterm is also expected to quadruple the memoryterm·var required. For example, for 15,011 16S rRNAs, FastTree required 0.56 GB of memoryterm·var, while RAxMLterm with required 2.6 GB.

Improvement of Likelihoodterm Over Time To compare the search strategies of FastTree and RAxMLterm more directly, we compared their improvement in likelihoodsterm over time for a nucleotide alignment of 4,114 16S rRNAs [11] and for seven protein alignments of COG families with 2,500 members. We ran both methods with the CATterm approximation and with either the generalized time-reversible (GTR) model of nucleotide substitution or the JTT model of amino acid substitution. We computed likelihoodsterm for intermediate and final trees with RAxMLterm, re-optimized branch lengthsterm, and . Figure 2 shows the running timeterm·var and log likelihoodterm for FastTree's minimum-evolutionterm and final tree, for RAxMLterm's initial parsimony tree and successive rounds of SPR moves, and also for RAxMLterm with FastTree's minimum-evolutionterm tree as its starting treeterm·var. These times do not include FastTree's support values or RAxMLterm optimizing branch lengthsterm under .

Likelihoodsterm over time for genuine alignments. Each line shows the time it takes a different tool to reach a given likelihoodterm. For the COG alignments, all times and likelihoodsterm are averages over the seven alignments. For FastTree, we show the time and the improvement in likelihoodterm for the minimum-evolutionterm topologyterm and the final (approximately-ML) topologyterm. For RAxMLterm, we show the maximum parsimony starting topologyterm, the first two rounds of SPR moves, and the final topologyterm (note the interrupted axis). For RAxMLterm with FastTree's (minimum-evolutionterm) starting treeterm·var, we show the starting topologyterm and RAxMLterm's first two rounds of SPR moves.

Given the same starting treeterm·var, FastTree's ML phase improved the likelihoodterm by roughly the same amount as one round of RAxMLterm's SPR moves, and in about 40% of the time ( Figure 2 ). FastTree's ML phase also performs about as well as one round of RAxMLterm's SPR moves in finding well-supported splits ( Figure S2 ). We obtained similar results for other large 16S alignments ( Table S1 ). Although this comparison shows that FastTree is initially faster than RAxMLterm, the RAxMLterm's first round of SPR moves is only a fraction of its run time. Most of the difference in speed between FastTree and RAxMLterm is because of RAxMLterm's more thorough search for a better topologyterm and because of RAxMLterm's branch lengthsterm.

Starting Treesterm·var: Minimum-Evolutionterm versus Maximum Parsimony RAxMLterm's parsimony phase was 4–17 times slower than FastTree's minimum evolutionterm·def phase, and generally slower than FastTree with ML NNIs. FastTree's speed advantage grows with larger alignments (data not shown), which is expected because FastTree should scale as O( ) and RAxMLterm's parsimony phase uses randomized stepwise addition, which scales as O( ), as well as limited parsimony-based SPR moves. There are faster implementations of parsimony, such as RAxMLterm 7.2.5 (which was released after we conducted the above experiments) or TNT [29] , but these still scale as O( ). For 15,011 16s RNAs, RAxMLterm 7.2.5's parsimony and FastTree's minimum evolutionterm·def phase take about the same time (data not shown).

As measured by likelihoodterm, FastTree's minimum-evolutionterm starting treesterm·var were much better than RAxMLterm's parsimony starting treesterm·var for the COG alignments, but much worse for large 16S rRNA alignments ( Figure 2 and Table S1 ). The differences in likelihoodterm reflects the criterion, and not merely differences in the search strategy: for the COG alignments, the RAxMLterm parsimony starting treesterm·var were more parsimonious than FastTree's minimum-evolutionterm trees (average parsimony scores of 281,237 and 283,125, respectively). Conversely, for the 16S alignment with 4,114 sequences, FastTree's minimum-evolutionterm tree was shorter than the parsimony tree (lengths of 43.0 and 44.6, respectively). For this alignment, the minimum-evolutionterm tree's log likelihoodterm was 2,705 worse than parsimony's, yet minimum evolutionterm·def found more of the strongly-supported splits in the final RAxMLterm tree: minimum evolutionterm·def found 826 of the 851 splits with a global bootstrapterm 90%, while parsimony found 814 of them. Thus, we are not sure if the difference in likelihoodterm is biologically meaningful.

Discussion We have shown that FastTree 2term computes accurate topologies in a reasonable amount of time for alignments with up to hundreds of thousands of sequences. FastTree is open source software and is available at http://microbesonline.org/fasttree . The C source code is extensively documented and contributions are welcome. FastTree trees for every microbial gene family, including families with tens of thousands of members such as ABC transporters, are available at MicrobesOnline ( http://microbesonline.org/ ), along with a “tree-browser” for examining these trees. These trees will be updated from FastTree 1 to FastTree 2term in the next release of MicrobesOnline.

Because DNAterm sequencing technology is improving rapidly, we expect to have alignments with millions of sequences soon. For these huge alignments, the most computationally demanding step will be the initial neighbor-joining phase. In FastTree 2term.0, which is described here, neighbor joiningterm takes O( ) time and O( ) space, while the other stages take at most O( ) time and O( ) space. For example, for 237,882 16S sequences, the neighbor-joining phase of FastTree 2term.0 already takes 10.8 of the 21.8 hours. In FastTree 2term.1, we have improved the scaling of time and memoryterm·var from O( ) to O( ), without affecting accuracy in our simulations (data not shown). FastTree 2term.1 also supports parallel execution of the key steps in the neighbor-joining phase. To improve scalability further, it might be possible to use a divide-and-conquer method to find clusters of closely related sequences in O( ) time, as in PartTree [30] . In our simulations, PartTree starting treesterm·var do not allow FastTree to reach the same accuracy as FastTree's neighbor-joining starting treeterm·var does (data not shown), but a divide-and-conquer approach might still suffice to obtain a partially resolved initial tree.

Such huge families also raise challenges for multiple sequence alignment. We have used profile alignmentterm to avoid the challenges of multiple sequence alignment on large families. This works well for 16S RNAs because Infernal takes advantage of highly conserved secondary structure [1] , but we are not sure that it gives accurate results for diverse protein families. In contrast, traditional progressive multiple sequence alignment methods are not scalable because their output grows as O( ): there are O( ) independent insertions, and each insertion requires a new column in the alignment and hence O( ) storage. However, Fast Statistical Alignment uses an O( ) representation, both internally and as an output format [31] . Combining this representation with fast guide tree construction, it should be possible to build progressive multiple sequence alignments with millions of sequences.

Finally, it is not clear how to assess the quality or reliability of such large trees. Different methods gave very different topologies and large differences in likelihoodterm, and yet few of the differences were well-supported by the bootstrapterm. In fact, a topologyterm with relatively poor likelihoodterm could have relatively good agreement with the best tree. This could indicate that higher-likelihood trees contain many improvements, but that few of the individual improvements are statistically significant. This is expected if there is limited phylogenetic signal. Alternatively, the bootstrapterm could be too conservative. Local support valuesterm do suggest a greater number of significant differences ( Table 3 ), but local support valuesterm may be biased upwards because they do not consider all of the alternate topologies. Further study of these questions is needed.

Materials and Methods Minimum-Evolutionterm “Linear” Subtree-Pruning-Regraftingterm To reduce the number of SPR moves considered from O( ) to O( ), FastTree does just two rounds of “linear SPRs.” For each node, FastTree does an exhaustive search for moves up to length two. It extends each of these moves up to a length of 10 along the best choice at each point along the way.

As suggested by Richard Desper and Olivier Gascuel, FastTree treats each potential SPR move as a sequence of NNIs. The change in tree length for the SPR move is then just the sum of the changes due to the NNIs, much as .

The change in tree length for an NNI from to , where A, B, C, and D may be subtrees rather than sequences, is estimated by . In FastMEterm, which introduced balanced minimum evolutionterm·def [12] , is a topologically weighted average of distances between the members of A and C. In contrast, in FastTree, is the log-corrected distance between the profiles for the subtrees A and C, and the profile of a subtree is derived from that of its children by . (Although FastTree 1 used weighted joins, as in BIONJ, FastTree 2term uses unweighted joins because they are faster, and the slight effect on accuracy is erased by the ML NNIs.) For nucleotide sequences, the log correction is the Jukes-Cantor correction , where is average dissimilarity of positions across profiles. For amino acid sequences, FastTree uses an empirical log correction similar to that of scoredist [32] , , where is based on an amino acid dissimilarity matrix derived from the BLOSUM45 similarity matrix.

In FastMEterm, the above formula for the change in tree length is exactly correct because the changes in other branch lengthsterm in the tree can be expressed as combinations of distances that cancel each other out [26] . In FastTree, however, the formula for the change in tree length is an approximation, because the log-corrected distances do not cancel in this way. Nevertheless, FastTree with NNIs and FastMEterm with NNIs give very similar results [2] , and computing the exact change in total tree length does not improve the accuracy of FastTree's SPRs (data not shown).

The Maximum-Likelihood Phase The key data structures for the maximum likelihoodterm·def phase are the tree topologyterm, the branch lengthsterm, and the posterior distributionsterm for each internal node. (FastTree stores the tree with a trifurcation at the root, but the placement of the root is not biologically meaningful and does not affect the likelihoodterm [5] .) The posterior distributionterm·def for an internal node describes the state of the corresponding ancestor, given the branch lengthsterm and the sequences beneath it. For example, for nucleotide data, it stores the probability that a given site was an A, C, G, or T. FastTree stores posterior distributionsterm for internal nodes (not for the root), and they require O( ) space each, where is the alignment's length and is the number of characters in the alphabet.

The key primitive operations are (1) to compute the joint likelihoodterm of two posterior distributionsterm, given the length between them, and (2) to compute the posterior distributionterm·def of a parent node given the posterior distributionsterm of its two children and their two branch lengthsterm. These suffice to compute the likelihoodterm of the tree [5] : for example, the likelihoodterm of the tree (A,B,(C,D)) is where AB and CD are posterior distributionsterm.

At the beginning of the ML phase, we have a minimum-evolutionterm topologyterm and branch lengthsterm. The steps for the maximum-likelihood phase are: Compute an approximate posterior distributionterm·def for each node, using the weighted averages of its children. Although the initial posterior distributionsterm are approximate, all future changes to the topologyterm or to the branch lengthsterm will update the posterior distributionsterm to their exact values.

Optimize all branch lengthsterm for one round, using a simplified model with no parameters (without CATterm, and with Jukes-Cantor instead of GTR if GTR was requested).

Perform one round of ML NNIs, using the simplified model.

If the GTR modelterm is being used, optimize the nucleotide transition rate parameters, switch from Jukes Cantorterm to the GTR modelterm and recompute posterior distributionsterm, and optimize all branch lengthsterm for one round with the new model.

If the CATterm model is being used, estimate rate categoriesterm for each site, recompute posterior distributionsterm, and optimize all branch lengthsterm for one round with the new model.

Perform additional rounds of ML NNIs, with subtree skippingterm and the star topologyterm test.

Perform a final round of ML NNIs without subtree skippingterm or the star topologyterm test.

Optimize all branch lengthsterm for one round.

Compute SH-like local supportterm·def values.

A round of ML NNIs During each round of NNIs, FastTree visits each node before it visits its parents (depth-first post-order traversal). At each node, it compares the likelihoodterm of the trees , , and , where A and B are its children, C is its sibling, and D is the rest of the tree. During this process, FastTree uses the posterior distributionsterm (or sequences) for A, B, and C and an “up-posteriorterm” D, which represents the rest of the tree. More precisely, the up-posteriorterm D is the posterior distributionterm·def of the node's parent N, given all of the nodes that are not children of N (see Figure 3 ). These up-posteriorsterm can be thought of as a way to temporarily reroot the tree at the current location. In particular, the likelihoodterm of the tree can be computed from the posteriors A, B, C, and D.

Traversing a tree with up-posteriorsterm. FastTree optimizes the tree near node N by analyzing the posterior distributionsterm for subtrees A, B, and C, as well as the “up-posteriorterm” D.

The up-posteriorterm for a node can be computed from its parent's up-posteriorterm and its sibling's posterior distributionterm·def. FastTree only stores these up-posteriorsterm for the path to the root from its current location in the tree, so they take O( ) space, where is the maximum depth of the tree. Because FastTree always visits children before their parents, the posterior and up-posteriorterm distributions it uses are up to date, even as the topologyterm changes.

When it visits each node, for each of the three alternate topologies around the node, FastTree optimizes the branch lengthsterm to maximize the likelihoodterm. For the topologyterm , the five initial branch lengthsterm are set from the current tree. For the other topologies, the branch lengthsterm to A, B, C, and D are maintained, as is the internal branch lengthterm. Given a quartet (say ), FastTree first optimizes the branch lengthterm between AB and CD, and then the branch lengthterm leading to A, B, C, and D. FastTree optimizes each branch lengthterm to an accuracy of 0.0001 or 0.1%, whichever is greater. These five optimizations define a round of optimization for the quartet. Within a round of optimization, FastTree reuses some of the internal posterior distributionsterm: it needs posterior distributionsterm for AB and CD so that it can optimize the branch lengthterm between AB and CD, and then it needs posterior distributionsterm for BCD, ACD, the new posterior distributionterm·def for AB given the new branch lengthsterm to A and B, and finally ABD and ABC.

By default, FastTree optimizes the branch lengthsterm within all three quartet topologies for one round. Any topologyterm that is significantly (5 log-likelihood units) worse than the current topologyterm is abandoned after the first round. If more than one topologyterm remains, then the remaining topologies are optimized for another round. After the rounds of optimization are complete, FastTree updates the topologyterm if necessary. In either case, it updates the branch lengthsterm to the re-optimized values and recomputes the posterior distributionterm·def for the node.

A difference of 5 in log likelihoodterm may seem like a small difference, so that the heuristicterm might miss a good change to the topologyterm. However, optimization of branch lengthsterm after the first round usually leads to small improvements in the log likelihoodterm. For example, if we analyze 40 randomly selected 16S rRNAs with FastTree and the GTR+CATterm model, and we increase the rounds of branch length optimizationterm·var to 4 (-mlacc 4), then the average improvement for any NNI is just 1.1 log-likelihood units in the second round of branch lengthterm optimization and just 0.04 in rounds 3 and 4 combined. To put these numbers in perspective, differences in log-likelihood of less than 2 are not statistically significant ( , likelihoodterm ratio test), and NNIs with much larger changes in likelihoodterm are common. For the simulated alignments with 5,000 protein sequencesterm, always optimizing for two rounds improved accuracy by a negligible amount (0.03%) and increased the running timeterm·var by 23%.

Optimizing model parameters After the first round of NNIs, FastTree optimizes any parameters in the model. First, if the GTR modelterm is being used, there are six relative rates to optimize, one for each nucleotide conversion. (The stationary distribution for the transition matrix is set to the empirical frequency of the four nucleotides.) FastTree optimizes the likelihoodterm of the tree (with fixed branch lengthsterm and topologyterm) by numerically optimizing each of the six parameters in the model in turn. With each change in the model, it recomputes all posterior distributionsterm. It then optimizes the six parameters a second time. This does not fully optimize the model parameters, but it gives acceptable results ( Table 2 ).

Second, unless the -nocat option is set, FastTree estimates the rate of evolutionterm at each site. Given the desired number of categories of relative rates , FastTree selects values that are logarithmically spaced between and . By default, , and the relative rates range from 0.05 to 20. For each of these relative rates, FastTree recomputes all posterior distributionsterm and recalculates the log likelihoodterm of the tree at each site. FastTree then uses a Bayesian approachterm to select which rate to use at each site: FastTree maximizes , where is a gamma-distributed prior. To avoid overfitting, we made the prior more peaked than real rate variation in alignments: the prior has a shape parameter of 3, a scale parameter of 1/3, and a mean of 1. After choosing the rate categoriesterm, FastTree scales the rates so that the average rate across all sites is 1.0.

We confirmed that the Bayesian approachterm to setting the rate categoriesterm prevents overfitting on small alignments. For example, on simulated protein alignments with just 10 sequences (from [2] ), adding the CATterm model improves FastTree's accuracy from 76.2% to 78.0%. (For comparison, PhyMLterm without or SPRs was 74.4% accurate [2] .) Conversely, on nucleotide simulations with 24 sequences that (unrealistically) do not contain any rate variation across sites (the fast-evolving alignments of [12] ), the CATterm model only reduces accuracy slightly, from 93.6% to 93.4%. (For comparison, PhyMLterm without or SPRs was 93.6% accurate [2] .) Completing the ML NNIs In later rounds of NNIs, FastTree uses the more accurate model and it uses two additional heuristicsterm “subtree skippingterm” and the “star topologyterm test,” which are described below. As discussed in the Results , these heuristicsterm have little effect on accuracy.

If no NNI leads to an improvement of more than 0.1 in the likelihoodterm of any quartet, then FastTree considers the NNIs to have converged. FastTree repeats rounds of NNIs until convergenceterm·def, up to a limit of rounds, which takes O( ) time. This is the slowest part of the ML phase. The limit on rounds ensures a predictable running timeterm·var, but FastTree usually converges before reaching the limit, even for huge alignments such as 237,882 16S rRNA sequences. We chose a limit so that a misplaced subtree could move all the way across a (roughly balanced) tree, and the factor of 2 is an arbitrary safety factor.

After convergenceterm·def, FastTree does one final round of ML NNIs with the subtree skippingterm and the star topologyterm test turned off, as in the first round. We view this as a safety valve for the heuristicsterm. Finally, FastTree does a final round of optimizing the branch lengthsterm and computes the SH-like local supportsterm·def.

Subtree skippingterm The intuition behind subtree skippingterm is that if a subtree has not changed during recent rounds of NNIs, then further attempts to optimize the subtree will be fruitless. Specifically, during ML NNIs, FastTree does not traverse into subtrees that have not seen any significant improvement in likelihoodterm (0.1 log likelihoodterm units) in either of the previous two rounds. Before skipping a subtree, FastTree also checks that none of the nodes adjacent to the parent node were affected by a significantly improving NNI in the previous round. The “subtree skippingterm” heuristicterm typically gives a 3-fold speedup, making it the most important of FastTree's ML heuristicsterm. Subtree skippingterm might be useful for SPR moves as well.

Star topologyterm test If the current topologyterm (A,B,(C,D)) is much better than the star topologyterm (A,B,C,D) then an NNI is unlikely to give an improvement. Specifically, if the current topologyterm is significantly (5 log-likelihood units) more likely than the star topologyterm (after optimizing the internal branch lengthterm), then FastTree does not optimize the other branch lengthsterm or consider the two alternate topologies. However, FastTree only uses this heuristicterm if the node that was unchanged in the last round of NNIs. To approximate the likelihoodterm of the star topologyterm, FastTree uses the likelihoodterm with the minimal internal branch lengthterm of 0.0001.

Branch lengthsterm To optimize all branch lengthsterm in the tree at the beginning and end of the ML phase and after optimizing the model parameters, FastTree again uses post-order traversal. At each node, it considers a three-node star topologyterm on the node's children and parent, using the posterior distributionsterm for the two children and the up-posteriorterm for itself. (At the root, it uses all three children instead.) It numerically optimizes these three branch lengthsterm in series for two rounds.

SH-like local supportsterm·def For each node, the local support is derived from the per-site likelihoodsterm for the current topologyterm and the two alternate (NNI) topologies. For the current topologyterm, FastTree uses the current (already optimized) branch lengthsterm. For the alternate topologies, FastTree optimizes branch lengthsterm for the quartets, as during the NNIs, for up to two rounds. Given the per-site likelihoodsterm for the three topologies, FastTree uses the SH testterm with 1,000 bootstrapterm replicates to estimate the confidence in the given split [19] . If there are poorly resolved nodes nearby, then the support values should be interpreted cautiously, because a high-likelihood alternate topologyterm might not have been considered.

Low-level optimization of likelihoodterm computations Whereas RAxMLterm stores likelihood vectorsterm (that is, the joint likelihoodterm of a subtree and of a given character at an internal node), FastTree stores posterior distributionsterm, which are normalized so that each site's values sum to 1. This may improve numerical stability for huge alignments. To reduce memoryterm·var usage, FastTree stores these vectors in single-precision floating pointterm. Log-likelihoods for the tree or for specific sites are stored with double precision.

Similar to RAxMLterm, FastTree stores the posterior distributionsterm in a rotated formterm, multiplied by the eigen-matrix of the transition matrix. (For the Jukes Cantorterm model, this is not necessary.) This reduces the time for likelihoodterm computations from O( ) per site to O( ), while leaving the cost of computing the posterior distributionterm·def at O( ) per site (but with a higher constant factor).

While computing the joint likelihoodterm for a pair of posterior distributionsterm, FastTree avoids performing a logarithm at every site by operating on likelihoodsterm instead of log likelihoodsterm. To prevent numerical underflow, FastTree rescales the likelihoodterm by a constant when necessary. It updates a separate (log-likelihood-based) counter whenever it does this. Similarly, when computing the tree's likelihoodterm at each site, for example while optimizing the rate categoriesterm, FastTree rescales each site's likelihoodterm if necessary after visiting each node.

FastTree uses SSE2 instructionsterm, a special feature of recent CPUs from Intel and AMD, to operate on 4 single-precision floating pointterm values with one instruction. This speeds up computations for protein alignments by up to 50% (data not shown).

Numerical optimization To find the parameters that optimize the likelihoodterm, FastTree uses Brent's methodterm, a numerical method that iteratively halves the interval it is searching within ( http://en.wikipedia.org/wiki/Brents_method ). Because Brent's methodterm only operates in one dimension, FastTree optimizes different parameters in turn, and then repeats the rounds of optimization (for example, it optimizes the first branch lengthterm, then the second, then the third, then repeats).

FastTree estimates the initial interval to search within from the initial guess (e.g., the previous length of the branch) and alternate values and . If is below the minimum value, it uses the minimum, , and instead. If the initial guess does not bracket the minimum (that is, the middle value is not better than the two endpoints), then FastTree expands the search interval until it does. However, the small interval is usually adequate. FastTree also terminates optimization if the parameter changes by a small amount or by a small proportion. Together, these modifications eliminate about a third of the evaluations of the likelihoodterm.

Biological and Simulated Alignments The simulated protein alignments and the genuine COG alignments were described previously [2] . The 16S alignment with 237,882 distinct sequences was taken from GreenGenes [33] ( http://greengenes.lbl.gov ). The 16S alignment with 15,011 distinct “families” is a non-redundant subset of these sequences ( identical). 16S alignments with 500 sequences are also non-redundant random subsets ( identical). Other large 16S alignments are from [11] .

For the 16S-like simulations with 78,132 distinct sequences, we used a maximum-likelihood tree inferred from a non-redundant aligned subset of the full set of 16S sequences ( % identity) by an earlier version of FastTree (1.9) with the Jukes-Cantor model (no CATterm). To ensure that the simulated trees were resolvable, which facilitates comparison of methods (but inflates the accuracy of all methods), branch lengthsterm of less than 0.001 were replaced with values of 0.001, which corresponds to roughly one substitution across the internal branch, as the 16S alignment has 1,287 positions. Evolutionary rates for each site were randomly selected from 16 rate categoriesterm according to a gamma distributionterm with a coefficient of variation of 0.7. Given the tree and the rates, sequences were simulated with Rose [34] under the HKY model and no transition bias. To allow Rose to handle branch lengthsterm of less than 1%, we set “MeanSubstitution = 0.00134” and multiplied the branch lengthsterm by 1,000.

Software Used We used FastTree 2term.0.0. We used the July 6 2009 release of the PhyMLterm 3.0 source code and modified BL_MIN from 1.e-10 to 1.e-8 to overcome numerical problems with some of the simulated protein alignments, as suggested by Stepháne Guindon. FastMEterm 2.06 was provided by Olivier Gascuel. RAxMLterm 7.0.4 and 7.2.1 were obtained from the author's web sites. RAxMLterm 7.2.1 was compiled with SSE instructions. NINJA was provided by Travis Wheeler and is available at http://nimbletwist.com/software/ninja/ . BIONJ was obtained from http://www.lirmm.fr/~w3ifa/MAAS/BIONJ/BIONJ.c . BIONJ was run with maximum-likelihood distances obtained with phylip's protdist ( http://evolution.genetics.washington.edu/phylip.htm ) and the JTT model (no gamma). Log-corrected distances were obtained with FastTree and the -makematrix option.

Supporting Information Branch lengthsterm for an alignment of 200 16S rRNA sequences vary systematically with the Γ approximation used. The CATterm lengths are from FastTree, and all Γ branch lengthsterm are from PhyMLterm with FastTree's topologyterm and with optimized shape parameters. The top panel shows that branch lengthsterm from the various models have a roughly linear relationship with each other, but they have different scales. The bottom panel shows how the total length of the tree varies with the number of categories (note log χ axis). The “Use Median” lengths are from running PhyMLterm with –use_median, which uses the median of each region, rather than the mean, to approximate the gamma distributionterm. The “Corrected” lengths are the “Use Median” lengths multiplied by the average posterior rates, which can be obtained by running PhyMLterm with –print_site_lnl (thanks to Stepháne Guindon for pointing this out). The corrected lengths converge to the correct value much more quickly than the other rates. The “CATterm/Gamma” tree length, from FastTree 2term.1 with -gamma, is also reasonably accurate. With this option, FastTree 2term.1 optimizes the Γ 20 likelihoodterm with a shape parameter and a rescaling parameter, using the site likelihoodsterm from FastTree's 20 relative rates and branch lengthsterm that were optimized under the CATterm model.

(0.13 MB PS) Click here for additional data file.

Total Splits or Strongly Supportedterm Splits that Disagree with RAxMLterm's Final Tree, versus Time. The 16S tree has 4,111 splits and the COG trees have 2,497 splits each. All values for the COG trees are averages over the 7 COGs.

(0.02 MB PS) Click here for additional data file.

Times and likelihoodsterm for large 16S rRNA alignments (0.02 MB PDF) Click here for additional data file.