Testing Multispecies Coalescent Simulators using Summary Statistics

TL;DR

This paper develops statistical tools to evaluate the accuracy of multispecies coalescent simulators, ensuring they produce valid samples, and provides an R package for easy testing of future simulators.

Contribution

It introduces new summary statistic-based tests for validating MSC simulators and uncovers flaws in existing tools, with implementation in an accessible R package.

Findings

Identified flaws in several published MSC simulators

Developed statistical tests for simulator validation

Provided an R package for easy performance checking

Abstract

As genomic scale datasets motivate research on species tree inference, simulators of the multispecies coalescent (MSC) process are essential for the testing and evaluation of new inference methods. However, the simulators themselves must be tested to ensure they give valid samples from the coalescent process. In this work we develop several statistical tools using summary statistics to evaluate the fit of a simulated gene tree sample to the MSC model. Using these tests on samples from four published simulators, we uncover flaws in several. The tests are implemented as an R package, so that both developers and users will be able to easily check proper performance of future simulators.

Tables

Table 1. Table 1. Rooted triple topology counts from 100,000 gene trees sampled for the species tree of Figure 3

	$((A,B),C)$	$((A,C),B)$	$((B,C),A)$	$p$-value	Int. branch
Expected	59564	20217	20217	-	0.5
Mesquite	66000	17556	16444	0	0.670
Hybrid-lambda	59564	20195	20240	0.975	0.499
Phybase	59395	20256	20349	0.492	0.495
SimPhy	59120	20397	20483	0.554	0.497
	$((A,B),D)$	$((A,D),B)$	$((B,D),A)$	$p$-value	Int. branch
Expected	70044	14977	14977	-	$0.83\overline{3}$
Mesquite	74074	12844	13082	4.286e-99	0.944
Hybrid-lambda	70850	14683	14466	0.207	0.827
Phybase	71072	14450	14478	0.940	0.834
SimPhy	71079	14465	14456	0.934	0.835
	$((A,C),D)$	$((A,D),C)$	$((C,D),A)$	$p$-value	Int. branch
Expected	52231	23884	23884	-	$0.3\overline{3}$
Mesquite	50533	24659	24808	6.16e-26	0.298
Hybrid-lambda	52230	23951	23818	0.830	0.3331
Phybase	52419	23606	23975	0.118	0.337
SimPhy	52141	23834	24025	0.579	0.331
	$((B,C),D)$	$((B,D),C)$	$((C,D),B)$	$p$-value	Int. branch
Expected	52231	23884	23884	-	$0.3\overline{3}$
Mesquite	51281	24389	24330	1.32e-08	0.313
Hybrid-lambda	52241	23960	23798	0.758	0.3335
Phybase	52164	23825	24011	0.635	0.331
SimPhy	52220	23889	23891	0.997	0.3330

Equations

$$((a\text{:}{\ell }_a,b\text{:}{\ell }_b)\text{:}{\ell }_v,c\text{:}{\ell }_c)$$

$$\mathbb{P}(((a,b),c))$$

$$\mathbb{P}(((a,c),b))$$

$$\mathbb{P}(((b,c),a))$$

$$Y=g_a+g_b+2X,$$

$$N^{\ast }(x)=N_{e_i}(x-m_{i-1})$$

$$c(x)$$

$$=\left(\prod _{j=1}^{i-1}{\eta }_j\right)\frac{\exp \left(-{\int }_0^{x-m_{i-1}}\frac{1}{N_{e_i}(\tau )}d\tau \right)}{N_{e_i}(x-m_{i-1})},$$

$${\eta }_i=\exp \left(-{\int }_0^{{\ell }_i}\frac{1}{N_{e_i}(\tau )}d\tau \right)$$

$$f(y)=\{\begin{array}{ll}0& \text{ for }y\le g_{ab},\\ & \\ \left(\prod _{j=1}^{i-1}{\eta }_j\right)\frac{\exp \left(-{\int }_0^{\frac{y-g_{ab}-2m_{i-1}}{2}}\frac{1}{N_{e_i}(\tau )}d\tau \right)}{2N_{e_i}(\frac{y-g_{ab}-2m_{i-1}}{2})}& \\ & \text{for }\{\begin{array}{l}{g}_{ab}+2m_{i-1}\le y<g_{ab}+2m_i,\\ 1\le i\le k,\end{array}\\ & \\ \left(\prod _{j=1}^k{\eta }_j\right)\frac{\exp \left(-{\int }_0^{\frac{y-g_{ab}-2m_k}{2}}\frac{1}{N_r(\tau )}d\tau \right)}{2N_{e_r}(\frac{y-g_{ab}-2m_k}{2})}& \\ & \text{ for }{g}_{ab}+2m_k\le y.\end{array}$$

$$P=(e_1,e_2,\dots e_i)$$

$$\mathbb{P}\left(((a,c),b)\right)=\mathbb{P}\left(((b,c),a)\right)=\frac{1}{3}\prod _{j=1}^i{\eta }_j.$$

$$\mathbb{P}\left(((a,b),c)\right)=1-\frac{2}{3}\prod _{j=1}^i{\eta }_j.$$

Full text

1. Introduction

With the increasing availability of genomic-scale datasets, comprised of sequences from many genetic loci, it has become clear that individual gene trees inferred for the same taxa are often discordant. While gene tree inference error may be a factor, there are also many biological processes that could lead to this discordance. Of these, incomplete lineage sorting as described by the multispecies coalescent (MSC) model is often viewed as the most fundamental. The MSC may thus be thought of as the null model of discordance, to be considered before further complications such as hybridization, lateral gene transfer, gene duplication and loss, and/or population structure, are invoked (Degnan, 2018). While inference of species trees under the MSC is still challenging, it can now be performed in a variety of ways (Liu, 2008; Heled and Drummond, 2010; Chifman and Kubatko, 2014; Vachaspati and Warnow, 2015; Zhang et al., 2018).

Evaluating the performance of any approach to species tree inference under the MSC, though, requires testing on simulated datasets, so that the true species tree is known and accuracy can be judged. Thus MSC software for simulating the formation of gene trees within a species tree or network, often paired with software for simulating the evolution of sequences along the gene trees, plays a critical role in advancing methodology. Unfortunately, as we demonstrate below, some of the available simulators do not produce valid MSC samples, although this may not be apparent even to knowledgeable users.

Testing a coalescent simulator for correctness is not simple, as the theoretical distribution of gene trees gives positive probability to all gene tree topologies, with the probability density for metric gene trees having a quite complicated dependence on the species tree branch lengths. “Straightforward” comparison of a sample to the theoretical distribution is simply not a practical approach. This has led some developers to validate their software by comparing output to other simulators, rather than to theoretical predictions (Mallo et al., 2016). Here we introduce several testing tools, based on summary statistics that capture either topological or metric information. Implemented in an R package, MSCsimtester, these can be applied to gene trees samples from any simulator to study whether its output is in accord with the MSC. Although examining such summary statistics cannot give an ironclad guarantee of correctness, we believe they are likely to uncover most problems.

When we applied the MSCsimtester tools to output from four well-known MSC simulators, we discovered that only one of the four, SimPhy (Mallo et al., 2016), behaved as we expected. The documentation for one, Phybase (Liu and Yu, 2010), was not sufficiently explicit on specifying its input and our initial interpretation was incorrect. A third, Hybrid-lambda (Zhu et al., 2015) passes our tests for gene tree topologies, but samples metric gene trees incorrectly. The fourth, Mesquite (Maddison and Maddison, 2018), produced samples with neither gene tree topologies nor metric properties in accord with the model. However, we found website documentation of only the Hybrid-lambda issue (which we had uncovered in a preliminary version of this work).

While we notified the authors of these simulators of the problems in advance of this publication, the larger community should be aware of the need to interpret results of previous simulation work with them cautiously. We also suggest that users of other simulators, and developers of new ones, test them with the MSCsimtester tools, which are available at the jarhodesuaf.github.io/software.html website.

2. New Approaches

We give an informal description of the MSC to introduce the summary statistics we focus upon.

Suppose first that 3 taxa are related by the species tree with topology $((a,b),c)$. Acknowledging that species are composed of populations, we depict this by a tree whose edges are ‘pipes’ as in Figure 1. The length of each pipe is elapsed time measured in generations, and the width of the pipe represents population size, which may vary over time and edge. When individual genes are sampled from the leaves of the species tree, they trace backwards in time within the species tree until they coalesce at a common ancestral individual. Coalescence is a random process, which can be thought of as individual genes on the lineage choosing their ‘parent’ uniformly at random from those existing in the previous generation, a panmictic viewpoint. Thus the only population detail of importance under the MSC is size. Importantly, there is a greater chance of coalescence when populations are small, but no requirement that lineages coalesce within any specified finite time.

Often, simplifying assumptions on population sizes are made by modelers and programmers, such as all populations at all times and on all branches throughout the tree are a constant $N$. More realistic is to at least allow different population sizes $N_{e}$ for each edge (pipe) $e$ of the species tree. While it would be highly desirable to be able to simulate gene trees under the MSC using arbitrary population size functions, current simulators do not makes this easy on a large tree. Nonetheless, our testing methods accommodate that generality.

When gene trees are produced under the MSC they have two aspects that can be viewed somewhat separately. One is the metric information, which is reflected in the distribution of pairwise distances between two fixed taxa across the gene trees. The second is topology, which is reflected in the distribution of rooted triple trees on three fixed taxa displayed on the gene trees. In testing the performance of a simulator, it is important both metric and topological aspects be examined.

2.1. The distribution of pairwise distances on gene trees

Consider now the species tree

[TABLE]

of Figure 1, with root $r$ and $v$ the most recent common ancestor of $a$ and $b$. Fix constant population sizes $N_{v}$ on the edge above $v$, and $N_{r}$ above the root $r$. For the purpose of illustration, assume $N_{r}$ is smaller than $N_{v}$. Tracking two genes $A$ and $B$ sampled from $a$ and $b$ backwards, their lineages cannot coalesce until the reach the population above $v$. On the branch above $v$ coalescence occurs by a Poisson process at a constant rate $1/N_{v}$. This leads to the time to coalescence above $v$ being exponentially distributed, at least for those times more recent than the root of the tree. At the root of the tree there is a discontinuous change in the coalescence rate, to the larger value $1/N_{r}$. Coalescence above the root $r$ again occurs with an exponential distribution, but at a faster rate. Putting together the 3 regions analyzed here (below $v$, between $v$ and $r$, above $r$, we find the distribution of distances between $a$ and $b$ on gene trees is a piecewise exponential, such as shown in Figure 2.

Note that the discontinuities in the distribution in Figure 2 occur due to 1) the impossibility of coalescence below $v$ and 2) the discontinuity in population size at $r$. For larger trees, with more edges leading from the MRCA of a pair of taxa to the root, there can of course be more discontinuities in the density even when the population size is constant on each edge. If population sizes vary on edges, the pieces need not come from exponential distributions, but are computable from the population size function. Discontinuities in the distribution arise whenever there is a discontinuity in the population size. Finally none of this depends on the species tree being ultrametric, which of course it need not be biologically since edge lengths are in generations. The precise form of the pairwise distance distribution, and its derivation, is given in the Methods section.

To test an MSC simulator using this distribution, one can produce a large sample of gene trees from a fixed species tree with population sizes, and then choosing some pair of taxa compare a histogram of the pairwise distances between these taxa across the gene trees to the theoretical distribution. This comparison can be done visually, as major deviations from the theoretical predictions will be obvious. One can further perform a statistical test, such as that developed by Anderson and Darling (1952), to compare the empirical distribution from the simulator to the theoretical one, giving a $p$-value.

2.2. The distribution of rooted triple topologies on gene trees

Returning to the species tree of Figure 1, we can ignore distances on gene trees and instead examine topological features. Supposing again we have constant population sizes $N_{v}$ and $N_{r}$, the chance that lineages from $a$ and $b$ fail to coalesce in the edge of length $\ell_{v}$ between $v$ and $r$ can be computed to be $e^{-x},$ where $x=\ell_{v}/N_{v}$ (Pamilo and Nei, 1988). Note that a longer edge (increasing $\ell_{v}$) and/or a smaller population (decreasing $N_{v}$) make coalescence more likely. The quantity $x$ here is the length of the edge in coalescent units, a convenient unit to address the confounded effects of population size and time. If the lineages fail to coalesce before the root, then lineages from $a,b,c$ will all be present above $r$, and all three rooted gene topologies $((a,b),c)$, $((a,c),b)$, and $((b,c),a)$ are equally likely to form. This leads to the gene tree probabilities

[TABLE]

which were derived by Pamilo and Nei (1988). More general formulae, accommodating changing population sizes, are given in the Methods section.

Given a large sample of gene trees from an MSC simulator, for any 3 fixed taxa one can tabulate the frequencies of the three rooted topologies displayed on the gene trees. The simulation can then be tested in two ways: First, one may find the maximum likelihood estimator of $x$ from the tabulated frequencies, and compare this to the theoretical value. Second, one can judge the fit of the empirical frequencies to the expected ones by calculating, for instance, a chi-squared statistic and judging it against the $\chi^{2}$ distribution with 2 degrees of freedom to determine a $p$-value. The second of these is sensitive to imbalances between the counts for the two topologies incongruent with the species tree, while the first is not.

3. Results and Discussion

We focused on simulators that 1) allowed easy input of a species tree either in Newick notation, or graphically, with branch lengths in generations and 2) allowed a population size to be assigned to each branch in the species tree independently. The four simulators we investigated were Mesquite (Maddison and Maddison, 2018), Phybase (Liu and Yu, 2010), Hybrid-lambda (Zhu et al., 2015), and SimPhy (Mallo et al., 2016). Our criteria ruled out several other well-known simulators, including ms (Hudson, 2002) (which does not provide for easy input of a species tree in Newick), and PhyloNet (Than et al., 2008) (which calls ms and allows for a Newick tree, but does not allow different population sizes on different branches).

Mesquite is a popular package with many phylogenetic tools accessible through a graphical interface. Phybase is an R package which also includes phylogenetic tools beyond coalescent simulation. Hybrid-lambda implements the standard coalescent model, and the more general $\lambda$-coalescent, on both trees and hybridization networks. SimPhy includes a very general model of gene tree generation, including gene duplication and loss. Our tests, therefore, should not be construed as complete tests of any of this software; we test only a standard MSC simulation on a species tree, with branch lengths in generations.

We performed tests with a number of species trees, with up to 6 taxa, using a variety of constant population sizes. Here we show only representative examples of this work, using the species tree and populations depicted in Figure 3. Additional test results are shown in the supplementary materials, in Figures S2-S21 and Tables S1-S4.

A sample of 100,000 gene trees was simulated with each of the programs, the values of $d(A,B)$ on them was extracted, and a histogram produced. These are shown in Figure 4, with the theoretical distribution superimposed on them. There is a good match for SymPhy, as was seen in all our simulations with this software. When our initial Phybase simulation did not match expectations, we learned that species tree branch lengths should be supplied as $\mu t$, where $t$ is in generations and $\mu$ is a mutation rate, while population sizes should be specified as $\theta=4\mu N_{d}$ where $N_{d}$ is diploid population size. Taking $\mu=1$ and $\theta=2N$ with $N$ the haploid population size, Phybase’s sample matched expectations well. Both the Hybrid-lambda and Mesquite simulations show pronounced deviations from the theoretical distribution. However, when the input tree is specified in coalescent units, Hybrid-lambda does give a sample of gene trees whose pairwise distances in coalescent units match the theory well; the poor fit occurs only when the species tree parameters involve separate time and population values.

To quantify deviation of the sampled $d(A,B)$ distribution from the expected distribution, we apply the Anderson-Darling test. Note that even small numerical errors in a simulation or its analysis may prevent the extremely close fit to theory that a large sample should exhibit. As a result, with very large samples such tests can produce misleadingly small $p$-values, leading to excessive rejection. We therefore divided each of our samples into 100 subsamples of size 1000, computing a $p$-value for each. A good fit is shown by a roughly uniform distribution of $p$-values for the subsamples. Figure 5 shows these $p$-value distributions for the samples from SimPhy and Mesquite, formally confirming the conclusions already described.

Topological features of samples were analyzed by tabulating counts of all rooted triple topologies across the sampled gene trees, and then performing both a chi-squared test and computing the MLE of the internal branch length on the species tree triple. Results are shown in Table 1. For all programs except Mesquite, no $p$-values were extreme enough to suggest poor fit. However, $p$-values for Mesquite are strongly suggestive of poor fit, being quite close to 0. Internal branch length estimates were also poorest for the Mesquite sample. Note that while Hybrid-lambda had poor metric performance, by this testing procedure it gave a good topological sample. This behavior is consistent with its pairwise distance performance when parameters are given in coalescent units; if the Hybrid-lambda algorithm is based on coalescent units, errors may occur in conversion into numbers of generations.

These results indicate that inadequate attention has been given previously to ensuring MSC simulators perform correctly. The tests based on metric and topological summary statistics implemented in the R package MSCsimtester can uncover errors in simulators, and in user input to simulators. We recommend these tests be routinely used by developers of such simulators and, until software has been fully vetted, by anyone performing coalescent simulations.

4. Methods: Derivations of summary distributions

Let $(S,\{\ell_{e}\},\{N_{e}\})$ be a metric species tree with population size functions, where each edge $e$ has length $\ell_{e}$ and population size $N_{e}:[0,\ell_{e})\to\mathbb{R}^{>0}$. Here $N_{e}(t)$ denotes the population size for a haploid organism $t$ generations above the child node of $e$. There is also an ‘above the root’ population size function $N_{r}:[0,\infty)$. (For diploid taxa, the population sizes should be doubled). For technical reasons, we assume $1/N_{e}(t)$ is integrable on finite intervals.

4.1. Pairwise distance distribution

Let $v$ be the most recent common ancestor of taxa $a$ and $b$ (that is, the node on $S$ where $A$ and $B$ lineages enter the same population for the first time), and let $P_{a}$ be the path in $S$ from $a$ to $v$, $P_{b}$ be the path in $S$ from $b$ to $v$, and $P_{v}$ be the path in $S$ from $v$ to the root $r$. Then $P_{v}=(e_{1},e_{2},...,e_{k})$, where $v$ is incident to $e_{1}$ and $r$ is incident to $e_{k}$. Finally, let $g_{a}=\sum_{e\in P_{a}}\ell_{a}$ and $g_{b}=\sum_{e\in P_{b}}\ell_{b}$ be the number of generations from $a$ and $b$, respectively, to $v$. Then the distance $d(A,B)$ is a random variable

[TABLE]

where $X$ is the random variable giving the time to coalescence of two lineages at $v$. Let $c(x)$ be the probability density function for $X$.

To compute $c(x)$, let $N^{*}:[0,\infty)\to\mathbb{R}^{>0}$ be the piecewise ‘union’ of the $N_{e}$ for $e\in P$ and $N_{r}$, which with $m_{0}=0$, $m_{j}=\sum_{i=1}^{j}\ell_{e_{i}}$ for $1\leq j\leq k$, $m_{k+1}=\infty$, and $N_{k+1}$ the population function ancestral to the root is given by

[TABLE]

for $x\in[m_{i-1},m_{i})$, $1\leq i\leq k+1$. Since the coalescent process for two lineages in the same population of size $N^{*}(x)$ occurs with instantaneous rate $1/N^{*}(x)$, the probability density function is (Allman et al., 2019)

[TABLE]

for $x\in[m_{i-1},m_{i})$, where

[TABLE]

is the probability that 2 lineages entering edge $e_{i}$ fail to coalesce on it.

Since $X=\frac{Y-g_{a}-g_{b}}{2}$, setting $g_{ab}=g_{a}+g_{b}$ this shows the probability density function for $Y$ is

[TABLE]

In the special case that a population size function $N_{e}(t)$ is constant, this shows that the corresponding piece of $f$ is a shifted, scaled, and possibly truncated exponential density.

4.2. Rooted triple frequencies

Suppose the rooted triple $((a,b),c)$ is displayed on $S$, and let

[TABLE]

denote the path from the most recent common ancestor of $a,b$ on $S$ to the most recent common ancestor of $a,b,c$. With the notation of the previous subsection, the probability that the $a$ and $b$ lineages fail to coalesce within $P$ is $\prod_{j=1}^{i}\eta_{j}.$ Note that the gene triplets $((a,c),b)$ and $((b,c),a))$ can only form if the $a,b$ lineages do not coalesce on $P$. Moreover, if $a,b$ have not coalesced on $P$ then by the exchangability of lineages in the same populations under the MSC, the probability that any particular pair of $a,b,c$ coalesce first is 1/3. Thus

[TABLE]

Since the probabilities of the three possible topologies sum to 1,

[TABLE]

In the special case of constant population sizes $\eta_{j}=\exp(-\ell_{e_{i}}/N_{e_{i}})$. More generally, the length of $e_{i}$ in coalescent units is $\int_{0}^{\ell_{i}}\frac{1}{N_{e_{i}}(\tau)}d\tau$, and $\prod_{j=1}^{i}\eta_{i}=\exp(-x)$ where $x$ is the length of $P$ in coalescent units.

5. Acknowledgements

This research was supported by the National Institutes of Health Grant R01 GM117590, awarded under the Joint DMS/NIGMS Initiative to Support Research at the Interface of the Biological and Mathematical Sciences.