Ancestral inference from haplotypes and mutations

TL;DR

This paper develops methods for inferring the evolutionary history of DNA sequences based on haplotype and mutation data, using coalescent theory and sampling techniques, with applications to human Y chromosome data.

Contribution

It introduces new theoretical results and sampling algorithms for ancestral inference from haplotypes and segregating sites, extending existing formulas and accommodating variable population sizes.

Findings

Effective inference methods demonstrated on human Y chromosome data

Extension of Ewens Sampling Formula for haplotype configurations

Implementation of rejection and importance sampling schemes

Abstract

We consider inference about the history of a sample of DNA sequences, conditional upon the haplotype counts and the number of segregating sites observed at the present time. After deriving some theoretical results in the coalescent setting, we implement rejection sampling and importance sampling schemes to perform the inference. The importance sampling scheme addresses an extension of the Ewens Sampling Formula for a configuration of haplotypes and the number of segregating sites in the sample. The implementations include both constant and variable population size models. The methods are illustrated by two human Y chromosome data sets.

Tables

Table 1. Table 1: Result of 10,000 runs for the constant population size coalescent model with a sample of size n = 1 , 544 𝑛 1 544 n=1,544 and s = 9 𝑠 9 s=9 segregating sites. Table shows average value of A n ​ ( t ) , S n ​ ( t ) subscript 𝐴 𝑛 𝑡 subscript 𝑆 𝑛 𝑡 A_{n}(t),S_{n}(t) given S n = 9 subscript 𝑆 𝑛 9 S_{n}=9 . Righthand columns give SE of the mean.

$t$	$S_n\left(t\right)$	SE	$A_n\left(t\right)$	SE
0	9.0		1544
0.1	3.18	0.031	19.47	0.051
0.5	1.11	0.023	3.72	0.023
1.0	0.38	0.015	1.78	0.015
1.5	0.12	0.009	1.25	0.009

Table 2. Table 2: Result of 1,000,000 runs for the constant population size coalescent model with a sample of size n = 1 , 544 𝑛 1 544 n=1,544 and 9 9 9 segregating sites, and haplotype frequencies given above. Table shows average values of K n ​ ( t ) , A n ​ ( t ) , S n ​ ( t ) subscript 𝐾 𝑛 𝑡 subscript 𝐴 𝑛 𝑡 subscript 𝑆 𝑛 𝑡 K_{n}(t),A_{n}(t),S_{n}(t) conditional on the haplotype frequencies and S n = 9 subscript 𝑆 𝑛 9 S_{n}=9 .

$t$	$K_n\left(t\right)$	$S_n\left(t\right)$	$A_n\left(t\right)$
0	10	9	1544
0.1	5.09	4.09	19.9
0.5	1.84	0.85	3.94
1.0	0.74	0.17	1.75
1.5	0.21	0.03	1.19

Table 3. Table 3: Relative errors for estimates from Tables 1 and 2

$t$	$S_n\left(t\right)$	$A_n\left(t\right)$
0.1	0.25	0.02
0.5	0.75	0.06
1.0	0.75	0.02
1.5	1.2	0.05

Table 4. Table 4: Extant haplotype counts for the Hammer data at different times in the past.

time $t$	Haplotype frequency
0	21	23	853	188	75	1	68	31	67	217
0.1	0.227	0.250	11.7	2.30	0.862	0.010	0.777	0.340	0.765	2.69
0.5	0.033	0.036	2.63	0.372	0.130	0.002	0.117	0.050	0.115	0.441
1.0	0.015	0.011	0.837	0.140	0.050	0.001	0.045	0.020	0.044	0.165
1.5	0.006	0.003	0.206	0.043	0.016	0.000	0.014	0.006	0.014	0.051

Table 5. Table 5: Average number of ancestral lineages in the Hammer data at time t 𝑡 t in the past.

t	12	13	14	15	16	17	18	19	20
0.1	0.001	0.003	0.009	0.023	0.048	0.083	0.121	0.149	0.159
	21	22	23	24	25	26	27	28
	0.142	0.111	0.074	0.044	0.022	0.010	0.004	0.001
0.5	1	2	3	4	5	6	7	8
	0.011	0.088	0.259	0.337	0.216	0.073	0.014	0.002
1.0	1	2	3	4	5
	0.426	0.414	0.143	0.016	0.001
1.5	1	2	3
	0.826	0.163	0.011

Table 6. Table 6: Allele multiplicities observed in the TBL1Y dataset

${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\alpha }_4$	${\alpha }_5$	${\alpha }_6$	${\alpha }_7$	${\alpha }_{14}$	${\alpha }_{32}$	${\alpha }_{50}$	${\alpha }_{61}$
107	12	6	1	1	2	1	1	1	1	1

Table 7. Table 7: Likelihoods with expansion in the TBL1Y data.

$\beta$	Replicate 1	Replicate 2	Average
0	2.0934e-61	3.4592e-62	1.2197e-61
0.5	1.2203e-60	1.0782e-61	6.6405e-61
1.0	4.1350e-60	1.1497e-60	2.6424e-60
1.5	9.2120e-61	2.8962e-61	6.0541e-61
2.0	1.5297e-61	5.9120e-61	3.7209e-61
2.5	4.9886e-62	3.7215e-62	4.3551e-62

Table 8. Table 8: Allele age within groups.

		Haplotype groups
${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\alpha }_4$	${\alpha }_5$	${\alpha }_6$
0.0123	0.0086	0.0222	0.0249	0.0421	0.0390
${\alpha }_7$	${\alpha }_{14}$	${\alpha }_{32}$	${\alpha }_{50}$	${\alpha }_{61}$
0.0290	0.1338	0.1331	0.1259	0.1304

Equations

$$\mathbb{P}(A_n^{\theta }(t)=k)$$

$$\mathbb{E}\big{[}A_{n}^{\theta}(t)_{[r]}\big{]}=\sum_{k=r}^{\infty}\rho_{k}^{\theta}(t)(2k+\theta-1){k-1\choose r-1}(\theta+k)_{(r-1)}\frac{n_{[k]}}{(n+\theta)_{(k)}};$$

$$f_{nl}^{\theta }(t)=\frac{2}{l(l+\theta -1)}\mathbb{P}(A_n^{\theta }(t)=l).$$

$$H_n(z)$$

$$\mathbb{P}(S_n=k)=(n-1){\int }_0^1(1-x{)}^{n-2}{e}^{-(-\theta \log x)}\frac{(-\theta \log x{)}^k}{k!}dx,$$

$$\mathbb{P}(S_{n}=k)=\frac{n-1}{\theta}\sum_{l=1}^{n-1}(-1)^{l-1}{n-2\choose l-1}\Bigg{(}\frac{\theta}{l+\theta}\Bigg{)}^{k+\theta}.$$

$$n(n+\theta -1)\mathbb{P}(S_n=s)=\theta \mathbb{P}(S_n=s-1)+n(n-1)\mathbb{P}(S_{n-1}=s),n\ge 2,$$

$$\mathbb{P}(K_n=k)={\theta }^k|S_k^n|/{\theta }_{(n)},1\le k\le n,$$

$$(s,l)\to \{\begin{array}{ll}(s+1,l)& \text{at rate }\theta l/2\\ (s,l-1)& \text{at rate }l(l-1)/2\end{array}.$$

$$\prod _{j=k}^n{e}^{\theta (z-1)j{t}_j/2}\left(\genfrac{}{}{0ex}{}{j}{2}\right)e^{-\left(\genfrac{}{}{0ex}{}{j}{2}\right)t_j/2}$$

$$=\frac{{\prod }_{j=k}^n\left(\genfrac{}{}{0ex}{}{j}{2}\right)}{{\prod }_{j=k}^n\frac{1}{2}j(j+\theta (1-z))}\cdot \prod _{j=k}^n\frac{1}{2}j(j+\theta (1-z))e^{-\frac{1}{2}j(j-1+\theta (1-z))t_j}.$$

$$\frac{{\prod }_{j=k}^n\left(\genfrac{}{}{0ex}{}{j}{2}\right)}{{\prod }_{j=k}^n\frac{1}{2}j(j+\theta (1-z))}\cdot \mathbb{P}(T_k^{\theta (1-z)}+\cdots +T_n^{\theta (1-z)}\le t)$$

$$=\frac{{\prod }_{j=k}^n\left(\genfrac{}{}{0ex}{}{j}{2}\right)}{{\prod }_{j=k}^n\frac{1}{2}j(j+\theta (1-z))}\mathbb{P}(A_n^{\theta (1-z)}(t)\le k)$$

$$\displaystyle=\prod_{j=k}^{n-1}\Big{(}1+\frac{\theta(1-z)}{j}\Big{)}^{-1}\mathbb{P}(A_{n}^{\theta(1-z)}(t)\leq k).$$

$$G_l(z;t)$$

$$\prod_{j=l-1}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}.$$

$$\frac{l(l-1)}{2}\int_{0}^{\infty}G_{l}(z;t)dt=\frac{l-1}{(l-1+\theta_{z})}\cdot\prod_{j=l}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}=\prod_{j=l-1}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}.$$

$$\prod_{j=l-1}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}=\frac{\Gamma(n)}{\Gamma(l-1)\Gamma(n-l+1)}\cdot\frac{\Gamma(l-1+\theta_{z})\Gamma(n-l+1)}{\Gamma(n+\theta_{z})}$$

$$B(l-1,n-l+1{)}^{-1}{\int }_0^1{x}^{l-2}(1-x{)}^{n-l}{e}^{-(-\theta \log x)}\frac{(-\theta \log x{)}^k}{k!}dx,$$

$$G^{*}_{l}(z,t)=\prod_{j=2}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}\mathbb{P}(A_{n}^{\theta_{z}}(t)=l).$$

$$\prod_{j=2}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1},$$

$$\mathbb{P}(A_n^0(t)=l).$$

$$\displaystyle\prod_{j=2}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}\mathbb{E}\big{[}A_{n}^{\theta_{z}}(t)\big{]}$$

$$\displaystyle~{}=\prod_{j=2}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}\sum_{k=1}^{n}\rho_{k}^{\theta_{z}}(t)(2k+\theta_{z}-1)\frac{n_{[k]}}{(n+\theta_{z})_{(k)}}.$$

$$M(k,z)=\prod_{j=2}^{n-1}\Big{(}1+\frac{\theta_{z}}{j}\Big{)}^{-1}(n+\theta_{z})_{(k)}^{-1}.$$

$$(n+\theta +k-1-\theta z)M(k,z)=M(k-1,z),$$

$$(n+\theta +k-1)a(r,k)=\theta a(r-1,k)+a(r,k-1).$$

$$b(r,k)=(2k+\theta -1)a(r,k)-\theta a(r,k-1)$$

$$c(r,k)=\sum _{m=0}^r{e}^{-k\theta t/2}\frac{(k\theta t/2{)}^m}{m!}b(r-m,k).$$

$$\mathbb{E}\big{[}A_{n}(t)\mid S_{n}=r\big{]}=\frac{\sum_{k=1}^{n}\rho_{k}(t)c(r,k)n_{[k]}}{\mathbb{P}\big{(}S_{n}=r\big{)}}.$$

$$(s,k,b,a^{\theta})\to\begin{cases}(s+1,k,b,a^{\theta})&\text{at~{}rate~{}}\theta b/2\\ (s+1,k+1,b+1,a^{\theta}-1)&\text{at~{}rate~{}}\theta a^{\theta}/2\\ (s,k,b-1,a^{\theta})&\text{at~{}rate~{}}\big{(}b(b-1)+2ba^{\theta}\big{)}/2\\ (s,k,b,a^{\theta}-1)&\text{at~{}rate~{}}a^{\theta}(a^{\theta}-1)/2\end{cases}$$

$$(a,a^{\theta })\to \{\begin{array}{ll}(a-1,a^{\theta })& \text{at rate }(b(b-1)+2b{a}^{\theta })/2\\ (a-1,a^{\theta }-1)& \text{at rate }{a}^{\theta }(a^{\theta }-1)/2\\ (a,a^{\theta }-1)& \text{at rate }{a}^{\theta }\theta /2\end{array}$$

Full text

We dedicate this paper to the memory of Paul Joyce, friend and collaborator.

1 Introduction

In this paper we study aspects of the ancestral history of a random sample of $n$ DNA sequences, conditional on features of the haplotype configuration obtained at the present time, labeled 0. Initially, we assume an infinitely-many-sites mutation model with constant population size where the sites are in a completely linked region of DNA where there is no recombination. We begin with some theory that describes the effects of mutations between time 0 and time $t$ in the past. We describe the distribution of the quantities $\widetilde{S}_{n}(t),A_{n}(t),\widetilde{K}_{n}(t)$ and $A_{n}^{\theta}(t)$, where $\widetilde{S}_{n}(t)$ is the number of mutations that have arisen in $(0,t)$, $A_{n}(t)$ is the number of ancestors at time $t$, $\widetilde{K}_{n}(t)$ is the number of distinct haplotypes present in those ancestors that are formed by the mutations in $(0,t)$, and $A_{n}^{\theta}(t)$ is the number of ancestors at time $t$ whose descendants have no further mutation in $(0,t)$.

Results such as these describe the effects of new variation, that arising from the present time 0 to time $t$ in the past. We then describe the effects of standing variation by providing a simulation approach for studying the joint distribution of $S_{n}$ and $K_{n}$, the number of segregating sites and haplotypes in the sample at time 0 resp., and $S_{n}(t),K_{n}(t)$, the number of segregating sites and haplotypes, resp., in the $A_{n}(t)$ ancestors formed by standing variation arising after time $t$ in the past. This approach allows for variable population size, as well as essentially arbitrary binary branching models such as the Yule process.

If the full haplotype configuration and the number of segregating sites is available, then there are many other questions that may be asked about the ancestral history of the current sample. A sequential importance sampling algorithm for studying the ancestral history of a sample of genes conditional on the haplotype configuration and the number of segregating sites is developed. If the complete pattern of mutations on haplotypes were known, then a perfect phylogeny, a genetree, could be constructed and ancestral inference, such as the ages of mutations and time to the most recent common ancestor made conditional on the genetree topology. An example appears in the analysis in [19].

Here we consider the case when just the haplotype frequencies and the number of mutations are known. The theory and computational algorithms are then much simpler. The implementation computes, inter alia, the probability of a sample configuration of haplotypes and segregating sites in a stationary population (this is an extension of the Ewens Sampling Formula), the average coalescence times, mutation times, allele loss times and allele ages back in time conditional on the current haplotype configuration and number of segregating sites, and the conditional average allele configuration and distribution of ancestors at time $t$ in the past. In the infinitely-many-sites model the haplotype configuration follows the infinitely-many-alleles model so the haplotype configuration in a sample is the same as the allele configuration and has a probability distribution of the Ewens Sampling Formula.

2 Ancestral distributions in the coalescent: theory

We begin by setting some notation. We let $A_{n}(t)$ denote the number of ancestors of a sample of size $n$ time $t$ ago, and $A^{\theta}_{n}(t)$ the number of those ancestors whose descendants in $(0,t)$ have no further mutation. We let $\widetilde{S}_{n}(t)$ be the number of segregating sites that arise as mutations in $(0,t)$, $\widetilde{K}_{n}(t)$ be the number of haplotypes in these ancestors that arise from mutations in $(0,t)$, and $S_{n}(t)$ the number of mutations in the sample arising after time $t$ .

Coalescence times between events are denoted by $T_{2},T_{3},\ldots$ In the constant size coalescent model these are independent exponential random variables with rates ${2\choose 2},{3\choose 2},\ldots$ [28]. When considering non-mutant lineages back in time it is appropriate to have lineages lost by either mutation or coalescence. Times between these events are denoted by $T_{2}^{\theta},T_{3}^{\theta},\ldots$. In the constant population size case these are independent exponential random variables with rates ${2\choose 2}+2\theta,{3\choose 2}+3\theta,\ldots$, $\theta$ being the scaled mutation rate. We denote by $f^{\theta}_{nk}(t)$ the density of $T_{n}^{\theta}+\cdots+T_{k}^{\theta}$.

The reader is referred to Ewens [7] for the distribution of haplotypes in a sample, the Ewens Sampling Formula; Griffiths [9], eq. (9), for an original derivation of the distribution of $A_{\infty}^{\theta}(t)$ from looking back in time in a diffusion process; Tavaré [34] for an original derivation of $A_{n}^{\theta}(t)$, connections with the Kingman coalescent and review of ancestral lineage distributions; Tavaré [34, 35] for an introduction to ancestral lineage distributions and ancestral inference in population genetics; and Griffiths [15] for a review and new representations for ancestral distributions.

The distribution of $A_{n}^{\theta}(t)$ is given by

[TABLE]

where $\rho_{j}^{\theta}(t)=e^{-\frac{1}{2}j(j+\theta-1)t}$; see [34], eq. (5.2). $\mathbb{P}(A_{n}(t)=k)$ is given by setting $\theta=0$. The formula also holds for $n=\infty$, where the interpretation is that of the whole infinite population coalescent.

The falling factorial moments of $A_{n}^{\theta}(t)$ are

[TABLE]

see [34], p.13. A simple rate argument establishes the correspondence between the distribution of coalescent times and the number of ancestors, namely

[TABLE]

In a stationary population the probability generating function (pgf) of the number of segregating sites in a sample of $n$ genes is

[TABLE]

Therefore in a stationary population,

[TABLE]

which is a Poisson mixture with mean $-\theta\log X$, where $X$ has density $(n-1)(1-x)^{n-2}$, $0<x<1$. A calculation in Tavaré [35] (5.3.6) using the first line of (4) shows that

[TABLE]

There is the simple recursion that

[TABLE]

with $\mathbb{P}(S_{1}=s)=\delta_{s0}$.

The distribution of the number of alleles in a sample of $n$ in a stationary population is

[TABLE]

where $\{S_{k}^{n}\}$ are Stirling numbers of the first kind.

2.1 The joint distribution of $\widetilde{S}_{n}(t)$ and $A_{n}(t)$

Measuring time back from the present, $(\widetilde{S}_{n}(t),A_{n}(t))$ is a Markov process beginning at $(0,n)$ at $t=0$ such that

[TABLE]

The pgf of the joint distribution of $A_{n}(t),\widetilde{S}_{m}(t)$ distribution is derived in Griffiths [10], eqns (8), (9). The approach taken in the present paper is new.

The joint pgf of the number of mutations arising while there are $n,\ldots,k$ ancestors and the density of $T_{n},T_{n-1},\ldots,T_{k}$ in an $n$-coalescent tree is

[TABLE]

Integrating over $t_{n}+t_{n-1}+\cdots+t_{k}\leq t$, the pgf of the number of mutations on edges of the coalescent tree during times $T_{n},\ldots,T_{k}$ and the distribution function of $T_{n}+\cdots+T_{k}$ is

[TABLE]

The joint probability that $A_{n}(t)=l$ and pgf of the number of mutations $\widetilde{S}_{n}(t)$ arising in $(0,t)$, which will be denoted by $G_{l}(z;t)$, is the probability that $T_{n}+\cdots+T_{l+1}=\tau<t$ (a coalescence necessarily occurs at $\tau$), there is no coalescence in $(\tau,t)$ and the pgf of the numbers of mutations on $T_{n},\ldots,T_{l+1}$ and the $l$ lines from $\tau$ to $t$. Therefore, with notation $\theta_{z}=\theta(1-z)$,

[TABLE]

The identity (3) is used in the calculation. As a check of (LABEL:pgfprob:0), if $z=1$ the probability that there are $l$ ancestors at $t$ is $\mathbb{P}(A_{n}^{0}(t)=l)$, as it should be.

The marginal pgf of the number of mutations arising while there are at least $l$ ancestors is given by

[TABLE]

To see this, calculate the pgf of the total number of mutations accumulating while there are at least $l$ ancestors as

[TABLE]

In a calculation similar to (5),

[TABLE]

is the pgf of a Poisson mixture with a rate $-\theta\log X$, where $X$ has a Beta$(l-1,n-l+1)$ distribution. Hence the probability of $k$ mutations, from the pgf (11), is

[TABLE]

for $k=0,1,\ldots$.

The pgf of the number of mutations arising in $(0,t)$ (not counting mutations when there is one ancestor) is $\sum_{l=2}^{n}G_{l}(z;t).$ The joint pgf for $S_{n}(t)$, $\widetilde{S}_{n}(t)$ and the probability of $l$ ancestors time $t$ ago is $Q_{l}(r)G_{l}(z;t)$, where $Q_{l}(r)$ is the pgf of the number of segregating sites in a sample of $l$ from the population at the initial time. In a stationary population $Q_{l}(r)=H_{l}(r)$, defined in (4).

The probability $\mathbb{P}(\widetilde{S}_{n}(t)=k)$ does not have a simple form, but if one considers standing variation in a stationary population as well as mutations in $(0,t)$ then the distribution of the number of mutations is as in a stationary population at time $t$ and the conditional distribution of the number of ancestors given the number of mutations is easier. The pgf/probability of the number of mutations and number of ancestors at time $t$ in the past is then

[TABLE]

Summing (14) over $l=2,\ldots,n$, we see that the marginal pgf of the number of mutations is

[TABLE]

and the distribution of the number of ancestors at time $t$ back, found by setting $z=1$ in (14), is

[TABLE]

Equations (15) and (16) recover earlier results. The joint pgf and expected number of ancestors time $t$ ago is

[TABLE]

Inversion of (17) is straightforward, if a little messy. We calculate \mathbb{E}\big{[}A_{n}(t)\mid S_{n}=r\big{]}. Let $a(r,k)$ be the coefficient of $z^{r}$ in

[TABLE]

Then by equating coefficients of $z^{r}$ in

[TABLE]

for $k=1,2,\ldots,n$, $r=0,1,\ldots$, we get

[TABLE]

Note that $a(r,0)=\mathbb{P}(S_{n}=r)$. Let

[TABLE]

and finally

[TABLE]

Then

[TABLE]

The pgfs (LABEL:pgfprob:0) and (14) and the formula (18) are new.

In a later section we give a rejection algorithm for simulating the distribution of $A_{n}(t)$ given $S_{n}=r$.

2.2 The joint distribution of $\widetilde{S}_{n}(t)$, $\widetilde{K}_{n}(t)$, $A_{n}(t),A_{n}^{\theta}(t)$

A method for computing the stationary joint distribution of $(S_{n},K_{n})$ is derived in [11], exploiting a diffusion generator in a model with $K$ alleles, then letting $K\to\infty$. It is found numerically that the joint distribution is strongly diagonal in that, approximately, $S_{n}=K_{n}-1$. It is also shown that $S_{n}-K_{n}+1$ has a proper limit distribution when $n\to\infty$, even though both $S_{n}$ and $K_{n}$ tend to infinity. It is always true that $\widetilde{S}_{n}(t)-\widetilde{K}_{n}(t)+1\geq 0$.

Here we consider the joint distribution of $(\widetilde{S}_{n}(t),\widetilde{K}_{n}(t),A_{n}(t),A_{n}^{\theta}(t))$ using a coalescent treatment different from that of Griffiths [11]. To obtain a Markov process, consider $(\widetilde{S}_{n}(t),\widetilde{K}_{n}(t),A_{n}(t),A_{n}^{\theta}(t)),t\geq 0$, beginning from $(0,0,n,n)$. Let $B_{n}(t)=A_{n}(t)-A_{n}^{\theta}(t)$. Then

[TABLE]

The total coalescence rate is $a(a-1)/2=(a^{\theta}+b)(a^{\theta}+b-1)/2$. The total mutation rate is $a\theta/2=\theta(a^{\theta}+b)/2$. This process counts mutations and alleles as they arrive back in time from $t$ in a sample of $n$ and its ancestors, by considering the two groups of lines $a-a^{\theta},a^{\theta}$ and in which groups mutations or coalescences occur.

The simpler process $(A_{n}(t),A_{n}^{\theta}(t))$ has rates

[TABLE]

with total rate of $\left(a(a-1)+a^{\theta}\right)/2$. The marginal transition functions are known explicitly from (LABEL:trnQ) when $\theta=0$ and $\theta>0$.

A slightly different approach, described in equation (2.9) of Griffiths [10], is to consider sample paths that have $i$ mutations and $j$ alleles at $t$ starting from fixed $m=a^{\theta}$, $r=a$ as being stationary distributions which satisfy the following recursive system:

[TABLE]

for $i=0,1,\ldots,a^{\theta};j=1,2\ldots\text{~{}and~{}}a\geq a^{\theta}\geq 2.$ The boundary probabilities satisfy

[TABLE]

Then we want $p(i,j;n,n)$, which can be computed from (21). These equations can be argued from a probabilistic perspective, though a different approach is taken in Griffiths [10], where the identification $m=a^{\theta}$, $r=a$ is not made. (21) is analogous, but more complex, to a similar recursion for the number of mutations which could be written as

[TABLE]

for $i=0,1\ldots,\>a=2,3,\ldots$.

A pgf version of (21) for

[TABLE]

appears in [11], eq. (2.8).

3 Simulation-based approaches

In this section we develop theory that will be used in the rejection sampling and importance sampling schemes. We begin by recalling a simulation method that generates stationary samples of haplotype counts together with the number of mutations in the tree. This works for constant population size coalescent models.

3.1 Growing a tree

A useful way to simulate the ancestral history of haplotype configurations in age order together with the mutations is to use a condensation of an algorithm for growing a gene tree whose nodes are mutations. The algorithm for the tree, described in Ethier and Griffiths [6],Theorem 5.4 and Griffiths [12], p.7, is the following.

Algorithm 1

Start with a tree of two leaves (individuals) as two edges joined at a node. 2. 2.

When there are $m$ leaves, select one of the $m$ leaves at random and duplicate from the same immediate mutation node with probability $(m-1)/(\theta+m-1)$; or add a mutation node on the chosen edge with probability $\theta/(\theta+m-1)$. 3. 3.

To get a sample of $n$, stop when there are first $n+1$ leaves and select the configuration just before the last leaf appeared.

Times between events can be added as exponential random variables $T^{\theta}_{2},\ldots,T^{\theta}_{n}$.

If we just look at haplotype frequencies in age order from the oldest and keep track of $s$, the number of accumulated mutations, the state space is $(m_{1},\ldots,m_{k};s)$ and transitions are Markovian. The condensed algorithm is the following.

Algorithm 2

Start with a configuration $(n_{1}=2;0)$ of two identical oldest haplotypes and no mutations. 2. 2.

When there are $m$ individuals make a transition $(m_{1},\ldots,m_{k};s)\to(m_{1},\ldots,m_{j}+1,\ldots,m_{k};s)$ with probability $(m_{j}/m)(m-1)/(m+\theta-1)$; or $(m_{1},\ldots,m_{k};s)\to(m_{1},\ldots,m_{l}-1,\ldots,m_{k},1;s+1)$ with probability $(m_{l}/m)\theta/(m+\theta-1)$. 3. 3.

To get a sample of $n$, stop when there are first $n+1$ leaves and select the configuration just before the last leaf appeared.

This algorithm is useful for simulation of an ancestral path forward in time which contains full information of haplotype count configurations and mutations.

Let $p^{\circ}(\bm{m})/\prod_{j=1}^{m}\alpha_{j}!$ be the probability of a non-age labelled configuration $\bm{m}$ (labelled in an arbitrary order), with $\alpha_{j}$ the number of allele frequencies equal to $j$, under the Markov chain without the stopping rule. Then

[TABLE]

The sample probability $p(\bm{n};s)/\prod_{j=1}^{n}\alpha_{j}!$ is such that

[TABLE]

A recursion is therefore

[TABLE]

The recursion in (25) is used in calculating the importance sampling weights in Section 4.

Note that the last sum includes the case when $i=l$, and it can be written as

[TABLE]

If $a,s$ are the number of haplotypes and mutations, respectively, with $\bm{n}$ then $s-a+1\geq 0$ so $p(\bm{n};s)=0$ if $s-a+1<0$ in the recursive equations. The Ewens Sampling formula [7] satisfies a similar equation to (25), summing over $s$, with

[TABLE]

The probability of the sample configuration is

[TABLE]

The joint distribution of $(\bm{n};s)$ is studied in Innan et al [26], where their derived equation (1) is essentially the same as (25). The authors then use these recursive equations to derive exact probabilities for small sample sizes. In Section 4 we use an importance sampling technique based on (25) which allows calculation for much larger sample sizes than are possible with exact calculation and we also carry out ancestral inference back in time.

3.2 The number of ancestors and mutations

The simplest ancestral question is to ask about the distribution of $A_{n}(t)$ conditional on $S_{n}=s$ segregating sites in the sample of size $n$. Rejection algorithms work well for problems like this, as was illustrated by Tavaré, Balding, Griffiths, Donnelly [36]. They took a general Bayesian approach in which $\theta$ was considered as a random variable and the times $T_{n},\ldots,T_{2}$ come from a variable population size coalescent model. The combinatorics for a general binary coalescent tree are studied in [22], where the sample frequency spectrum and the mean age of a mutation known to have $b$ descendents in a sample of $n$ genes are calculated.

We illustrate by describing how to generate observations from the conditional distribution of $(\theta,A_{n}(t))$ conditional on $S_{n}=s$. To this end, let $W_{1}=T_{n},W_{2}=T_{n}+T_{n-1},\ldots,W_{j}=T_{n}+\cdots+T_{n-j+1},\ldots,W_{n-1}=T_{n}+\cdot+T_{2}$. $W_{n-1}$ is the height of the coalescent tree. Define the total length of the tree as $L_{n}=nT_{n}+\cdots+2T_{2}$. Form a set of bins as follows:

[TABLE]

Let $J$ be the bin that covers $t$. Then the number of ancestors at time $t$ is

[TABLE]

and the length of the coalescent tree from 0 back to $t$ is

[TABLE]

Finally, define the length of the coalescent tree from time $t$ to the most recent common ancestor as $L_{n}(t)=L_{n}-\widetilde{L}_{n}(t)$. Let $\widetilde{S}_{n}(t)$ be the number of mutations arising in $(0,t)$ and let $S_{n}(t)$ the number of mutations from $t$ to the most recent common ancestor. Then, conditional on $T_{n},\ldots,T_{2}$ and $J$,

[TABLE]

where ${\rm Po}(\lambda)$ denotes the Poisson distribution with parameter $\lambda$, such that

[TABLE]

and $\widetilde{S}_{n}(t)$ and $S_{n}(t)$ are conditionally independent. The total number of segregating sites in the sample at time 0 is $S_{n}=\widetilde{S}_{n}(t)+S_{n}(t)$. The simplest algorithm gives the probability distribution of $J$, and therefore the distribution of $A_{n}(t)$, conditional on $S_{n}=s$.

Algorithm 3

Simulate $\theta$ from the prior, $\pi(\cdot)$ 2. 2.

Simulate $T_{n},\ldots,T_{2}$ from an appropriate coalescent model 3. 3.

Compute $J$, $A=A_{n}(t)=n-J+1$, $L_{n}$ 4. 4.

Accept $(\theta,A)$ as an observation from the posterior with probability

[TABLE]

We also note that the same rejection approach may be used to approximate conditional distributions for many other ancestral variables. For example, to study the distribution of the number of mutations $S_{n}(t)$ present in the ancestors at time $t$, which is a measure of the standing variation at time $t$, we can use the following.

Algorithm 4

Simulate $\theta$ from the prior, $\pi(\cdot)$. 2. 2.

Simulate $T_{n},\ldots,T_{2}$ from an appropriate coalescent model. 3. 3.

Compute $J$, $A=A_{n}(t)=n-J+1$, $L_{n}(t)$ and $\widetilde{L}_{n}(t)$, defined in (27). 4. 4.

Accept $(\theta,A)$ with probability

[TABLE]

and else return to 1. 5. 5.

Simulate $S$ from a Binomial($s,L_{n}(t)/L_{n}$) distribution, and return $(\theta,A,S)$ as an observation from the posterior of $(\theta,A_{n}(t),S_{n}(t))$ given $S_{n}=s$.

We may treat $\theta$ as fixed in this approach (that is, as having a degenerate prior), the approach then addressing the problems studied in the first section of the paper.

Blum and Rosenberg [4] construct a rejection algorithm for maximum likelihood estimation of the number of ancestral lineages at time $t$ back based on the sample frequency spectrum. The algorithm is constructed differently from Algorithms 3 and 4.

3.3 Hammer et al. example

As an illustration we consider the Y chromosome data of 1544 sequences from Hammer et al. [24]. In this paper a perfect phylogeny was constructed from the sequence data and the program GENETREE was used to find the TMRCA and ages of mutations in the ancestral tree shown in Figure 7 of [24] with $\theta=2.5$. There were 9 segregating sites and 10 haplotypes observed in the data. The unconditional expected height of the coalescent tree is 2 time units, and the expected height conditional on $s=9$ segregating sites is 1.21 units, this latter found from the method in [36].

We generated 10,000 repetitions of the previous algorithm for a series of times $t$, obtaining the information in Table 1 for the conditional expectations of $S_{n}(t)$ and $A_{n}(t)$.

The simple rejection schemes illustrated here are not as useful for considering more complicated summaries of the data. In the next section we show how to exploit an importance sampling approach to derive conditional distributions given the haplotype frequency distribution and the number of mutations.

4 Importance sampling

Sequential importance sampling for ancestral inference in population genetics has a long history, illustrated by Griffiths and Tavaré [18, 19, 23], Felsenstein et al. [8], Stephens and Donnelly [32], Griffiths [14], De Iorio and Griffiths [5]. The technique can be described as constructing a proposal distribution for events back in time, simulating back in time, then correcting for the approximate proposal distribution by calculating the exact probability of the path forward in time and taking the ratio of the probability of the forward path divided by the approximate probability of the backward path as the importance weight. If there are $r$ simulation runs then an empirical ancestral history is returned as $(\widehat{p}_{1},{\cal H}_{1}),\ldots,(\widehat{p}_{r},{\cal H}_{r})$ where $\{\widehat{p}_{j}\}$ are the importance weights scaled to add to unity and $\{{\cal H}_{j}\}$ are the histories. A general reference to sequential importance sampling is Liu [29].

Choosing a proposal distribution is an art. We use the principal of choosing a lineage which can be involved in a transition back in time uniformly. This has a theoretical justification, described in [32, 5]. Sequential importance sampling for a haplotype configuration can be regarded as a simplification of the technique used for a complete genetree, constructed as a perfect phylogeny from the pattern of mutations on DNA sequences. The simpler scheme counts different haplotypes of the sequences, the extra information being the number of mutations back to the most recent common ancestor. Time information, such as coalescence times, ages of mutations and time to the most recent ancestor can be included. Details of how to include time in an importance sampling algorithm are in Griffiths and Tavaré [21]. We develop a new importance sampling approach for the Kingman coalescent models conditional on an observed configuration $\bm{n},s$.

The proposal distribution $\widehat{p}$ for reverse transitions in a haplotype history is detailed in the following equations. Suppose a current configuration is $\bm{n}=(n_{1},\ldots,n_{k})$, the number of mutations to the most recent common ancestor is $s$, and the number of singletons is $q$.

For $k>2$, if $s-k+1>0$,

[TABLE]

or if $s-k+1=0$ then

[TABLE]

The first factors in (29) involve a choice of either mutations that define allele types and those which appear on lineages between defined alleles. Importance weights for transitions back in time are therefore

[TABLE]

When $k=2$ we have to consider the following possible cases:

(a) If $n_{1}>1,n_{2}>1$

[TABLE]

(b) if $n_{1}>1,n_{2}=1,s>1$,

[TABLE]

(c) if $n_{1}>1,n_{2}=1,s=1$,

[TABLE]

(d) similarly when $n_{1}=1,n_{2}>1$,

(e) if $n_{1}=1,n_{2}=1,s>1$,

[TABLE]

(f) if $n_{1}=1,n_{2}=1,s=1$

[TABLE]

Importance weights are:

[TABLE]

4.1 Implementation

Our implementation provides

The probability of a sample configuration of haplotypes and number of segregating sites. This is an extension of the Ewens Sampling Formula, which is the probability of the configuration of haplotypes.

The next calculations are conditional on the configuration of haplotypes and segregating sites at time 0.

The average coalescence times in the past.

• 2.

The average mutation times in the past.

• 3.

The average times when alleles are lost in the past. (The time of loss of the last haplotype is truncated at the TMRCA if not lost by mutation.)

• 4.

The average allele ages in the past.

• 5.

The distribution of ancestor lines and the average allele configuration at a given time in the past.

The program also implements a variable population size option with exponential growth. Coalescent times then have a distribution that depends on the time when they occur. We do not go into detail here, but refer the reader to Griffiths and Tavaré [20]. The analogue of the Ewens Sampling Formula in this case is derived in Griffiths and Lessard [16].

Accuracy of the implementation was checked by ensuring that for smaller sample sizes the equation

[TABLE]

was approximately satisfied, where $k$ is the number of alleles in $\bm{n}$ and $B$ is a suitable upper bound. The simulation variance, starting with different seeds, was observed to be small for sample sizes such as in example 4.2 that follows.

4.2 Hammer et al. example, continued

We continue with the example started in Section 3.3. The Y haplotype data of $n=1,544$ sequences from Hammer et al. [24], had 10 haplotypes and $s=9$ segregating sites. We continue to use their value of $\theta=2.5$ for illustration. The 10 haplotype frequencies are

[TABLE]

in the lineage order shown in Figure 7 of [24]. The average values in the tables below are conditional on the configuration and number of segregating sites, thereby extending the results of Section 3.3.

The Appendix describes the input for the implementation of the method. A command line of

esf_stl HammerHap.dat 10 9 2.5 1000000 93849 -a

in which the input file HammerHap.dat contains the haplotype frequencies in the order above, produces the output described below; the average coalescence times are not shown.

Two runs with different seeds gave identical output to three significant places, showing some confidence in the output. The probability of obtaining the sample configuration and $s=9$ segregating sites was $1.4785\times 10^{-19}$. As a comparison the probability of the sample configuration, calculated from the Ewens Sampling Formula, was $1.1722\times 10^{-18}$. The mean TMRCA, conditional on the data, in coalescent units was $1.15$, which may be compared to the value of $1.21$ obtained in Section 3.3.

4.2.1 Stationary properties

Here we record some information about the sample at time [math].

The conditional expected mutation times in increasing time order are

[TABLE]

while the conditional expected haplotype loss times in increasing time order are

[TABLE]

Most of the tree structure has developed by an average time of less than 1.00 coalescent time unit.

The conditional expected haplotype ages, in the order they are listed above, are

[TABLE]

These are monotonic in the number of copies of the haplotype in the sample, confirming the intuition that common haplotypes tend to be older.

4.2.2 Time-varying properties

At a given time $t$ in the past, the distribution of the configuration, number of mutations, and number of ancestral lineages conditional on the current configuration and number of mutations can be calculated by the importance sampling program. We illustrate this by considering time points $t=0.1,0.5,1.0,1.5$ and taking averages at those times.

We begin by comparing the conditional distribution of $A_{n}(t)$ and $S_{n}(t)$ with the analogous results in Table 1. Additionally, Table 2 shows the average of the number of haplotypes, $K_{n}(t)$, present at time $t$.

The results in the second and third columns should be compared with those in Table 1; they show qualitatively the same results. Table 3 shows the relative errors for estimates of $S_{n}(t)$ and $A_{n}(t)$. For example, letting superscripts 1 and 2 denote estimates from Tables 1 and 2, the relative error for $S_{n}(t)$ is $|S_{n}^{1}(t)-S_{n}^{2}(t)|/\frac{1}{2}(S_{n}^{1}(t)+S_{n}^{2}(t))$.

Table 4 gives the number of haplotype counts at different times $t$ in the past. Haplotypes decrease because of coalescence and types are eventually lost when their defining mutation takes place. In Table 5 we give the distribution of the number of ancestral lines at different times $t$ in the past, conditional on the current haplotype configuration and the number of segregating sites. The most interesting time configurations are when $t\leq 0.1$; afterwards the number of lineages and number of haplotypes decrease rapidly.

4.3 A 1000 Genomes Y chromosome dataset

A larger Y chromosome data set comes from the 1000 Genomes Project. An analysis of these data is made in Poznik et al. [30], where a phylogeny is constructed. The paper concludes that the data show evidence of expansion. As an example of our approach, we consider a subset of this data set consisting of sequences in the A, B and E haplotype groups. These are the three oldest groups in the phylogeny, and are composed of 334 sequences. There may be explanations other than expansion for the data configuration, such as a non-random choice of individuals, however we will assume a random sample for this example analysis. For illustration we focus on the TBL1Y gene, composed of some 180,000bp, and containing 278 biallelic SNPs. The haplotype configuration, with $\alpha_{j}$ equal to the number of alleles of multiplicity $j$, is given in Table 6.

Watterson’s estimate of $\theta$ [38] based on the number of segregating sites $s$ in a constant-sized population is

[TABLE]

The maximum likelihood estimate $\widehat{\theta}_{E}$ based on the Ewens’ sampling formula uses $k=\sum_{j=1}^{n}\alpha_{j}$, a sufficient statistic for $\theta$. $\widehat{\theta}_{E}$ satisfies

[TABLE]

In the TBL1Y dataset $\widehat{\theta}_{E}=82$. The large number of singletons $\alpha_{1}$ in the data suggests exponential growth in the population. Growth produces a star shaped coalescent tree, which leads to a greater number of singleton sequences. The mean number of singletons in the constant size population setting is

[TABLE]

If $\theta=82$, then $\mathbb{E}\left[\alpha_{1}\right]=66$, which is much less than the observed $\alpha_{1}=107$. Tajima’s $D$ [33] is given by

[TABLE]

where $\pi$ is the average number of pairwise differences, an unbiased estimate of $\theta$. This may be used to test for population growth or other departures from the coalescent model with no growth. Large negative values of $D$ indicate population growth. In our data $\pi=6.49$ and $D=-2.6$, consistent with expansion. In the Appendix we describe another statistic for testing the no-growth model based on the frequency spectrum, and particularly on the number of singletons, when $\theta$ is large. Achaz [1] develops neutrality tests based on the frequency spectrum which generalize tests based on the number of segregating sites and Tajima’s $D$.

In this large data set it is difficult to obtain a very precise estimate of $\theta$ and growth rate $\beta$ because there is a large amount of variation in the importance sampling scheme due to the number and length of the sequences. This is not so much an importance sampling issue, but due to the size and structure of the data. Random subsets of the data could be chosen, but the main feature of the data is the large number of singletons. Ancestral inference of the number of ancestral lines at $t$ back would also be difficult to interpret for subsamples. We try a large value $\theta=100$ with different growth rates $\beta$. Growth decreases the variation in the sample, but increases the proportion of singletons, because the coalescent lengths are shortened and the tree is star shaped. Increasing $\theta$ with growth keeps the variation as well as increasing the number of singletons. The likelihood of the allele configuration and number of segregating sites was calculated for $\theta=100$ with several values of the growth rates $\beta$. Two different replicates each with 10 million runs gave the results in Table 7.

A plausible maximum likelihood estimate when $\theta=100$ is $\widehat{\beta}=1.0$. If $\theta$ and $\beta$ are increased together, it is possible that the likelihood estimates fall on a ridge. In the first replicate with these values of $\theta$ and $\beta$ the TMRCA was $1.461$ and the average of ages within halpotype groups are shown in Table 8.

With the large value of $\theta=100$, times where mutation creates an allele are close to the leaves of the tree.

The scale of modern molecular datasets points out the difficulty of exact inference techniques, and highlights the need for alternative approaches. Among these are the Approximate Bayesian Computation (ABC) and other advanced simulation methods, to which Paul Joyce made several contributions. Two examples of Paul’s research are in estimating evolutionary rates of trait evolution by ABC [31], and a perfect simulation method for simulation from a non-neutral high dimensional distribution of allele frequencies [27].

5 Acknowledgements

We thank Chris Tyler-Smith, Yali Xue and two referees for their helpful comments about the 1000 Genomes Y chromosome data.

Code for the rejection and importance sampling methods may be obtained from the authors.

Appendix

Importance sampling code

The input for the importance sampling method illustrated in Section 4.2 is:

esf_stl configfile k [#alleles] m [#mutations] theta replicates seed

ΨOptions Ψ-g beta [exponential growth] Ψ-a [age information] Ψ-t time [Configuration at time] Ψ ΨBob GriffithsΨ4 May 2017, Version 1.7

Poisson approximation for large $\theta$

Motivated by the discussion in Section 4.3 we discuss the behaviour of the Ewens sampling Formula for large values of $\theta$ and $n$. The Ewens Sampling Formula gives the distribution of the number of haplotypes and their frequencies in a sample taken from a constant-size population. Writing $\alpha_{j}$ for the number of haplotypes with frequency $j$, the distribution is

[TABLE]

where $\theta_{(n)}:=\theta(\theta+1)\cdots(\theta+n-1)$, and $\alpha_{1}+2\alpha_{2}+\cdots+n\alpha_{n}=n$. The formula (31) shows that, were it not for the condition that $\alpha_{1}+2\alpha_{2}+\cdots+n\alpha_{n}=n$, the $\alpha_{j}$ would be independent Poisson random variables with mean $\theta/j$. Indeed, for fixed $\theta$ it is known that for any $b=o(n)$ as $n\to\infty$, the total variation distance between the distribution of $(\alpha_{1},\ldots,\alpha_{b})$ and that of $(Z_{1},\ldots,Z_{b})$, for independent Poisson random variables with $\mathbb{E}\left[Z_{j}\right]=\theta/j$, is $O(b/n)$ as $n\to\infty$. See Arratia, Barbour and Tavaré [2], Theorem 5.2.

Here we consider the case in which $\theta\to\infty$ with $n$, and we show that for fixed $b$, $(\alpha_{1},\ldots,\alpha_{b})$ has asymptotically the distribution of $(Z_{1}^{\theta},\ldots,Z_{b}^{\theta})$, where

[TABLE]

To see this, consider the joint falling factorial moments of $\alpha_{1},\ldots,\alpha_{b}$, given by Watterson [37] as

[TABLE]

where $m:=r_{1}+2r_{2}+\cdots+br_{b}$. The term on the right gives the falling factorial moments of $(Z_{1}^{\theta},\ldots,Z_{b}^{\theta})$, and the result follows from the method of moments.

In practice, different limit laws are obtained depending on the way $\theta$ varies with $n$. For example, if $\theta\sim\eta n$, then $\alpha_{1}$ has approximately a Poisson distribution with mean $\theta/(1+\eta)$. For the data in Section 4.3, with $\theta=82$, the number of singleton haplotypes has mean $82\cdot(334/416)\approx 65.84$. Since the probability of observing 107 or more singletons is then $\approx 1.92\times 10^{-6}$, we conclude that the constant-size model does not provide an adequate fit. In a similar spirit, $\alpha_{1}+\alpha_{2}$ has approximately a Poisson distribution with mean 92.27. We observed $\alpha_{1}+\alpha_{2}=119,$ the probability of a larger value being $\approx 0.0043$; once more, this suggests the constant-size model is not a good fit.

6 References