An algorithm for computing the gene tree probability under the multispecies coalescent and its application in the inference of population tree.

MOTIVATION: Gene tree represents the evolutionary history of gene lineages that originate from multiple related populations. Under the multispecies coalescent model, lineages may coalesce outside the species (population) boundary. Given a species tree (with branch lengths), the gene tree probability is the probability of observing a specific gene tree topology under the multispecies coalescent model. There are two existing algorithms for computing the exact gene tree probability. The first algorithm is due to Degnan and Salter, where they enumerate all the so-called coalescent histories for the given species tree and the gene tree topology. Their algorithm runs in exponential time in the number of gene lineages in general. The second algorithm is the STELLS algorithm (2012), which is usually faster but also runs in exponential time in almost all the cases.
RESULTS: In this article, we present a new algorithm, called CompactCH, for computing the exact gene tree probability. This new algorithm is based on the notion of compact coalescent histories: multiple coalescent histories are represented by a single compact coalescent history. The key advantage of our new algorithm is that it runs in polynomial time in the number of gene lineages if the number of populations is fixed to be a constant. The new algorithm is more efficient than the STELLS algorithm both in theory and in practice when the number of populations is small and there are multiple gene lineages from each population. As an application, we show that CompactCH can be applied in the inference of population tree (i.e. the population divergence history) from population haplotypes. Simulation results show that the CompactCH algorithm enables efficient and accurate inference of population trees with much more haplotypes than a previous approach. AVAILABILITY: The CompactCH algorithm is implemented in the STELLS software package, which is available for download at http://www.engr.uconn.edu/ywu/STELLS.html CONTACT: ywu@engr.uconn.edu SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 Introduction

Consider n gene lineages sampled from one population. When we trace these lineages backward in time, sooner or later two of these lineages will find a common ancestor. When this occurs, we say these two lineages coalesce or the coalescent of these two lineages occurs. Gene lineages coalesce in a stochastic way, which is influenced by multiple population genetic parameters, such as population sizes. Therefore, coalescents may potentially reveal important aspects of population evolution. In 1982, Kingman (1982) introduced the coalescent theory, which provides the analytical foundation for the study of coalescents. Since then, coalescent theory has quickly become a very active research subject in population genetics. There are numerous theoretical results in coalescent theory, along with many coalescent-based software tools. See Wakeley (2008) and Hein ) for an overview of the growing field of coalescent theory.

While mathematically appealing, coalescent is known to be challenging computationally. One of the most important problems on coalescents is computing the coalescent likelihood. That is, we want to compute the probability of observing some population variation data under a coalescent model. However, coalescent likelihood is well known to be difficult under most formulations [see e.g. Hein ; Wakeley, 2008]. At present, there are few polynomial time algorithms for computing coalescent likelihood except under very restricted conditions. Almost all existing approaches for coalescent likelihood computation are based on statistical techniques such as Markov chain Monte Carlo (MCMC) or importance sampling (Wakeley, 2008). These statistical approaches can be very useful in practice. However, from the computational point of view, these approaches cannot compute the exact coalescent likelihood. Moreover, approaches using MCMC tend to be very slow. Since coalescent likelihood is usually used in the inference of population evolution, accurate and efficient coalescent likelihood computation is highly desirable.

Since the work by Kingman (1982), the basic coalescent has been extended to address more aspects of population evolution. The multispecies coalescent is one such extension, and is the focus of this article. Multispecies coalescent concerns the coalescents among gene lineages from multiple populations where coalescents may cross the population boundary. Multispecies coalescent is fundamental to population evolution (Degnan and Salter, 2005; Rosenberg, 2002). Multispecies coalescent has also found applications in other domains, such as species tree inference [see, e.g. Heled and Drummond, 2010; Mirarab ; Wu, 2012].

A key computational problem on multispecies coalescent is the computation of the gene tree probability (Degnan and Salter, 2005; Wu, 2012, 2015). Here, we are given a gene tree topology (without branch lengths) and a species tree (with branch lengths in the standard coalescent units). Note that the gene tree may not be bifurcating and there can be more than one gene lineages for a population (or taxon) in the gene tree. The gene tree probability is the probability of the gene tree topology as a result of the multispecies coalescent within the given species tree. There are only two existing algorithms for the exact gene tree probability computation. The first such algorithm is due to Degnan and Salter (2005). However, their algorithm has the exponential running time and there are no known cases of species and gene trees for which their algorithm runs in polynomial time. A much faster algorithm is given in Wu (2012), which is implemented in the program STELLS. However, the STELLS algorithm computes the gene tree probability in polynomial time only for some very special types of gene trees and species trees. We are not aware of any statistical approaches for computing the gene tree probability. The closest related approach is the MCMC approach in Heled and Drummond (2010), which is also based on multispecies coalescent. However, the likelihood computed in Heled and Drummond (2010) is for DNA sequences, not for gene trees.

In this article, we present a new algorithm, called CompactCH, for computing the exact gene tree probability. The main advantage of CompactCH is that it runs in polynomial time for any gene tree topologies and species trees when the number of taxa (populations) is fixed to be a constant. That is, the CompactCH algorithm is efficient for gene tree topologies with large number of gene lineages if the number of populations is small (i.e. there can be many gene lineages from each population). Note that although the number of gene lineages in many phylogenetic studies tends to be relatively small, the number of gene lineages in large-scale genetic studies can be large. Therefore, CompactCH may become more useful when more gene lineages are sampled from a population which is the case in genetic projects such as the 1000 Genomes Project (2015). To the best of our knowledge, CompactCH is the first algorithm for computing the gene tree probability in polynomial time that allows arbitrary gene trees and species trees when the number of populations is fixed to a constant. Note that if gene tree topologies are considered to be observed, then the gene tree probability is the likelihood of the observed gene tree topologies given the species (population) tree. In this regard, the coalescent likelihood of gene tree topologies can be computed in polynomial time under the assumption of small number of populations. We show through theoretical analysis and empirical simulations that CompactCH outperforms the existing STELLS algorithm when the number of populations is small and the gene tree contains multiple gene lineages (alleles) per taxon. We also present a new method for the inference of population divergence history (called the population tree) from large population samples, which is based on the CompactCH algorithm.

2 Background

2.1 Multispecies coalescent and incomplete lineage sorting

In this article, we use population and species interchangeably. Multispecies coalescent process (or just multispecies coalescent) concerns the coalescents of gene lineages that may occur outside the species boundary. See Figure 1 for an illustration. Multispecies coalescent lies at the intersection of population genetics and phylogenetics, and has been widely used in phylogenetic and population inference [see e.g. Heled and Drummond, 2010; Wu, 2012]. We use a rooted bifurcating tree T (with branch lengths) to represent the population divergence history. We use a rooted tree topology T to represent the genealogical history of gene lineages at a single locus. It is important to note that T does not have branch lengths. Also T may be multifurcating if the gene genealogy cannot be fully determined from the given genetic data. A common observation is that T and T do not always have the same topology. This phenomenon is called incomplete lineage sorting. We define the gene tree probability to be the probability of observing T for a given T under the multispecies coalescent model. Computing the probability of gene genealogies under multispecies coalescent is central both to phylogenetic inference of a species (or population) tree and to population genetic inference of gene flow between populations [e.g. Degnan and Salter, 2005; Heled and Drummond, 2010; Hudson, 1983; Rosenberg, 2002; Wu, 2012, 2015]. When T and T are small, simple closed-form formulas are known to calculate the gene tree probability (Hudson, 1983; Rosenberg, 2002; Takahata and Nei, 1985). When T and T become larger, computation by hand is no longer feasible and an algorithm is needed.

Fig. 1.

A gene tree (in thin lines) in the species tree (in thick lines) shown in (a). Gene lineages a1, a2 and a3 originate from species A, b1 and b2 from B and c1 from C. The species tree T is shown separately in (b), so is the gene tree T in (c). Internal nodes of both T and T are labeled. Coalescents: internal nodes of T. Branches of trees are represented by their lower nodes

2.2 Computing gene tree probability using coalescent history

Degnan and Salter (2005) gave the first algorithm for computing the gene tree probability of a gene tree topology T and a species tree T. Their algorithm is based on the notion of coalescent history. We give parts of technical details of their algorithm because our new algorithm is an improvement over their algorithm. In the following, we use the lower (i.e. closer to the leaves) node of an edge to refer to the branch in a tree T. For example, when we say branch g in Figure 1c, we mean the branch with lower node g in T. Each internal node v in T corresponds to a coalescent event between two lineages. For simplicity, we use v to also refer to the coalescent at v. For example, in Figure 1, the node g corresponds to a coalescent between lineages a1 and a2. A coalescent history determines, for each coalescent v in T, on which species tree branch v occurs. In Figure 1a, for example, the embedded T within T corresponds to a coalescent history. In this history, coalescent g occurs on species tree branch A, p on D and h, q and r on E (the root branch). Note that the precise positions of these coalescents within T are not important since T does not have branch lengths and thus branches of T may be stretched or shrunk as long as the coalescents occur within the same species tree branches.

Degnan and Salter’s algorithm enumerates all possible coalescent histories for the given T and T. We denote the set of all coalescent histories for T and T as . We let m be the number of branches in T and use an integer to refer to a branch in T. We let λ be the branch length of the branch b in T. Then.

Here, we call the number of gene lineages at the bottom (respectively top) of a species tree branch b the lower (respectively upper) lineage count of b and denote as u (respectively v). Recall that a coalescent history specifies along which species tree branch each coalescent occurs. Thus, when is given, both the upper and lower lineage counts of each species tree branch are known. So we use and to refer u and v specified by . For example, consider the branch A and let be the coalescent history shown in Figure 1a. Then and . is the probability of u lineages coalesce into v lineages within time T (where T is the in the standard coalescent units). Equation (2) is a classic result (Takahata and Nei, 1985; Tavarè, 1984; Watterson, 1984) [also see Wakeley (2008) and Rosenberg (2002)] in coalescent theory which gives a closed form formula for calculating .

The reason why there are terms in Equation (1) is that does not impose the order of coalescents as specified in T. For example, consider the coalescent history shown in Figure 1a. We have and . Then, only requires four lineages coalesce to a single lineage within time λ, but does not impose the condition that g and b1 must coalesce and also p and c1 must coalesce. That is, allows, for example, the coalescent of g and c1, which violates the topology of T. Degnan and Salter define as the number of all possible ways of coalescing lineages into lineages. Since any pairs of lineage can coalesce, we have (Degnan and Salter, 2005):

The term is a little more complex. It is equal to the number of ways of coalescents on branch b as specified by where these coalescents match the topology of T. Note that there are coalescents on branch b. We denote these coalescents as . For convenience, we use to denote the number of coalescents along a species tree branch b. There are ways of ordering these coalescents, but only a subset of these ways match T. For each coalescent , let be the number of coalescents in that are within the subtree rooted at node u in T. Here, u is included in the subtree. For example, we again consider the coalescent history in Figure 1a. Then and . Note only coalescents within branch E are considered here and thus the coalescents p and g are not counted for and . Note that among all possible permutations of coalescents (with u being one of the coalescents), we want those with u be placed after the coalescents. There is precisely one out of permutations with this property. Then, we have (Degnan and Salter, 2005):

As an example, for the coalescent history in Figure 1a,

The exact gene tree probability can be computed if all coalescent histories are enumerated for T and T. However, as shown in Degnan and Salter (2005), the number of coalescent histories grows rapidly when the sizes of T and T grow. In fact, when T and T have matching topology [assumed in most analytical results on coalescent history; see e.g. (Rosenberg, 2013)], there exists no known cases where lead to polynomial number of coalescent histories.

Recently, an algorithm called STELLS was developed in Wu (2012). The STELLS algorithm is based on dynamic programming on a data structure called ancestral configuration (or AC). The STELLS algorithm is much faster than the algorithm by Degnan and Salter. Moreover, it is known that the number of ACs is bounded by a polynomial when T and T have matching topology and the topology is the maximal-asymmetric (i.e. caterpillar) tree. So in this special case, STELLS computes the gene tree probability in polynomial time. However, except the simple caterpillar trees (or trees very similar to caterpillar), there are no known cases where STELLS can compute the gene tree probability in polynomial time. Usually the STELLS algorithm becomes slow when the sizes of T and T increase.

3 The CompactCH algorithm for computing gene tree probability

We now present a new algorithm, called CompactCH, for computing the exact gene tree probability. CompactCH builds on Degnan and Salter’s algorithm, and is very different from the STELLS algorithm. The key idea of the CompactCH algorithm is combining multiple coalescent histories into a compact coalescent history and thus making the computation in Equation (1) more efficient. For the ease of exposition, we assume gene trees are bifurcating unless otherwise stated. We will later extend to the multifurcating gene tree case.

3.1 Gene tree probability computation using compact coalescent history

Since our algorithm involves a number of notations, we provide a list of notations in the Supplemental Materials. We first recall that Equation (1) is a summation of products over coalescent history . Each product consists two terms for each species tree branch b: the coalescent factor term and the coalescent probability term . The coalescent probability term only depends on the upper and lower lineage counts of the branch b. Suppose we group the coalescent histories with the same upper and lower lineage counts for every species tree branch into a compact coalescent history . Then all these histories will have the same coalescent probability terms. More precisely, a compact coalescent history is a list of upper lineage counts at all species tree branches (except the root branch). For example, consider the coalescent history shown in Figure 1a. We arrange the upper lineage counts for species tree branches in the order of A, B, C and D. Then the corresponding compact coalescent history is represented as: {2, 2, 1, 3}. That is, there are two, two, one and three gene lineages at the top of A, B, C and D, respectively. We denote the set of all possible compact coalescent histories for the given T and T as . A compact coalescent history can be viewed as the combination of one or multiple coalescent histories. As an example, consider T and T in Figure 1. There are total seven compact coalescent histories as shown in Table 1. The compact history {2,2,1,3} combines two coalescent histories: the history shown in shown in Figure 1a and the history with h as the only coalescent occurring within branch D and p occurring within the branch E.

Table 1.

The list of all compact coalescent histories for the trees T and T in Figure 1

	CCH	#histories
ch_1	{3,2,1,5}	1
ch_2	{3,2,1,4}	2
ch_3	{3,2,1,3}	2
ch_4	{3,2,1,2}	1
ch_5	{2,2,1,4}	1
ch_6	{2,2,1,3}	2
ch_7	{2,2,1,2}	1

CCH: compact coalescent history (the numbers are the numbers of upper lineage counts for species tree branches A, B, C and D). For each compact history, we give the number of coalescent histories that are merged into this compact history (denoted as #histories).

While a compact coalescent history only specifies the upper lineage counts, Lemma 3.1 shows that both upper and lower lineage counts as well as the number of coalescent events on all species tree branches (including the root branch) are all fully determined by a compact coalescent history.

The lower and upper lineage counts, and the numbers of coalescent events along all species tree branches are determined by a compact coalescent history.

Recall a compact coalescent history specifies the numbers of gene lineages on top of each species tree branch (except the root branch). Note that for the root r, the upper lineage count v is always one. This is because all gene lineages must coalesce into one lineage within the root branch.

Consider a species tree branch b. If b is a leaf branch (i.e. one of the nodes of b is a leaf in T), its lower lineage count is equal to the number of sampled linages from this leaf, which is given as part of the input. This is because there is no time for coalescent to occur at a leaf of T.

Let b be an internal branch. Let Desc(b) be the set of descendant branches of T (i.e. these branches are the outgoing edges from the lower node of b). For example, in Figure 1, and . Then, . This is because the lineages at the top of the descendant branches of b enter b from below. Thus, the set of lineages at the lower node of b is simply the merged set of the sets of lineages at the top of all b’s descendant branches. Note that this merging process is instantaneous and there is no time for coalescent to occur.

Note that the number of coalescents on b is simply equal to the difference between the lower and the upper lineage counts at b. Therefore, the number of coalescent events on b is fully determined by the compact history for all b.

From Lemma 1, we use and to denote the lower and upper lineage count at branch b, and as the number of coalescents on b, which are specified by the compact coalescent history .

We now consider Equation (1) again. Note that each coalescent history can be mapped to a compact coalescent history. For a compact coalescent history , let be the set of coalescent histories combined in . We now group the coalescent histories in together for each when applying Equation (1). Note that for all histories in , (respectively ) are the same. Thus they have identical coalescent probability terms terms in the summation over . Moreover, we have the same terms. This is because according to Equation (3), only depends on and , which are identical for histories in . This allows us to extract the common product of the terms out of the summation over . Therefore,

In Equation (5), . We call the coalescent coefficient of . We now show for a given , its coalescent coefficient can be computed efficiently.

3.2 Efficient computation of coalescent coefficient

We first note in Equation (4), . Thus, for all in , can be written as the same . So,

In Equation (6), . Note that the species tree branch b is implicit in . In , the term depends on specific because the numbers of coalescents along species tree branches depend on . When the number of is large, direct summation of these terms is inefficient. We now show can be computed efficiently using a recurrence.

For the remaining part of this section, we consider a specific . Thus, we omit [e.g. by writing as C1 and as ]. We consider each coalescent (internal node) c in T. We denote the species tree branch on which c occurs as b. We denote the number of coalescent events that are within the subtree of T rooted at c and occur on b as n. We denote the set of species tree branches below (and including) b as . That is, . Here, c itself is considered to be within the subtree rooted at c and thus is included in n. We denote S as the set of species tree branches where the coalescent c may occur. For example, In Figure 1, and . b=, b=, n = 3 (including coalescents r, h and q) in the coalescent history shown in Figure 1a. We denote the list of the numbers of coalescents on each species tree branch as . We define unit vector as the vector where and when . For a specific integer list of length m, we say is a sublist of if for all . We denote the set of all possible sublists of as . For example, suppose . Then a possible sublist is [1,1,1,2].

For a coalescent c, a branch b of T and an integer list , we define:

That is, restricts C1 to the subtree rooted at c where c occurs on b and the numbers of coalescents within the subtree rooted at c is . We impose the following constraints for . If any of the constraints is violated, we have .

(1). only specifies coalescent numbers for species tree branches in (i.e. either at b or within the subtree below : for any .

(2). There is at least one coalescent at b: . Moreover, .

We let b be the root branch of T and c be the root of T. Note that we assume there is at least one gene lineage from each species. Then c must occur on b. We have:

Here, refers to the list of upper lineage counts implied by the compact history . We initialize the computation of for the coalescent c, where there is no coalescent within the subtree of T rooted at c except c (i.e. c is the lowest coalescent in :

If c is not the lowest coalescent, then c may have one or two children that are also internal nodes. Recall that we assume bifurcating gene trees in this section. We have two cases.

Case 1: c1 is the only descendant internal node. Note that if there are still more coalescents along b (i.e. ), c1 must also occur along b. If there is only a single coalescent on b, c1 has to occur in some branch in . Then:

Case 2: both immediate descendants c1 and c2 of c are internal nodes. In this case, we need to split into two parts, one for c1 and one for c2. Then, we may sum over all such partitions and all possible species tree branches where c1 and c2 may occur.

Remark: Equations (8) and (9) can be written in a more compact form because many terms contained in these two equations are zero. This is because many partitions of lead to violations to constraints 1 and 2. The current exposition is a simplification to avoid overly complicated formula.

Example. To illustrate the process of computing the C2 terms, we show how to calculate these terms for the compact history in Figure 1a. We have . Starting from the more recent (i.e. close to the extant time) coalescents,

. Then, is equal to . Here, we omit all terms in C2 calculations that are zero.

3.3 The CompactCH algorithm for computing the gene tree probability

The CompactCH algorithm for computing the gene tree probability is as follows.

Enumerate and compute all coalescent coefficients for T and T as in Equations (8) or (9).

Enumerate all compact coalescent histories for T and T. The enumeration can be performed in a bottom-up way. Compute the terms and for each species tree branch b. Obtain using Equations (6) and (7).

Compute the gene tree probability by summing over all compact coalescent histories using Equation (5).

3.4 Analysis of the CompactCH algorithm

We now analyze the running time of the CompactCH algorithm in Section 3.3. It can be shown that CompactCH runs in exponential time when the number of populations is unbounded. See the Supplemental Materials for more details. We now show that CompactCH runs in polynomial time in the number of gene lineages when the number of populations is fixed to be a constant.

The gene tree probability computation in Equation (5) depends on the number of compact coalescent histories for the given T and T, which in turn depends on the number of gene lineages n and the number of species tree branches m. The main advantage of our new algorithm is that it runs in polynomial time when there are constant number of species. The algorithm is more efficient than the existing algorithms when there are multiple gene lineages per species. Thus, we assume m is fixed to a constant. We first show that is polynomial in terms of n when m is a constant. We then show that coalescent coefficient is polynomial time computable for each compact history.

First, recall that a compact coalescent history consists of a list of the upper lineage counts, one for each species tree branch. For a T with m species tree branches, the length of this list is m – 1, and at each position of the list the upper lineage count is at most n – 1. So the number of choices for each position is n. Therefore, there are at most compact histories. When m is fixed to be a constant, this is polynomial in n. For each compact history, Equation (5) takes O(m) time. So, the total time is .

We now show coalescent coefficient computation in Equation (6) runs in polynomial time for a given compact coalescent history in a bifurcating T when m is a constant. These coalescent coefficients are pre-computed before Equation (5) is evaluated. By Equation (7), we need to show can be computed efficiently. Note that is computed using a recurrence over . Here, c is a coalescent in T, b is a species tree branch and is the list of coalescent counts (i.e. the number of coalescents) for the branches with the subtree of T rooted at b. The number of such is bounded by , which is a polynomial of n when m is a constant. Each can be computed by Equations (8) or (9) in time. So, all coalescent coefficients can be computed in time. Therefore, we have:

Theorem 3.2. The CompactCH algorithm runs in time for a bifurcating T, which is polynomial in n when the number of species in T is a constant.

Remark. The CompactCH algorithm computes the gene tree probability runs in polynomial time in n when m is fixed to be constant. At the first glance, the CompactCH algorithm appears to be very slow: when m  =  3 (i.e. two populations), the running time is . In practice, CompactCH appears to be significantly faster than the STELLS algorithm when m is small. This is because we may overestimate the number of needed steps (e.g. for computing the coalescent coefficients). Table 2 shows that when there are two populations, the CompactCH algorithm can remain practical when n is as large as 200 where the STELLS algorithm becomes very slow. See the Section 5 for empirical performance of this algorithm. It can be shown that STELLS runs in exponential time in the number of gene lineages when there two populations. See the Supplemental Materials for details.

Table 2.

Running time of CompactCH (outside the parenthesis) and the STELLS algorithm (inside parenthesis) for computing gene tree probability for 500 simulated gene trees

mp	1	2	5	10	15	20	50	100
2	<1 (<1)	<1 (<1)	1 (1)	3 (7)	20 (453)	94 (21 516)	21 974 (—)	1 166 613 (—)
3	<1 (<1)	1 (<1)	11 (5)	2062 (18747)	83 272 (—)	— (—)	— (—)	— (—)
4	<1 (1)	1 (1)	1634 (218)	— (—)	— (—)	— (—)	— (—)	— (—)

Results are not given if it takes longer than 15 days. m: number of populations. Columns: g, the number of gene alleles per population. Time: in seconds.

3.5 Multifurcating gene trees

So far, gene trees are assumed to be bifurcating. The original STELLS algorithm (Wu, 2012) and the algorithm in Degnan and Salter (2005) also assume bifurcating gene trees. Multifurcating gene trees, however, may be preferred for gene tree probability computation in some cases. For example, it is possible that some splits in the gene trees do not have sufficient support and thus multifurcating gene trees may be used to allow uncertainty in gene tree topologies. The STELLS algorithm was extended to allow multifurcating gene trees in Wu (2015). The STELLS algorithm is slower for multifurcating trees than bifurcating trees.

Our new algorithm can be extended to compute the gene tree probability for multifurcating gene trees. Compact coalescent histories remain the same as the bifurcating gene tree case. For each species tree branch b, we specify the number of gene tree lineages at the top of b. Note that at a multifurcating gene tree node v of out-degree d, there are d − 1 coalescents at v. For a given compact coalescent history and a species tree branch b, the terms and can all be easily computed as before. This is because these terms only depend on the numbers of gene lineages at specific positions of T, which are fully determined by . These terms do not depend on the topology of T. The main difference of the multifurcating gene tree case is on the coalescent coefficient computation. At a gene tree internal node c with out-degree d, recurrences in Equations (8) and (9) need to be modified. Let D be the set of c’s children in T that are internal nodes. We let be the set of all possible proper subsets of D. We define for , a branch b of T and an integer list in the same way as when we treat the lineages in S form a new split in T. Then, following the same reasoning as in Equations (8) and (9):

Here, and are two integer lists which combine to . Thus, coalescent coefficients for multifurcating T can be computed in a recursive way, similar to the bifurcating case. The running time for coalescent coefficient computation depends on the maximum degree dmax in T. Since there are non-empty proper subsets for dmax gene lineages, the algorithms becomes slow when dmax is large. If dmax is bounded by a constant, the gene tree probability for a multifurcating tree can be computed efficiently when the number of species is small.

4 Inference of population tree from haplotypes from pairwise population distance

In Wu (2015), we demonstrate that gene tree probability can be used in the inference of population tree (i.e. the population divergence history) when the given data is in the form of haplotypes. Briefly, haplotypes contain the alleles (states) at closely linked genetic variation sites. See e.g. Wu (2015) for some background on haplotypes. Assuming population haplotypes satisfy the infinite sites model of mutations (Kimura, 1969) and no intra-locus recombination, we may infer the underlying genealogical tree topologies, although these trees are usually multifurcating [see, e.g. Gusfield, 1991]. For a single locus, there is a unique gene tree topology and there are mutations on the branches of the gene tree. Under the infinite sites model of mutations and with no recombination, the probability of haplotypes is exactly equal to the probability of the unique gene tree with mutations [see, e.g. Wakeley, 2008]. For the sake of computational efficiency, the approach in Wu (2015) uses the probability of the gene tree topology (i.e. ignoring mutations) to approximate the probability of haplotypes. The probability of gene tree topology is easier to compute than that of haplotypes and can be used to infer the population tree by maximum likelihood (Wu, 2015). See Wu (2015) for more details. We use the same approximation in this article. In Wu (2015), we demonstrate that population tree inference with gene tree probability performs well when compared with the program TreeMix (Pickrell and Pritchard, 2012). STELLS infers population trees by searching the tree space to find the maximum likelihood estimate of the population tree. However, one main computational challenge in Wu (2015) is the computational efficiency: population tree inference becomes slow when the numbers of taxa and hapotypes increase (see the Section 5). Genetic studies now routinely involve ten or more populations, and multiple haplotypes are collected from each population. For example, the 1000 Genomes Project (2015) recently released over 1000 genome-scale haplotypes from more than 20 populations. Such large-scale genetic data imposes huge challenges for inference.

In this article, we develop a new population tree inference method based on the CompactCH algorithm. This method is distance based, and does not perform maximum likelihood inference. The main idea of our method is (i) first infer the distance of pairs of two populations using haplotypes from the two populations, and (ii) construct the population tree using neighbor joining from the inferred pairwise population distances. The key is in the population distance estimate. There is exactly one population tree topology (i.e. population divergence history) of two populations A and B. With the haplotypes from A and B (and the implied gene genealogies), we can infer the branch lengths d and d. Here d (respectively d) is the length of the branch connecting A (respectively B) and the common ancestor of A and B in the species tree of A and B. This can be done by optimizing d and d to maximize the gene tree probability of the genealogies. See Wu (2012, 2015) for more details. Once d and d are estimated, the pairwise distance between A and B is estimated to be . Then these pairwise distances are then used to infer the population tree of all populations by neighbor joining. The CompactCH algorithm enables fast pairwise population distance estimate from large number of haplotypes. In the Section 5, we demonstrate that this simple method can give accurate population tree inference using more haplotypes per population than the original method in Wu (2015) can handle.

5 Results

We have implemented our new CompactCH algorithm as part of the STELLS program (http://www.engr.uconn.edu/ywu/STELLS.html). We compare with the STELLS algorithm, which is currently the best algorithm for computing the gene tree probability. We first show that CompactCH outperforms the original STELLS when the number of species is small and there are multiple gene lineages per species. Then, as an application, we show that CompactCH allows fast and accurate inference of population trees from haplotypes when there are multiple haplotypes from each population.

5.1 Performance of CompactCH in gene tree probability computation

We evaluate the efficiency of CompactCH in computing the gene tree probability. We first set the number of populations m to be two, three and four. We randomly generate species trees. We let the number of gene copies g to be 1, 2, 5, 10, 20, 50 and 100. Note the total number of haplotypes n is g. We generate 500 bifurcating gene trees for a species tree using the program ms (Hudson, 2002). The running times of using CompactCH and STELLS to compute the probability of these 500 gene tree topologies for each simulated species tree are shown in Table 2. Both algorithms are very slow on many settings where both m and g are large. So no results are collected on these settings. Here, we say the computation is too slow if the computation does not finish within 15 days on a 3192 MHz Intel Xeon workstation running Linux (with 32 GB memory). Our second simulation aims to investigating the efficiency with larger number of populations. We fix g to be one, and then run CompactCH and STELLS on the data. The results are shown in Table 3. As expected, the STELLS is more efficient than CompactCH for larger number of populations. Overall, CompactCH is much faster than STELLS when m is relatively small and g is relatively large.

Table 3.

Running time of CompactCH (outside the parenthesis) and the STELLS algorithm (inside parenthesis) for computing gene tree probability for 50 simulated gene trees and g = 1 (i.e. 1 allele per population)

2	3	4	5	10	15
<1 (<1)	<1 (<1)	<1 (<1)	<1 (<1)	3 (1)	208 (1)

Columns: number of populations. Time: in seconds.

5.2 Using CompactCH in population tree inference from haplotypes

We run CompactCH and STELLS to infer population trees on simulated data. The design of simulations is similar to that in Wu (2015). Briefly, we generate simulated datasets as follows. We use randomly generated populations trees with m populations to model the population divergence history. Here, we let m to be eight. The population trees are the same as in Wu (2012). The length of each branch in the tree is assigned to a length that is uniformly chosen between 0.0 and 1.0. We then scale the population trees so that the trees satisfy clock property and the trees have fixed heights (being either 0.1 or 0.5). Note that CompactCH does not require clock-like population trees. The clock property is mainly for the ease of simulation, where we can reduce the number of parameter combinations. Then we simulate haplotypes for 100 loci from a given population tree using the program ms Hudson (2002), where there are g gene lineages per species. Mutation parameters are fixed to be 20. We assume constant population size, no recombination within each locus and no migration between populations. Recall that the inferred genealogical tree topologies from haplotypes are usually multifurcating, and CompactCH becomes slow when the gene tree topologies have large degree at internal nodes. Thus, we discard gene trees with node degree larger than 10. We use the normalized Robinson–Foulds (RF) distance to measure the inference error in the topology. That is, for the inferred population tree topology and the true tree with n taxa, the normalized RF distance is . Here, is the number of splits in but not in .

It has been demonstrated in Wu (2015) that the inference of population trees from these inferred genealogical tree topologies can give accurate results. In Table 4, we show that more accurate inference results can be obtained by using the pairwise population distance approach (denoted as CompactCH) in Section 4 than the STELLS approach (denoted as STELLS), when more haplotypes are used. These simulations are conducted on a computer cluster. STELLS is more accurate than CompactCH when only two haplotypes are used for each population. The main advantage is that CompactCH is much more scalable than STELLS: STELLS takes more than 5 h on average for each dataset with two alleles per population. CompactCH only takes seconds in this case. Note that STELLS can be more accurate than the results by CompactCH because STELLS uses the maximum likelihood approach to search in the tree space while CompactCH builds trees from pairwise distances. The accuracy of CompactCH increases when more haplotypes are provided. For example, when 20 or more haplotypes per population are used for the population trees of height 0.1, the inference error of CompactCH is only about one-third of that of STELLS which uses two alleles per population. It appears CompactCH performs relatively better for the species tree of height 0.1 than those of height 0.5. We also note that CompactCH becomes slow if a large number of haplotypes are analyzed. For comparison, we also run TreeMix on these data. Our results indicate that CompactCH outperforms TreeMix in most of the cases, especially for the case of shorter species tree height.

Table 4.

Accuracy and time for inferring population trees using pairwise population distances

Ht	Accuracy/Time	gc
		2	4	10	20	30
0.1	Inf. error	0.38 (0.47) STELLS: 0.30	0.20 (0.34)	0.20 (0.24)	0.11 (0.18)	0.11 (0.16)
	Time	8 s (44 h,18 m, 8 s)	33 s	21 m, 49 s	7 h, 2 m, 4 s	36 h, 42 m, 49 s
0.5	Inf. error	0.23 (0.22) STELLS:0.14	0.20 (0.15)	0.10 (0.12)	0.12 (0.14)	0.12 (0.13)
	Time	6 s (25 h, 30 m, 37 s)	29 s	7 m, 18 s	1 h, 33 m, 5 s	5 h, 56 m, 55 s

Average over 50 replicates. Eight populations. Inference error: normalized RF distance. Population tree height (Ht): 0.1 or 0.5. 100 loci. g: number of haplotypes per population. Time: in seconds (s), minutes (m) and hours (h). Results for TreeMix are inside the parentheses. The original STELLS is only feasible for g = 2 and so only the results for the g = 2 case are provided for the original STELLS.

5.2.1 Inference with the 1000 Genomes Haplotypes

We use CompactCH to infer population trees with haplotypes from the 1000 Genomes Project (2015). We use haplotypes from the following ten populations in the 1000 Genomes Project: Han Chinese in Beijing, China (CHB), Japanese in Tokyo, Japan (JPT), Southern Han Chinese (CHS), Utah Residents with Northern and Western European ancestry (CEU), Toscani in Italia (TSI), Finnish in Finland (FIN), British in England and Scotland (GBR), Iberian population in Spain (IBS), Yoruba in Ibadan, Nigera (YRI), and Luhya in Webuye, Kenya (LWK). We use ten diploid individuals (i.e. 20 haplotypes) for each population. We choose the loci where there are few recombinations, and then construct gene trees from haplotypes within these loci. See Wu (2015) for details on how the loci are picked. We infer one gene tree topology from the haplotypes at each locus. These gene trees are then used to infer the population trees. Here, we discard gene trees that have out-degrees at internal nodes of nine or larger because gene trees with large degrees greatly slow down the computation. We first infer the population tree for four populations: CEU, CHB, JPT and YRI. As expected, CHB and JPT are close siblings and in the population tree, and the ancestral population of CHB and JPT is the sibling of CEU. We then infer the population tree for all ten populations (which takes less than 2 h). The inferred neighbor joining tree is shown in Figure 2. Note that the tree should be viewed as an un-rooted tree. As expected, African, European and Asian populations cluster on the tree. The tree agrees mostly with the result in Wu (2015) with small differences (e.g. the location of GBR).

Fig. 2.

The inferred population tree from ten populations in the 1000 Genomes Project using 20 haplotypes from then individual per population. Branch length shown is the estimated time in coalescent units

6 Discussion and conclusions

In this article, we present the CompactCH algorithm for computing the gene tree probability under the multispecies coalescent model. While CompactCH is much faster than the STELLS when the number of populations is small, computing the exact gene tree probability for more populations remains a challenging task (Table 2). Nonetheless, we show in this article that efficient computation of gene tree probability for small number of (say two) populations can find applications in the inference of population history. We believe the CompactCH algorithm is a step forward for more efficient computation on coalescent models. The key idea of CompactCH is using compact coalescent histories. While merging multiple entities is a common idea in algorithm design, designing a working algorithm for coalescent computation based on this high-level idea is not trivial, as we show in this article. Our work suggests that coalescent computation can indeed be made more efficient by a well-designed algorithm [see Wu (2010) for an algorithm that speeds up coalescent computation under a different formulation].

We note that coalescent histories (in particular the mathematical properties of coalescent histories) have been actively studied recently [see e.g. Rosenberg 2013]. Our result presented here demonstrates that computation based on coalescent history can also be improved computationally. The CompactCH approach may also lead to speedup in coalescent computation in other formulations, where more complex coalescent models are used and computational efficiency is highly desirable. For example, the coalescent likelihood computed by MCMC in Heled and Drummond (2010) considers sequences, not the gene trees inferred from the sequences. Our techniques may be applied to speed up the computation of such likelihood.