Constrained incremental tree building: new absolute fast converging phylogeny estimation methods with improved scalability and accuracy

The inference of phylogenies from molecular sequence data is generally approached as a statistical estimation problem, in which a model tree (equipped with a model of sequence evolution) is assumed to have generated the observed data, and the properties of the statistical model are then used to infer the tree. Various statistical approaches can be applied for this estimation, including maximum likelihood, Bayesian techniques, and methods that operate by computing a distance matrix and then computing the tree from the distance matrix.

Many stochastic sequence evolution models have been developed, starting with the Cavender–Farris–Neyman [1, 2] symmetric two-state model (referred to henceforth as “CFN”) and including increasingly complex molecular sequence evolution models (with four states for DNA, 20 states for amino acids, and 64 states for codon sequences). However, typically the theory that can be established under the CFN model can also be established under the more complex molecular sequence evolution models used in phylogeny estimation.

Under standard sequence evolution models, many methods are known to be statistically consistent (meaning that they will provably converge to the true tree as the sequence lengths increase), including maximum likelihood [3] and many distance-based methods [4, 5]. However, with the increasing availability of sequence datasets (e.g., the SILVA database has several million RNA sequences [6]), the need for phylogeny estimation methods that can be highly accurate on ultra-large datasets has increased. Therefore, a key challenge is to have methods that can scale to large datasets while maintaining good accuracy. Hence, the most desirable methods are those that run in polynomial time and that recover the true tree with high probability from short sequences (i.e., sequences that do not have very many sites). In this respect, methods that are “absolute fast converging” [7–9] (i.e., methods that recover the true tree with high probability from polynomial length sequences) are the most promising.

There are several methods that have been established to be absolute fast converging (AFC) under the CFN model, including maximum likelihood [3] (if solved exactly) and various distance-based methods [7–14]. Some of these algorithms achieve a poly-logarithmic sample complexity but require a balanced model tree and an upper bound on g, the maximum edge weight (defining the expected number of changes of a random site) in the CFN model. Specifically, methods based on reconstruction of ancestral sequences provide the best sample complexity bounds but cannot handle the case where g is larger than what is known as the Kesten–Stigum threshold, which is \(\ln (\sqrt{2})\) [15] for the CFN model. The “short quartet” methods were the earliest AFC methods (which are AFC in the regime where g is unbounded), but these are designed to either return the true tree or else fail to return anything [7–9].

Two of the fastest AFC methods are the Harmonic Greedy Triplets (HGT + FP) method by Csűrös [16] and a method developed by King et al. [17]; these methods use \(O(n^2)\) time and are based on quartet trees, and have the desirable property that they always return a tree for every input. Another good approach is \(\texttt {DCM}_{NJ},\) which uses a divide-and-conquer technique [9]. In the first phase, \(O(n^2)\) trees are computed, each based on dividing the sequence dataset into overlapping subsets, constructing trees on each subset using the polynomial time distance-based method neighbor joining (NJ) [18], and then combining the subset trees using a supertree method. This approach differs from other AFC methods in its use of neighbor joining to construct trees on subsets, whereas the other AFC methods in essence construct the tree by independently constructing quartet trees, and then assembling the quartet trees together using a quartet amalgamation method.

Very few AFC methods have been implemented; however, a study [10] comparing \(\texttt {HGT}+\texttt {FP}\) [16] and \(\texttt {DCM}_{NJ}\) [9] showed that \(\texttt {DCM}_{NJ}\) had better accuracy. Since the theory does not predict this, the results on simulated data suggest that the use of \(\texttt {NJ}\) to construct trees on subsets and then combine the trees using a supertree method may be empirically advantageous compared to methods that combine quartet trees that are estimated independently. Unfortunately, the reliance on a supertree method to combine subset trees means that \(\texttt {DCM}_{NJ}\) is unlikely to scale to ultra-large datasets, because no current supertree method has shown the ability to maintain good accuracy and reasonable running times on large datasets [19].

The purpose of this paper is to describe a new polynomial time AFC phylogeny estimation method that should improve on \(\texttt {DCM}_{NJ}{:}\) it is designed to have comparable accuracy to \(\texttt {DCM}_{NJ}\) but also to be able to analyze ultra-large datasets (i.e., more than 100,000 sequences). The basic approach of this new method is similar to \(\texttt {DCM}_{NJ}\) in that it uses divide-and-conquer, applies neighbor joining to subsets of the sequence dataset, and then merges the subtrees together. However, it differs from \(\texttt {DCM}_{NJ}\) in a few important ways, which we describe below. Most importantly, it divides the taxon set into disjoint subsets and then merges the subset trees without relying on any supertree method; thus, it avoids the challenge of relying on existing supertree methods, none of which are likely to scale to large datasets. We present the results here initially for the CFN model, and then extend the results for the Generalized Time Reversible (GTR) model [20]. Our arguments for correctness are very similar to those in [16, 17].

We begin by describing the \(\texttt {INC}\) method in the unconstrained condition (i.e., when there are no constraint trees), and prove that it is AFC.

High-level description of \(\texttt {INC}\)

The input to \(\texttt {INC}\) is a dissimilarity matrix d. We present a high-level description of how tree \(t = \texttt {INC}(d)\) is computed.

• Find the insertion ordering \(\sigma = x_1, x_2, \ldots , x_n.\)

• Initialize t as the three-leaf tree on the first three taxa in the ordering.

• For \(i=4\) up to n.

  • Determine the set of valid quartets to use for placing \(x_i\) into tree t.

  • Compute quartet trees for each valid quartet, and let them vote for where to place \(x_i.\)

  • Pick the edge e in t that has the most votes; if this is a tie, pick an edge at random from the set of edges with the most votes.

  • Insert \(x_i\) into edge e.

• Return the resulting tree t.

Thus, at a high-level, the \(\texttt {INC}\) algorithm operates by greedily growing a tree t based on a computed sequence addition ordering. Yet, many of the details of the algorithm are unspecified (e.g., we do not say how we calculate the insertion ordering, how we determine the set of valid quartets, how we compute trees on valid quartets, and how the quartet trees vote). Below, we provide details for each of the steps for this algorithm.

Computing the sequence addition ordering

We begin by constructing a minimum spanning tree S of \(TG(d,\infty ),\) where d is the input dissimilarity matrix. Once the spanning tree S is computed, we choose an arbitrary leaf in S to be the starting vertex in the ordering. Then we order the vertices according to order of traversal in a BFS (or DFS) from the starting vertex. See Fig. 1 for an example.

How we compute the growing tree t

We initialize our tree t by choosing the first three taxa according to our insertion ordering and forming the three-leaf tree with an internal node that connects to all three taxa. In order to define how we insert the remaining taxa into t, we need to formally define the “valid quartets” and how we compute quartet trees for valid quartets.

Definition 4

Let \(q_0\) be the maximum weight of an edge in S and let \(q = 8q_0,\) A quartet of leaves is valid if its maximum pairwise distance (i.e., diameter) is at most q. The quartet tree for the valid quartet is computed using the Four Point Method.

As we will show later, the restriction to just the “valid quartets” allows us to develop a tree construction method that runs in polynomial time and that is AFC. Note also that restricting the diameter of the valid quartets to small values has a mixed effect: if the maximum permitted diameter is too small then even correct quartet trees will not be sufficient to reconstruct the tree, but if the maximum permitted diameter is too large then some computed quartet trees are more likely to be incorrect. We show that this setting for the maximum diameter q to be \(8q_0\) is sufficient to allow us to prove that the algorithm is AFC. However, we did not try to optimize this constant, and hence the choice of the constant 8 is likely not optimal (i.e., smaller constants might give better theoretical results).

There may be many valid quartets, but we only need to examine a linear number of these, as we now show. Suppose we wish to add a vertex x into t. Given an internal node u of t, because t is binary the removal of u splits t into three non-empty components which we will refer to as \(t_1,t_2,t_3\) (with internal nodes included). Because the leaves of t form a connected induced subtree of S, we can find taxa \(u_i\) in \(t_i\) and \(v_i \in V(t)\setminus t_i\) such that \((u_i,v_i)\) is an edge in S, for \(i=1,2,3.\) We will associate the triplet \(u_1,u_2,u_3\) to node u, so that the addition of leaf x defines a quartet, and we can check to see if the quartet is valid. We summarize this discussion with the following definition.

Definition 5

Let u be an internal node of a binary tree t and let \(u_1,u_2,u_3\) be leaves in the three components of t upon removing u such that there exists \(v_i\) with \((u_i, v_i)\) an edge in S where \(v_i\) is not in the same component as \(u_i.\) Then, for some choice of \(q \in \mathbb {R}^+,\) a quartet query on \(\{u_1,u_2,u_3,x\}\) is q-valid iff the d-diameter (maximum pairwise distance with respect to the input matrix d) of \(\{u_1,u_2,u_3,x\}\) is less than q. We use the term valid when q is clear from context.

To determine how to place leaf x into the growing tree t, we compute a tree for each valid quartet query (which includes x) using the Four Point Method. For example, if the tree computed on valid quartet query \(u_1,u_2,u_3,x\) (with \(u_i \in t_i,\) and \(t_1,t_2,t_3\) the three subtrees off internal node u) returns \(xu_1|u_2 u_3,\) then this implies that x should be placed in the subtree of t induced on \(t_1 \cup {u}\) hence, each edge in that subtree will receive a vote (Fig. 2).

Thus, each valid quartet adds a vote to each edge in some non-empty subtree of the tree. We define the support of an edge in the tree t to be the number of valid quartet trees that voted for that edge. We then choose an edge e in the tree t that has the largest total support (and if there is a tie, we pick an edge e at random with the maximum support). We then subdivide e and then make x adjacent to the node created. See Fig. 3 for an illustration of how this voting procedure operates.

Furthermore, as we will shortly show, when the sequences are long enough, then there will be a unique edge on which all queries agree, and so the algorithm will correctly add x into t (in such a way that it agrees with the model tree T), and so inductively the greedy algorithm will construct the model tree.

We show how \(\texttt {INC}\) can be modified to take an arbitrary set of disjoint constraint trees. We present several variants of this method: one that optimizes speed while still guaranteeing that the resultant algorithm is AFC, and other variants that are designed for improved empirical accuracy but with an increase in running time and potentially a loss of the AFC property.

Constrained \(\texttt {INC}\)

The input to the constrained version of \(\texttt {INC}\) includes a set of leaf-disjoint trees. Therefore, the constraint set is compatible (i.e., a compatibility tree exists). We will describe a straightforward modification to \(\texttt {INC}\) so that it never violates the topological information in the constraint trees, by which we mean only that the final output tree is required to induce each constraint tree, when restricted to that specific set of leaves. Thus, the constraint trees are not required to be clades in the output tree.

Before proceeding, we define the induced tree topology, \(t^{A},\) on a tree t and a subset of leaves A as follows. Consider the minimal subtree of t that contains the leaves A, and then suppress all nodes of degree two (i.e., replace all maximal paths of degree two nodes by a single edge). The endpoints of an edge \(e \in t^{A}\) correspond to vertices in t and so e corresponds to a (possibly length 1) path in t; we denote the corresponding endpoints of e as \(e_t(u)\) and \(e_t(v).\) For an induced tree \(t^{A},\) the component in t corresponding to an edge \(e \in t^A\) is the set of edges and vertices that can be reached by a walk starting from \(e_t(u)\) and ending at \(e_t(v)\) without using \(e_t(u)\) or \(e_t(v)\) multiple times. Note this will be a subgraph of t that includes \(e_t(u)\) and \(e_t(v)\) as leaves. See Fig. 4.

When working with constraint trees, we alter the insertion algorithm to account for the constraints. Recall that the constraint trees are on disjoint sets of taxa. Thus, when inserting a taxon, we identify the corresponding constraint tree \(t_c\) that includes that new taxon. Let A be the set of taxa shared between \(t_c\) and t and note that \(t_c^{A} = t^{A}\) since the trees are consistent. Thus, in \(t_c^{A\cup \{x\}},\) the taxon x can be viewed as being attached to some edge e in \(t^{A}.\) The eligible edges are in the component in t corresponding to e. We then simply run \(\texttt {INC}\) on the component to find the desired edge of insertion.

In this section, we extend our AFC convergence guarantees to the generalized time reversible (GTR) Markov model [20], which is the most commonly used site evolution model used in phylogenetics. We will show that \(\texttt {INC}\) is AFC for the GTR model, drawing on established techniques and theory which we now present and summarize.

The GTR model can be seen as a special case of the general Markov model [25], which we first describe. Under the general Markov model on \(m \ge 2\) states, we have a rooted binary tree T with some distribution of states \(\pi > 0\) at the root of the tree and m-by-m stochastic transition matrices M(e) for each edge e of T. Starting with a random sequence drawn i.i.d. from \(\pi,\) each site of the sequence evolves i.i.d. down the tree, according to the transition matrices specified by the model.

Note that \(\det (M(e))\) takes the values 1 or – 1 precisely if M(e) is a permutation matrix. Also, for the CFN model, \(det(M(e)) = 1 - 2p(e),\) where p(e) is the substitution probability on edge e; thus we know that \(\det (M(e))>0\) and \(\det (M(e)) \rightarrow 0\) as \(p(e) \rightarrow 0.5,\) which correctly supports the intuition that information loss is large when p(e) is close to 0.5. In general, \((1/2)[1 -\det (M(e))]\) plays the role of p(e) in the general model and \(-\log (\det (M(e)))\) is the corresponding w(e).

Thus, the natural extension of \(CFN_{f,g}\) model to the GTR model, \(GTR_{f,g},\) is to enforce \(f \le -\log (\det (M(e))) \le g\) for all edges in all model trees in \(GM_{f,g}.\)

Now, we apply INC with the extended version of the empirical distance matrix and we define TG(d, q), \(q_0,\) \(\epsilon (q)\) analogously for our generalized distance measure. Note that by our Markov property and time-reversibility, if P is the unique path between leaves i, j in T, then the true distance matrix is additive and satisfies \(D_{ij} = \sum _{e \in P} - \log (\det (M(e)))\) up to a constant shift [8]. Thus, if \(\epsilon (q) < f/2,\) then applying a quartet query via the Four Point Method will return the accurate quartet tree. Therefore, the claim that \(\texttt {INC}\) is AFC for the GTR model follows immediately after we establish a slightly extended concentration bound.

This paper presents a novel algorithmic technique for constructing phylogenetic trees, which allows statistically consistent and highly accurate methods to be used on subsets of the taxa in a divide-and-conquer framework. We proved that several of these variants are absolute fast converging (AFC) under the CFN and GTR sequence evolution models, and that some of these methods achieve \(O(n^2)\) time and space (where n is the number of sequences) once the dissimilarity matrix relating the sequences is computed. More generally, tree estimation methods could be scaled to large datasets using this approach, provided that it is possible to compute a matrix of pairwise distances that is guaranteed to converge to an additive matrix as the amount of data increases.

Many of the ideas in the algorithmic design of \(\texttt {INC}\) are derived from prior algorithms for similar problems. For example, the node query insertion procedure was explored earlier in [26, 27] to produce a fast \(O(n \log n)\) tree-growing algorithm, but analyzed under a different model of quartet tree error than we use here; furthermore, a direct implementation of their model does not produce an AFC algorithm. The idea of using a distance matrix to combine disjoint trees was originally presented in [28] in the context of multi-locus species tree estimation; that technique (called NJMerge) is a modification of the neighbor joining method [18] to allow for constraint trees. Since neighbor joining is not AFC [29], it seems unlikely that NJMerge can be used to produce an AFC method.

The goal of this work was improved empirical performance, compared to prior AFC methods. One of the interesting aspects of \(\texttt {INC}\) is the potential to use it with disjoint constraint trees. As noted in this paper, maximum likelihood (if solved exactly) is AFC under the standard sequence evolution models, and although it is NP-hard to solve exactly there are many seemingly good heuristics for maximum likelihood (e.g., RAxML [30]). One possible technique that could result in excellent empirical accuracy and scalability is to use a maximum likelihood heuristic to compute constraint trees, and then combine the constraint trees using constrained \(\texttt {INC}.\) The value of such an approach, however, depends on the decomposition strategy, as not all decompositions will produce AFC methods. Consider, for example, the impact of selecting quartets from a large phylogeny: if these quartets are selected arbitrarily, it is possible to produce quartets whose true trees fall in the Felsenstein Zone, where methods such as maximum parsimony are not statistically consistent and maximum likelihood—although consistent—will tend to have poor accuracy unless very long sequences are available [31, 32]. Therefore, the use of constrained \(\texttt {INC}\) with maximum likelihood, or any other tree estimation method, will only be AFC for some decomposition strategies (and not for arbitrary ones), and may not provide accuracy advantages except when combined with carefully designed decomposition strategies.

Thus, additional work is needed to develop provably AFC divide-and-conquer strategies that enable this kind of approach to be used to the greatest empirical advantage on large challenging phylogenetic datasets.