The large-sample asymptotic behaviour of quartet-based summary methods for species tree inference.

methods seek to infer a species tree from a set of gene trees. A desirable property of such methods is that of statistical consistency; that is, the probability of inferring the wrong species tree (the error probability) tends to 0 as the number of input gene trees becomes large. A popular paradigm is to infer a species tree that agrees with the maximum number of quartets from the input set of gene trees; this has been proved to be statistically consistent under several models of gene evolution. In this paper, we study the asymptotic behaviour of the error probability of such methods in this limit, and show that it decays exponentially. For a 4-taxon species tree, we derive a closed form for the asymptotic behaviour in terms of the probability that the gene evolution process produces the correct topology. We also derive bounds for the sample complexity (the number of gene trees required to infer the true species tree with a given probability), which outperform existing bounds. We then extend our results to bounds for the asymptotic behaviour of the error probability for any species tree, and compare these to the true error probability for some model species trees using simulations.

Introduction

A central problem of phylogenetics is the inference of the evolutionary history of a group of species in the form of a species tree, where leaves represent extant species and internal vertices represent speciation events. Likewise, the gene families contained within the genomes of the species have their own evolutionary histories, which can also be represented by trees that can be considered as evolving ‘within’ the species tree. This implies that the gene trees are heavily dependent on the species trees.

Historically, species trees have been inferred using so-called ‘concatenation’ methods, where the sequences of individual genes are concatenated together to form a set of ‘species sequences’, which are then used as input to a tree-building method of choice, often maximum-likelihood based (e.g., RAxML, Stamatakis 2014). However, this carries the implicit assumption that the individual gene trees are identical to the species tree. This is not necessarily the case, due to a number of processes which cause discordance between the gene and species trees.

One such process is incomplete lineage sorting (Maddison 1997), which occurs when gene lineages do not coalesce immediately (going backwards in time) within a single species. This can cause branch length differences between the gene and species trees, and if alleles are maintained over a sequence of rapid speciations, even topological differences. Thus in recent years there has been a rapid growth of ‘summary’ methods for species tree inference, where gene families are identified and a gene tree built from the sequences of each individual family. The set of gene trees are then summarised in some way to produce a species tree. This paradigm acknowledges and incorporates the fact that gene trees may differ from each other and the species tree.

A popular summary method is ASTRAL (Mirarab et al. 2014; Mirarab and Warnow 2015), which works by solving the ‘maximum quartet support species tree’ (MQSST) problem. Each gene tree is decomposed into a number of gene quartets, with topologies determined by the gene tree. ASTRAL then uses a dynamic programming algorithm to find the species tree which maximises the number of gene quartets that agree with it. As this problem is NP-hard (Lafond and Scornavacca 2019), ASTRAL has two running modes: one which solves this problem exactly, and a faster heuristic which restricts the search space but does not guarantee optimality. We consider only the first mode in this paper.

The dominant statistical model for incomplete lineage sorting is the multispecies coalescent (MSC, Rannala and Yang 2003), which considers each species branch as a separate population in which Kingman’s coalescent (Kingman 1982) is run. The lineages at the top of each branch (going backwards in time) then become input to the coalescent at the bottom of its parent branch. This model provides exact probabilities of each gene tree (and gene tree topology) given a specified species tree. It has been shown for both unpartitioned and partitioned methods (Roch and Steel 2015; Roch et al. 2019) that under this model, there are species trees where concatenation methods are not statistically consistent: they do not reconstruct the true species tree with probability approaching 1 as the number of input gene trees becomes arbitrarily large. Conversely, ASTRAL is statistically consistent (for all species trees) under the multispecies coalescent (Mirarab et al. 2014), when input gene trees are considered to be correctly reconstructed from sequences.

In addition, it is well known that incomplete lineage sorting is not the only cause of discordance between gene and species trees. Discordance can also arise due to other genetic evolutionary processes such as gene duplication, gene loss, and lateral genetic transfer. Models of gene evolution incorporating duplication and loss only (DL) are well-established, dating back to Goodman (Goodman et al. 1979), as are models with duplication, loss, and transfer (DTL; see Doyon et al. 2011, for a review). These generally consider duplications and losses (and sometimes transfers) to occur as independent linear birth-death processes within each species branch. More recently, some unified models (DLCoal, Rasmussen and Kellis 2012, and MLMSC, Li et al. 2021) have been developed that account for both DTL and the population-level processes that cause incomplete lineage sorting.

For each of these models, it is of interest to determine if ASTRAL (and other summary methods) is statistically consistent. This becomes progressively more difficult as the models increase in complexity, but it has been proven that this is indeed the case for the DL model (Legried et al. 2021) and the DLCoal model (Markin and Eulenstein 2020). A numerical study (Yan et al. 2021) also looked at the practical accuracy of ASTRAL under the DLCoal model. The consistency of ASTRAL under the MLMSC (or, indeed, any models with transfers) remains unknown.

We approach this problem in a general way by not restricting our consideration to the MSC or any other particular model of gene evolution. Instead, we consider the asymptotic behaviour of the probability of inferring the wrong species tree (the error probability) as the number of input gene trees, N, becomes large. Rather surprisingly, while there has been considerable interest on whether the error probability goes to 0 or not (i.e., statistical consistency), no-one (to the best of our knowledge) has studied the asymptotic behaviour of the error probability when it does go to 0. Our results below show that the decay is exponential in nature, which is unsurprising but does not seem to have been theoretically proven before.

A very closely related work is that by Shekhar et al. (2017), who studied the sample complexity (the number of gene trees required to infer the true species tree with a given probability) of ASTRAL. They derived bounds for this quantity under the multispecies coalescent, and recently Hill et al. (2020) derived corresponding results under the DLCoal model. The focus of these papers is on the behaviour of the sample complexity as the length of the shortest branch of the species tree goes to 0, for a fixed (but small) error probability. Conversely, our results study the asymptotic behaviour of the error probability as the number of input gene trees becomes large, for a fixed species tree. However, their bounds are comparable in some situations to our results.

In this paper, we study the asymptotic behaviour of the error probability as . We first study the case of a 4-taxon species tree, and derive a closed form for the asymptotic behaviour in terms of the probability that the gene evolution process produces the correct topology. We also provide bounds for the error probability that hold for all N. This enables us to numerically derive corresponding bounds on the sample complexity for ASTRAL under the MSC, which outperform the bounds from Shekhar et al. (2017) in practice. We then extend our results to species trees with arbitrary numbers of taxa, resulting in bounds (both asymptotic and for all N) for the error probability for any species tree. Using simulations, we compare our bounds with the true error probability for some sample topologies.

4-taxon species tree

We first study the case where the species tree S has 4 taxa. We suppose that N unrooted gene trees are simulated from the species tree under an unspecified gene evolutionary process that has probability to produce the correct topology S, and equal probability to produce either of the other two possible topologies. (Although it is possible that the incorrect topologies may be produced with unequal probability, many models, such as the MSC, will produce them with equal probability.) These trees are then used as input to a quartet-based summary method to infer a species tree. Although we refer to this method as ASTRAL throughout this paper, we only require that it solves the MQSST problem. (In the case where the input gene trees have an equal maximum number of two or more different topologies, we assume that the method randomly chooses one to infer as the species tree.) Since , we know that this method is statistically consistent.

We begin by defining a function that we will use repeatedly throughout this paper.

Definition 1

Let f(p) be the functionwhere

We illustrate the behaviour of this function in Fig. 1. Important aspects to note are that it is defined (for our purposes) on the domain , is always negative on this domain, and monotonically decreases to an asymptote at .

Fig. 1

The function f(p)

We now derive the central result of this paper: we show that the asymptotic behaviour of the error probability goes to 0 exponentially with N, and derive a closed form for the exponential constant that only depends on the probability p of producing the correct topology.

Theorem 1

Given a 4-taxon species tree, let be the probability that the gene evolution process produces the correct topology. As the number of gene trees , the probability of inferring the wrong species tree with ASTRAL behaves as

Proof

Let be the number of gene trees produced with the three topologies, with being the correct one. (Obviously, .) Then has a multinomial distribution with parameters N and .

The correct species tree will be inferred if . Consider now the case for some . In this case, the wrong species tree will be inferred with probability . The probability of this event is:(For brevity, we leave the out in the following, as it is implied by .)

Using Stirling’s formula ,We now choose c to maximise the (log-)probability of this event:Setting this to 0 yieldsNote that there is no guarantee that cN is an integer for this value of c; however, as , it will become arbitrarily close to such a value.

We now claim that the probability of inferring the wrong species tree has the same asymptotic behaviour as . Since this event will produce a wrong species tree with probability , we haveTo show the opposite bound, consider that there are possible combinations of which add to N (we can choose such that , and then is determined). Of these, produce the wrong species tree by symmetry. For each of these combinations , we assume without loss of generality that . Then one of the following cases occurs:We can recursively apply the latter two relations until , then use the first case. Thus each of these combinations has a probability which is bounded above by .

, so by construction of c;

and , so

, so

Hence (with a slight abuse of notation, as the inequality only holds in the asymptotic limit)Since the two last terms are o(N), the claim is proved, as is the result.

Note. An initial attempt at deriving this result involved using a normal approximation to the multinomial; rather surprisingly (to us), this failed because the approximation of multinomial tail probabilities became inaccurate faster than the normal approximation became accurate.

This theorem leads to an immediate result on the Robinson-Foulds (RF) accuracy of the inferred species tree.

Corollary 2

Under the conditions of Theorem 1, as the number of gene trees , the average RF distance between S and the tree inferred by ASTRAL behaves as

This follows immediately from Theorem 1 and the relation .

Theorem 1 can be extended to derive lower and upper bounds on the error probability, as we show in the following theorem. These bounds apply for all N, not just in the limiting case.

Theorem 3

Under the conditions of Theorem 1, the probability of inferring the wrong species tree with ASTRAL is bounded below bywhere is the nearest number tosuch that is an integer. Likewise, it is bounded above by

The proof uses a bounded version of Stirling’s formula and is given in Appendix A.

Finally, we note that, depending on the gene evolution process, it is possible that the two incorrect topologies may not be generated with equal probability. This cannot happen for the MSC, but in models which include duplication it may occur for species trees that are not ultrametric, or if duplication rates are allowed to vary between branches. This can, and more easily, also happen in models including lateral gene transfers or introgression. In these cases, Theorem 1 can be generalised to the following result.

Theorem 4

Given a 4-taxon species tree, let the probabilities of producing the three quartet topologies be , where is the probability of producing the correct topology and . As the number of gene trees , the probability of inferring the wrong species tree with ASTRAL behaves aswhere

We sketch a proof in Appendix B, and also state a generalisation of Theorem 5 there. Note that, in models including lateral gene transfers or introgression, we can have extreme cases where the topology matching the species tree is not the most common one. In these cases, even our generalisations do not hold.

General species tree

We now consider the case of a species tree S with an arbitrary number n of taxa. As before, we start by assuming that for any quartet, the gene evolution process produces the two incorrect topologies with equal probability.

To derive a lower bound on the asymptotic behaviour of the error probability, we require an assumption on the gene evolution process. Consider an internal branch x of the species tree which divides the extant species into four clades , , , and . Now consider two species quartets and , such that one species is drawn from each clade for both quartets. We say that these two species quartets are positively correlated if the events of producing each of the three topologies are positively correlated between the two quartets; that is,and likewise for the other two topologies. We then say that the gene evolution process is positively correlated if this property holds for all such pairs of quartets for all internal species branches. It is easy to see that the MSC is positively correlated, as there is a positive probability of the gene lineages from and coalescing before they reach the root, and likewise for the remaining three pairs of species.

Theorem 5

Given a species tree S, assume the gene evolution process is positively correlated. Let be the minimum probability among all quartets in S of producing the correct topology. For any internal species branch , let , , , and be the four clades that it divides the leaves of S into, and letThen the probability of inferring the wrong species tree with ASTRAL has the limiting behaviour (as )where the (logarithm of the) growth constant is bounded by

To derive the upper bound, a sufficient (but not necessary) condition to infer the correct species tree is that each of the species quartets have a plurality of gene quartets with the correct topology. (We use the word ‘plurality’ to mean that more gene quartets have the correct topology than any other, not that more than half the gene quartets have the correct topology.) Let be the event that species quartet i has a majority of gene quartets with the correct topology. ThenClearly, the probability of is maximised for the species quartet that has the minimum probability of producing the correct topology (the smaller the probability of producing the correct topology, the smaller the probability that the majority of quartets have the correct topology). The upper bound then follows from Theorem 1.

To derive the lower bound, consider an internal species branch x, and consider the species quartets that contain one species from each of the clades . There are such species quartets, which generate gene quartets. A necessary (but not sufficient) condition to infer the correct species tree is that the majority of these gene quartets have the correct topology. To see this, suppose that the correct species tree is inferred, but the majority of these quartets have the wrong topology. In this case, performing a nearest-neighbour interchange move on x to this wrong topology will result in a species tree that agrees with more of these quartets, and not change the agreement with any other quartets. This improves the quartet score, a contradiction.

Let be the event that the majority of these gene quartets have the same wrong topology.

Now although gene quartets from different gene families are independent, this is not the case for gene quartets from the same gene family. Since the gene evolution process is positively correlated, the probability of is bounded below by the probability that the majority of the gene quartets have the same wrong topology if they were all generated independently. Thus we haveWe now apply the same argument as in Theorem 1, with the difference that instead of N gene quartets, we have quartets. Since we can choose the x which maximises the resulting bound, this produces the lower bound.

Finally, because both the lower and upper bounds are O(N), we can say that the error probability also has this limiting behaviour as .

Note in particular that the upper bound does not depend on n, the size of the species tree, whereas the lower bound does.

While the derivation of the upper bound is very similar to that in Shekhar et al. (2017), the lower bound is qualitatively different. In that paper, they show the existence of a species tree for which the error probability is bounded below by a specified value. Our result is more general, in that it bounds the error for any species tree. (On the other hand, this means that their lower bound can be tighter than ours.)

Note also that under the MSC, the probability of producing the correct topology for a quartet depends only on the length of the internal branch; thus all species quartets (a, b, c, d) in (2) have exactly the same chance of producing the correct topology for a branch x, and there is no maximisation required to calculate .

It is worth considering what sort of branch x will tend to maximise the lower bound. Since f(p) is a negative and decreasing function, it is maximised for small p. However, we also wish to simultaneously minimise ; this generally occurs when x is a branch near the leaves, so that three of the four clades are small. The smallest possible value for this quantity is either (if x is the internal branch of a subtree ( (a, b), c)) or (if no such subtree exists and x is one of the internal branches of a subtree ( (a, b), (c, d) )). Note that the first case in particular corresponds to the construction used to prove Claim 2.1 in Shekhar et al. (2017). Therefore, the branch that maximises the lower bound will usually be a short branch close to the tips of the species tree; this is the kind of branch that is most difficult to infer correctly.

Using similar arguments to Theorem 3, we can bound the error probability for general trees for all N; we state this result without proof in Appendix D.

Sample complexity under the multispecies coalescent

The results in the previous sections apply under any model of gene evolution, requiring only the probability p of producing the correct topology. The most-studied model for this purpose is the multispecies coalescent, under which p depends only on the length of the internal branch l (in coalescent units), given by the relation . (We change the notation here from Shekhar et al. 2017, where the length of the internal branch is denoted f, but we have already used this notation for the asymptotic growth constant function.) Substituting this relation into Theorem 1 gives the following result for the sample complexity of ASTRAL under the MSC on a 4-taxon species tree.

Theorem 6

Under the multispecies coalescent, the minimum number of genes required to correctly reconstruct a 4-taxon species tree with ASTRAL with probability at least grows asymptotically (as ) aswhere g(l) is the functionAs ,

The first expression follows immediately from setting in Theorem 1.

To derive the asymptotic behaviour, we expand (1) as a Maclaurin series for c in l:Substituting this and into the expression for f(p) and expanding as a Maclaurin series in l gives the asymptotic behaviour for g(l). A more detailed derivation is given in Appendix C.

Corresponding bounds for the sample complexity are given in Shekhar et al. (2017). In the case , their bounds are (in our notation)where is the normal cdf. We see that the asymptotic behaviour of our expression agrees with the upper bound, as they are both . The lower bound is not directly comparable, as it only applies for values of l below a threshold depending on . Thus it is only suited to analysis of the case for fixed , not for fixed l.

We can derive bounds on the error probability for all N under the MSC from Theorem 3, by substituting the relation for p. In Fig. 2, we compare the asymptotic behaviour for the error probability against our bounds and the upper bound of Shekhar et al. (2017). We can see that our upper bound outperforms that of Shekhar et al. (2017) in the limit , and in general for large l; however theirs works better for small N and l. We can also see that the asymptotic behaviour lies between our bounds, which is expected in the asymptotic limit but appears to hold for all N anyway.

Fig. 2

The logarithm of the error probability (for a 4-taxon species tree) vs N, for fixed p (fixed l under the MSC). We show the asymptotic behaviour (black), our bounds (blue), and the upper bound of Shekhar et al. (2017) (red) (color figure online)

We can also derive bounds for the sample complexity from Theorem 3 by solving for N numerically (an analytical inversion is likely impossible). We show the results of this in Fig. 3. We can see that our bounds, particularly the upper bound, outperform the bounds of Shekhar et al. (2017). This is because we require very large N to reach the desired error probability, so the better asymptotic properties of our bounds result in tighter bounds than that of Shekhar et al. (2017), which outperform ours for low N but have worse asymptotic properties.

Fig. 3

The sample complexity of ASTRAL (for a 4-taxon species tree) vs , for fixed . We show the asymptotic behaviour (black), our bounds (blue), and the bounds of Shekhar et al. (2017) (red) (color figure online)

Finally, we state bounds for the sample complexity under the MSC for the case of a species tree with an arbitrary number n of taxa.

Theorem 7

Under the multispecies coalescent, the minimum number of genes required to correctly reconstruct a species tree with ASTRAL with probability at least grows asymptotically (as ) aswherewhere is the length of the shortest internal species branch and is the branch that maximises the lower bound in (3).

Simulations

We follow the simulations of Shekhar et al. (2017) and study three model species trees with species, shown in Fig. 4. The caterpillar tree is maximally unbalanced, while the balanced tree is maximally balanced; both have all internal branches with length l. The ‘double-quartet’ tree has the same topology as the balanced tree, but with a very long internal branch (separating the two quartets) and all other branches of length l.

Fig. 4

The three model trees. The branch that provides the lower bound is highlighted in red (color figure online)

The branch that provides the best lower bound in Theorem 5 for each tree is highlighted in Fig. 4. Because all internal branches (save the long branch in the double-quartet tree) have the same length (and therefore the same probability of producing the correct topology under the MSC), this branch is the one that minimises . The value of this quantity is 5, 8, and 8 for the caterpillar, balanced, and double-quartet trees respectively, which also gives the factors by which the lower bound is (negatively) greater than the upper bound. For these trees, the upper bounds are identical.

We simulate N independent gene trees under the MSC model for a range of N from 10 to 60000. The gene trees are then input into ASTRAL to infer a species tree, which is then compared with the true species tree. This is done for a range of l from 0.005 to 0.1 (with evenly spaced as in Shekhar et al. 2017), with varying numbers of replicates (always at least 2000) for each value of N and l. We then fit a linear regression to against N to calculate the growth constant .

In Fig. 5, we show the growth constants for the three trees against , together with our asymptotic bounds from Theorem 5 and the upper bound of Shekhar et al. (2017). (As mentioned above, the lower bound of Shekhar et al. 2017, cannot be compared here.) We can clearly see that our formulas do indeed bound the asymptotic growth constants, and moreover that we provide a tighter upper bound on the asymptotic behaviour than Shekhar et al. (2017).

Fig. 5

Growth constants for for the three model trees, with our asymptotic bounds (blue) and the upper bound of Shekhar et al. (2017) (red) (color figure online)

In general, the true growth rate is closer to our upper bound than our lower bound. This is unsurprising considering that the lower bound is derived by considering all quartets corresponding to the same internal branch as independent, a large relaxation for the highly positively-correlated MSC. We also observe the pattern observed in Shekhar et al. (2017) that ASTRAL performs worse (the growth rate of the log-error probability is smaller negatively) on the double-quartet tree than the other two trees; see Shekhar et al. (2017) for a discussion of why this is so.

Discussion

In this paper, we have derived the asymptotic behaviour of the error probability (the probability of inferring the wrong species tree) of ASTRAL in the large-sample limit (as the number of input gene trees becomes large). In particular, we show that the error probability goes to 0 exponentially with respect to N, and derive a closed form for the growth constant of this exponential behaviour exactly for the simplest case of a 4-taxon species tree. We calculate bounds on the growth constant for species trees with arbitrary number of taxa, and extend our results to rigorous bounds for any N for both the 4-taxon and general species tree. Our results improve previous bounds on the sample complexity from Shekhar et al. (2017).

It is true that in practice each species only has a finite number of genes (where a gene is simply a non-recombining locus in a genome); thus N has a natural boundary which it cannot go beyond. On the other hand, the increasing amount of whole-genome data available does allow a greater and greater number of genes (and therefore gene trees) to be used in species tree inference, and so an exploration of the behaviour of these methods in the large-N limit is of practical importance. In particular, our results are good enough to be useful in the question of determining the sample complexity, which is of great practical interest as it allows one to estimate a priori how many gene trees need to be extracted from the genome for use in species tree inference.

Our results are quite general, in the sense that although they are expressed in relation to ASTRAL, they apply for any species tree inference method which solves the maximum quartet support species tree problem. Nor are they limited to the multispecies coalescent model; although this is the dominant model at the moment, this may change with the recent rise of unified models that also include DTL processes, such as DLCoal or MLMSC. Our results (possibly apart from the lower bound for n species, which additionally requires positive correlation) apply equally well to these models, and to future models that may additionally include other processes. We require only the ability to calculate the probability of the gene evolution process generating a correct quartet topology. While this is not always an easy task theoretically under the more complex models, it can often be estimated very accurately via simulations on 4-taxon species trees (vastly more efficient than conducting simulations over many species).

One restriction on the gene evolution process that we require for our lower bound is that it must be positively correlated. This seems a natural requirement, and it is indeed easy to show that it holds for the MSC (as discussed above). However, it is not trivial to extend this argument to more complicated models such as DLCoal or MLMSC, and it remains an open question whether these models are indeed positively correlated or not.

Finally, we note that the important practical consideration of whether the gene trees are reconstructed with or without error does not affect our theoretical results in Sects. 2 and 3. We only require the probability p of producing the correct quartet topology (or in the more general case, the probabilities of the three topologies). If there is gene tree error, we can simply consider p as the probability of reconstructing, with error, the correct topology (likewise in the more general case). For our results to hold, it only matters that , so that we are in the statistically consistent regime. As discussed in Shekhar et al. (2017), this will hold under a simple model where the gene tree error proportion is less than or equal to , and errors are unbiased and independent between quartets.

For Sects. 4 and 5, the presence of gene tree error does affect the formula , and so the asymptotic behaviour for the sample complexity becomes a lower bound in practice.