Terrace Aware Data Structure for Phylogenomic Inference from Supermatrices

Partial and Full Phylogenetic Terraces

Let $k⩾2$ denote the number of genes, loci, or codon positions in a supermatrix. In the following we use “partition” to generally refer to any subset of genomic positions. Denote by $Y_1,Y_2,\dots ,Y_k$ the species sets for the $k$ partitions and $X=Y_1\cup Y_2\cup \dots \cup Y_k$ the set of all species. The $S_1,S_2,\dots ,S_k$ denote the corresponding alignments and $S$ is the concatenated alignment (supermatrix) of $S_1,S_2,\dots ,S_k$. Stretches of gaps are added to $S$ if a species has no sequence for some partition (i.e., when $Y_i\subset X)$. Only these gaps are referred to as missing data.

For a species tree $T$ and a partition $Y_i$, the associated induced partition tree, denoted $T|Y_i$, is the tree obtained from $T$ by pruning species with no sequence for partition $Y_i$ (i.e., missing sequences). Hence, for every species tree there is a corresponding set of $k$ induced partition trees.

If for two species trees $T_1$ and $T_2$ there exists a set of indices $J\subseteq \{1,\dots ,k\}$ such that $\forall j\in J$ the corresponding induced partition trees $T_1|Y_j$ and $T_2|Y_j$ are identical then $T_1$ and $T_2$ belong to one partial terrace (Chernomor et al. 2015). In this context, a phylogenetic terrace coined in Sanderson et al. (2011) is a special case when $J=\{1,\dots ,k\}$. For clarity we call this case a full terrace. $\forall j\in J$ scores (likelihood or parsimony) of $T_1|Y_j$ and $T_2|Y_j$ are equal. Therefore, if during the tree search we identify a full terrace, we need to compute the score of $T_1|Y_j$ or $T_2|Y_j$ only once $\forall j\in J$.

How to Identify Partial Terraces during the Tree Search

When searching for the optimal species tree we move from one species tree $T$ to another $T_{NEW}$ by means of some topological rearrangements. In phylogenetic software the most common topological rearrangements are Nearest Neighbor Interchange (NNI), SPR, and Tree Bisection and Reconnection (TBR) (Felsenstein 2004).

We now illustrate the necessary condition presented in Chernomor et al. (2015) for the NNI to change the topology of a given induced partition tree. Let $T$ be a species tree and $e$ an interior edge of $T$. We denote by $e_1,e_2,e_3,e_4$ edges adjacent to $e$ and by $A,B,C,D$ the taxon sets leading from them, respectively (Fig. 1). Now, let a new tree $T_{\text{NNI}}$ be obtained from $T$ via NNI. Then Proposition 1 (Chernomor et al. 2015) states that

For a partition with a taxon set $Y$ the topologies of $T|Y$ and $T_{\text{NNI}}|Y$ are different if and only if $Y$ has at least one representative taxon in each subset $A,B,C,D$.    (C.1)

Condition (C.1) simply means that for an NNI to change the topology of a partition tree $T|Y$ all $A\cap Y,B\cap Y,C\cap Y$, and $D\cap Y$ must be non-empty.

In the following we discuss the three partition models implemented in IQ-TREE and used to examine the performance of the terrace aware data structure proposed here.

Partition Models

Partition models allow for different evolutionary scenarios for each partition. More formally, each partition $Y_i$ is assumed to evolve under its own substitution model $M_i$ and the model parameters of each $M_i$ are optimized separately on the corresponding partition tree.

The log-likelihood of a species tree $T$ under a partition model $M$ is the sum of partition tree log-likelihoods

Here, the log-likelihood $\ell$ of the partition trees depends on the topology and the edge lengths of each $T|Y_i$. The relation of edge lengths between species tree and partition trees is what distinguishes partition models.

The most general EUL model, denoted by $M_{\text{EUL}}$, allows each partition tree $T|Y_i$ to have its own set of edge lengths that are optimized separately per partition tree. We denote a length of $e$ on $T|Y_i$ by ${\lambda }_i(e)$. The EUL model implies that the species tree $T$ has no defined edge lengths. However, one can display the edge lengths for $T$, for example, as the weighted average of corresponding edge lengths on partition trees.

In contrast to EUL, the more restrictive EL models assume a relationship between edge lengths of the species tree and all partition trees. The EL-proportional model, denoted by $M_{{\text{EL}}_{\text{prop}}}$, assumes one set of edge lengths, $\lambda (.)$, for the species tree $T$ and rescales the partition trees with specific positive rates $r_1,r_2,\dots ,r_k$ such that the weighted average rate is 1 (i.e., $\frac{{\displaystyle {\sum }_{i=1}^k{w}_i}{r}_i}{{\displaystyle {\sum }_{i=1}^k{w}_i}}=1$, where $w_i$ is the length of the partition alignment $S_i)$. The partition rates $r_i$ are optimized separately for each partition tree.

The second EL model, EL-equal, denoted by $M_{{\text{EL}}_{\text{equal}}}$, is a special case of $M_{{\text{EL}}_{\text{prop}}}$ with all the partition rates equal to 1 (i.e., $r_1=r_2=\cdots =r_k=1$). This simply means that the edge lengths of partition trees are equal to the lengths of corresponding edges on the species tree. In contrast to $M_{\text{EUL}}$, which optimizes the edge lengths per partition tree, EL models optimize edge lengths, $\lambda (.)$, on the species tree.

Map from the Species Tree onto Partition Trees

Let $E$ denote the edge set of $T$ and $E_i$ the edge set of $T|Y_i$. We represent each edge $e\in E$ by its split $e=A|B$, where $A$ and $B$ are disjoint complementary non-empty subsets of the leaf set $X$ with $|X|=n$. Let $ϵ$ denote the empty edge or no edge, then for every partition $Y_i$ we introduce the map

Basically, $f_i(e)$, if not equal to $ϵ$, is an edge on $T|Y_i$ corresponding to $e$.

The collection of all maps $F=\left\{f_1,\dots ,f_k\right\}$ together with the trees $\left\{T,T|Y_1,\dots ,T|Y_k\right\}$ forms the PTA data structure for partition model analyses.

An Efficient Algorithm for Building F

We now describe a dynamic programming algorithm to build $F$ for unrooted bifurcating trees. It first assigns $f_i$ for external edges and then proceeds toward internal edges of the tree. More specifically, let $e=\left\{x\right\}|X\backslash \{x\}\in E$ be an external edge then

Obviously, Equation (3) follows directly from Equation (2). Now, let $e$ be an internal edge with two adjacent edges $e_1,e_2$ (e.g., see Fig. 2a), where $f_i(e_1)$ and $f_i\left(e_2\right)$ were already computed, we assign $f_i(e)$ as follows:

where $e_i^{\ast }$ is an edge on $T|Y_i$ adjacent to $f_i(e_1)$ and $f_i(e_2)$ (Fig. 2b). Note that for an edge $e$ where all four neighboring edges have their maps computed, assigning $f_i(e)$ from either side of $e$ will have the same result. We provide the proof of correctness of Equation (4) in Appendix 1.

The algorithm recursively computes $f_i$ using Equations (3 and 4) by a single post-order tree traversal. Thus, the time complexity of computing $f_i$ is linear in the number of taxa.

The described algorithm builds $F$ in $O(nk)$ time, where $n$ and $k$ are the number of species and partitions, respectively. Note that $F$ has to be recomputed each time the topology of the species tree changed. This sounds expensive at the first glance, but we will now show how to efficiently update the PTA when an NNI is applied to the species tree.

Identifying Unchanged Partition Trees with PTA

For a species tree $T$ and a partition $Y_i$ condition (C.1) provides a check whether an NNI applied to $e$ changes the topology of $T|Y_i$. Using the map notation, (C.1) is equivalent to

For a partition with a taxon set $Y_i$ the topologies of $T|Y_i$ and $T_{NNI}|Y_i$ are different if and only if all $f_i(e_1),f_i(e_2),f_i(e_3),f_i(e_4)\ne ϵ$. (C.2)

Condition (C.2) follows directly from (C.1) and the definition of $f_i(.)$.

Therefore, when an NNI is applied to $e$, one updates each partition tree $T|Y_i$ using two rules:

1. If all $f_i\left(e_1\right),f_i\left(e_2\right),f_i\left(e_3\right),f_i\left(e_4\right)\ne ϵ$, then $T_{NNI}|Y_i$ will result from $T|Y_i$ by swapping the corresponding edges. For example, if $e_1$ and $e_3$ are swapped on $T$, then $f_i\left(e_1\right)$ and $f_i\left(e_3\right)$ are swapped on $T|Y_i$.

2. Otherwise, the topologies of $T_{NNI}|Y_i$ and $T|Y_i$ are identical, $T_{NNI}$ and $T$ belong to a partial terrace. Thus, we keep the tree topology of $T|Y_i$ and have to update $f_i(e)$ according to Equation (4).

Under the EUL partition model when computing the log-likelihood of $T_{NNI}$ we only have to compute the log-likelihood of $T_{NNI}|Y_i$ when rule 1 applies. Otherwise the optimal log-likelihoods of $T|Y_i$ and $T_{NNI}|Y_i$ are equal. This advantage comes from the fact that under the EUL model the edge lengths of partition trees are optimized independently. Therefore, if the topologies of $T_{NNI}|Y_i$ and $T|Y_i$ are identical, there is no need to optimize edges. Thanks to this, the computing time for the EUL model greatly benefits from partial terraces. In fact, when $T_{NNI}$ and $T$ belong to one full terrace, that is, when for all partitions condition (C.2) is not satisfied, no likelihood recomputation is necessary.

More restrictive EL models require additional care. Since under EL models the edge lengths of the species tree and partition trees are linked, together with the topological changes of $T|Y_i$ we also have to account for changes in edge lengths. Using the map $f_i(.)$, the edge lengths of the partition tree $T|Y_i$ are computed as

where $r_i$ is the rate of partition $Y_i$. Therefore, each time map $f_i(.)$ changes, the edge lengths on $T|Y_i$ will also change. If the topology of $T|Y_i$ is not changed by the NNI (i.e., condition (C.2) is not satisfied), four cases are possible, which lead to varying computational time savings (see Appendix 2).

For EL models PTA helps to save computation time during edge length optimization. To speed up the optimization of each edge $e\in E$, in Equation (1) one only has to sum the log-likelihoods over those partitions $Y_i$ where $f_i(e)\ne ϵ$. This advantage is a result of using induced partition trees instead of complete partition trees (i.e., with a full set of species instead of $Y_i)$.

Since NNI changes one split of the species tree, namely, the split associated with the edge NNI is applied to, $F$ is recomputed in $O\left(k\right)$ time. Thus, the extra computations needed to maintain $F$ are negligible compared with the expensive likelihood computations.

Although we only reformulated the condition for NNI, we note that a similar reformulation applies to SPR and TBR (see Appendix 3). Thus, the PTA can be employed with all common topological rearrangements.

CPU Time Comparison

For each partition model and each alignment we computed: (i) the average CPU time for 10 runs for each program; and (ii) the speedup of ${\text{IQ-TREE}}_{\text{PTA}}$ compared with other implementations: the ratio between the average CPU time of each program and that of ${\text{IQ-TREE}}_{\text{PTA}}$.

Under the EUL model ${\text{IQ-TREE}}_{\text{PTA}}$ was the fastest program for 11 out of 12 alignments (Table 3a and Fig. 3a). Averaging over all alignments ${\text{IQ-TREE}}_{\text{PTA}}$ runs about three times faster than the standard IQ-TREE and 3.79 and 3.69 times faster than RAxML and ${\text{RAxML}}_{\text{optU}}$, respectively. For four alignments (DNA5, DNA7, DNA8, and DNA9), IQ-TREE was slower than RAxML and ${\text{RAxML}}_{\text{optU}}$ (Table 3a). However, ${\text{IQ-TREE}}_{\text{PTA}}$ improved the runtimes for all these alignments compared with IQ-TREE and ${\text{IQ-TREE}}_{\text{PTA}}$ was faster than RAxML except for DNA5.

Under the EL-equal partition model ${\text{IQ-TREE}}_{\text{PTA}}$ was on average 2.0 times faster than IQ-TREE (Table 3b and Fig. 3b). The smallest speedups correspond to alignments with the lowest percentages of missing data (DNA1, DNA5, AA10, AA11, AA12). ${\text{IQ-TREE}}_{\text{PTA}}$ was slower than RAxML and ${\text{RAxML}}_{\text{optU}}$ for DNA5 and DNA7, equally fast for DNA9 and faster for the remaining nine alignments.

Finally, since the EL-proportional model is not available in RAxML it was only examined with IQ-TREE and ${\text{IQ-TREE}}_{\text{PTA}}$. ${\text{IQ-TREE}}_{\text{PTA}}$ is on average 2.0 times faster than the standard IQ-TREE (Table 4 and Fig. 3c).

Figure 3 shows that ${\text{IQ-TREE}}_{\text{PTA}}$ is typically faster than the other three programs. We also computed the speedup between pairs of runs with the same starting trees (Online Appendix Fig. S1, available on Dryad at http://dx.doi.org/10.5061/dryad.v02t3). While ${\text{IQ-TREE}}_{\text{PTA}}$ found ML trees faster than IQ-TREE for all the cases, except runs for DNA5 under EL-proportional model, the running times of RAxML and ${\text{RAxML}}_{\text{optU}}$ were not significantly different.

Log-Likelihood Comparison

The results discussed in this section are based on the log-likelihoods reported by the respective programs. Note that after a major bug fix in RAxML included in version 8.1.24, the log-likelihoods reported by RAxML and IQ-TREE are directly comparable. That means that for the same tree topology with given branch lengths and same model parameters, the two programs return virtually identical log-likelihoods (Stamatakis A., personal communication).

Under the EUL model the log-likelihoods of the best-found trees obtained by different programs are quite similar with maximal difference of 50 (Fig. 4a) except for DNA1, DNA3, and DNA6, where the ${\text{IQ-TREE}}_{\text{PTA}}$ ML trees have log-likelihoods, which are 200, 190, and 100 units higher than RAxML trees, respectively. The log-likelihoods from 10 ${\text{IQ-TREE}}_{\text{PTA}}$ runs have consistently small variance, whereas the other programs sometimes show large variances (RAxML for DNA3 and DNA6, IQ-TREE for DNA9).

Under the EL-equal model the log-likelihoods of the best-found trees obtained by different programs are similar for most alignments (Fig. 4b) except for DNA3. The DNA3 alignment showed the largest variance of log-likelihoods over 10 runs: $\pm 150$ units for RAxML and $\pm 100$ units for ${\text{IQ-TREE}}_{\text{PTA}}$.

Under the EL-proportional model the average log-likelihoods for ${\text{IQ-TREE}}_{\text{PTA}}$ and IQ-TREE runs agree on nine alignments (Fig. 4c) while not for DNA3, DNA6, and DNA9. For DNA3 and DNA6 ${\text{IQ-TREE}}_{\text{PTA}}$ found trees that have on average 100 units higher log-likelihoods than IQ-TREE, whereas for DNA9 the ${\text{IQ-TREE}}_{\text{PTA}}$ trees showed 250 units smaller log-likelihoods. Finally, for DNA3 IQ-TREE showed large variance of $\pm 100$ units in log-likelihoods obtained from 10 runs.

The tree searches with RAxML and ${\text{RAxML}}_{\text{optU}}$ given the same starting tree resulted in identical ML trees. While the tree searches with IQ-TREE and ${\text{IQ-TREE}}_{\text{PTA}}$ with the same set of starting trees resulted in topologically different trees, but with sometimes identical log-likelihoods (Online Appendix Figs. S2–S4, available on Dryad). The search algorithm in both versions is the same, but re-optimization of some edge lengths can take place in IQ-TREE when an induced partition topology is not modified by the topological rearrangement, whereas ${\text{IQ-TREE}}_{\text{PTA}}$ leaves some or all edge lengths unchanged. This can create small numerical differences in likelihoods and therefore at some point different NNIs can be preferred by IQ-TREE and ${\text{IQ-TREE}}_{\text{PTA}}$, which results in different search paths. As a consequence this might lead to topologically different trees obtained by the two versions. These findings indicate that accounting for partial and full terraces influences the tree search and should not be ignored. The influence of partial/full terraces deserves further investigation and the development of new terrace-aware tree searches.