The transition distribution of a sample from a Wright-Fisher diffusion with general small mutation rates

TL;DR

This paper derives an approximation for the transition distribution of a sample from a Wright-Fisher diffusion with small mutation rates using a coalescent approach, addressing a previously unknown transition density.

Contribution

It provides the first known approximation of the transition density for Wright-Fisher diffusion with general small mutation rates.

Findings

Derived a coalescent-based approximation for the transition distribution.

Established an explicit form for the transition density in the small mutation regime.

Enhanced understanding of genetic variation dynamics in population genetics.

Abstract

The transition distribution of a sample taken from a Wright-Fisher diffusion with general small mutation rates is found using a coalescent approach. The approximation is equivalent to having at most one mutation in the coalescent tree of the sample up to the most recent common ancestor with additional mutations occurring on the lineage from the most recent common ancestor to the time origin if complete coalescence occurs before the origin. The sampling distribution leads to an approximation for the transition density in the diffusion with small mutation rates. This new solution has interest because the transition density in a Wright-Fisher diffusion with general mutation rates is not known.

Equations

$$\mathcal{L}=\frac{1}{2}\sum _{i,j=1}^d{x}_i({\delta }_{ij}-x_j)\frac{{\partial }^2}{\partial x_i\partial x_j}+\sum _{i,j=1}^d{\gamma }_{ji}{x}_j\frac{\partial }{\partial x_i},$$

$$m(\mathbf{n},n,\mathbf{p})=\left(\genfrac{}{}{0ex}{}{n}{\mathbf{n}}\right)\prod _{j=1}^d{p}_j^{n_j},$$

$$p(\bm{n};t,\bm{x})=\mathbb{E}\big{[}m(\bm{n};n,\bm{X}(t))\mid\bm{X}(0)=\bm{x}\big{]}$$

$$\mathcal{L}{x}_1^{n_1}\cdots x_d^{n_d}=$$

$$\frac{1}{2}\sum _{i=1}^d{n}_i(n_i-1)x_1^{n_1}\cdots x_i^{n_i-1}\cdots x_d^{n_d}-\frac{1}{2}\sum _{i,j=1}^d{n}_i(n_j-{\delta }_{ij})x_1^{n_1}\cdots x_d^{n_d}$$

$$+\frac{\theta }{2}\sum _{i,j=1}^d{P}_{ji}{n}_i{x}_1^{n_1}\cdots x_i^{n_i-1+{\delta }_{ij}}\cdots x_j^{n_j+1-{\delta }_{ij}}\cdots x_d^{n_d}-\frac{\theta }{2}\sum _{i=1}^d{n}_i{x}_1^{n_1}\cdots x_d^{n_d}.$$

$$\frac{d}{dt}p(\mathbf{n};t,\mathbf{x})$$

$$p(n{\mathbf{e}}_a;\theta )$$

$$p(n_a{\mathbf{e}}_a+n_b{\mathbf{e}}_b;\theta )$$

$$p(\mathbf{n};\theta )$$

$$q_{nl}(t):=P(A_n(t)=l)=\sum _{j=l}^n{\rho }_j(t)(-1{)}^{j-l}\frac{(2j-1)(j+l-2)!}{l!(l-1)!(j-l)!}\frac{n_{[j]}}{n_{(j)}},$$

$$\begin{array}{r}{n}_{(j)}=n(n+1)\cdots (n+j-1),\\ n_{[j]}=n(n-1)\cdots (n-j+1),\end{array}$$

$$f_{nl}(t)dt$$

$$q_{nl}(t)={\left(\genfrac{}{}{0ex}{}{l}{2}\right)}^{-1}{f}_{nl}(t).$$

$$r_{\mathbf{n}\mid \mathbf{l}}$$

$${\mathcal{D}\mathcal{M}}_{\mathbf{\theta },n}(\mathbf{n})$$

$${\mathcal{D}}_{\mathbf{\theta }}(\mathbf{x})=\frac{\Gamma (\theta )}{\Gamma ({\theta }_1)\cdots \Gamma ({\theta }_d)}{x}_1^{{\theta }_1-1}\cdots x_d^{{\theta }_d-1},$$

$$r_{\mathbf{n}\mid \mathbf{l}}=\left(\genfrac{}{}{0ex}{}{n-l}{\mathbf{n}-\mathbf{l}}\right)\cdot \frac{(l-1)!}{{\prod }_{l_a>0}(l_a-1)!}\cdot \frac{{\prod }_{l_a>0}(n_a-1)!}{(n-1)!}.$$

$$\mathbb{P}(n_i,n_j\mid \mathbf{l})$$

$$p_{nk}(c)$$

$$\displaystyle x_{a}^{2}e^{-\theta t}e^{-t}+\big{(}2x_{a}^{2}P_{aa}+\sum_{j\neq a}2x_{a}x_{j}P_{ja}\big{)}\frac{\theta}{2}te^{-\theta t}e^{-t}+o(\theta)$$

$$\displaystyle~{}~{}=x_{a}^{2}(1-\theta t)e^{-t}+\big{(}\sum_{j\in[d]}2x_{a}x_{j}P_{ja}\big{)}\frac{\theta}{2}te^{-t}+o(\theta).$$

$$\sum _{j\in [d]}{x}_j{\int }_0^t{\left(e^{(t-w)\gamma }\right)}_{ja}{e}^{-\theta w}(1+P_{aa}\theta w)e^{-w}dw+o(\theta )$$

$$p(2{\mathbf{e}}_a;\theta )={\pi }_a(1-\theta (1-P_{aa}))+o(\theta ).$$

$$p(2{\mathbf{e}}_a;t,\mathbf{x})=x_a^2(1-\theta t)e^{-t}+\left(\sum _{j\in [d]}2x_a{x}_j{P}_{ja}\right)\frac{\theta }{2}t{e}^{-t}+x_a(1-e^{-t})+$$

$$2x_a{x}_b{e}^{-\theta t}{e}^{-t}+x_a^2\frac{\theta }{2}t{e}^{-\theta t}2P_{ab}{e}^{-t}+\sum _{j\ne a}2x_a{x}_j\frac{\theta }{2}t{P}_{jb}{e}^{-\theta t}{e}^{-t}$$

$$+x_b^2\frac{\theta }{2}t{e}^{-\theta t}2P_{ba}{e}^{-t}+\sum _{j\ne b}2x_b{x}_j\frac{\theta }{2}t{P}_{ja}{e}^{-\theta t}{e}^{-t}+o(\theta )$$

$$=2x_a{x}_b(1-\theta t)e^{-t}+\theta t{e}^{-t}\sum _{j\in [d]}\left(x_a{x}_j{P}_{jb}+x_b{x}_j{P}_{ja}\right)+o(\theta ).$$

$$\sum _{j\in [d]}{x}_j{\int }_0^t{\left(e^{(t-w)\gamma }\right)}_{ja}{e}^{-\theta w}{P}_{ab}\theta w{e}^{-w}dw+a\leftrightarrow b+o(\theta )$$

$$p({\mathbf{e}}_a+{\mathbf{e}}_b;\theta )=\theta ({\pi }_a{P}_{ab}+{\pi }_b{P}_{ba})+o(\theta ).$$

$$p({\mathbf{e}}_a+{\mathbf{e}}_b;t,\mathbf{x})$$

$$\gamma=\left(\begin{array}[]{rrr}-0.00056&0.00032&0.00024\\ 0.00080&-0.00360&0.00280\\ 0.00100&0.00100&-0.00200\end{array}\right),$$

$$q_{nl}(k;t)$$

$$\omega_{jk}(t)=\begin{cases}te^{-{k\choose 2}t}&\text{~{}if~{}}j=k\\ \frac{1}{{j\choose 2}-{k\choose 2}}\Big{(}e^{-{k\choose 2}t}-e^{-{j\choose 2}t}\Big{)}&\text{~{}if~{}}j\neq k.\end{cases}$$

Full text

††journal:

1 Introduction

Burden and Tang (2016, 2017) find an approximation for the stationary distribution in a $d$-allele neutral Wright-Fisher diffusion with general small mutation rates. The model has a connection to boundary processes in which only transitions involving non-segregating and bi-allelic states are considered. Schrempf and Hobolth (2017),Vogl and Bergman (2015); De Maio, Schrempf, Kosiol (2015) study this boundary mutation property.

The Kingman coalescent process is dual to the Wright-Fisher diffusion in describing the ancestral history of a sample of $n$ individuals in a population back in time. In Burden and Griffiths (2019) a coalescent approach with mutations in the tree is used to find an approximate sampling formula for small mutation rates which then leads to an approximation for the stationary distribution in the diffusion. They show that finding approximate formulae for small rates is equivalent to considering at most one mutation in a coalescent tree.

In the current paper we find the sampling distribution in the Wright-Fisher diffusion with general mutation rates at time $t\geq 0$ with small rates by using a coalescent approach. The probability distribution of a gene tree at time $t$ is found in a similar way. A gene tree is equivalent to a sample of sequences in the infinitely-many-sites model.

The details of the sampling distribution in the general mutation model are more complex than in a stationary model because the coalescent tree back in time is constrained by the origin. The sampling distribution leads to an approximation for the transition density in the general mutation model which is of interest as a solution to the diffusion process, in population genetics, and more widely as an approximate differential equation solution. Vogl (2014), Theorem 3, considers a spectral expansion of the transition function in a two-allele model when mutation rates are small. The eigenfunctions are known in a two-allele model, but are unknown if there are more than two alleles, so this approach is not used here.

Burden and Griffiths (2018) also find an approximation to the stationary density in a two island, two allele model when mutation and migration rates are small by using a flux argument in the Wright-Fisher diffusion process as well as a coalescent argument.

2 Preliminaries

A $d$-dimensional Wright-Fisher diffusion process $\{\bm{X}(t)\}_{t\geq 0}$ with general mutation rates has a generator

[TABLE]

where $\gamma=(\gamma_{ij})=\frac{\theta}{2}(P-I)$, where $\theta/2$ is the overall mutation rate away from each type and $P$ a $d\times d$ transition probability matrix for type changes which has non-negative entries and row sums $1$. There are different possible pairs $(\theta,P)$ corresponding to a fixed rate matrix $\gamma$ with non-negative off-diagonal entries and row sums [math]. Note that for a choice $(\theta,P)$ with $\gamma$ fixed $P=\frac{2}{\theta}\gamma+I$, so $\theta$ has to be chosen so the off-diagonal entries of $P$ are non-negative. That is, choose $\theta$ large enough so that $\frac{2}{\theta}\max_{1\leq i\leq d}(-\gamma_{ii})\leq 1$. The minimal value of $\theta$ in choosing a pair $(\theta,P)$ for a given $\gamma$ is $\theta=2\max_{1\leq i\leq d}(-\gamma_{ii})$. This may define an appropriate pair $(\theta,P)$ for us to use in this paper because we are considering small mutation rates to order $\theta$. It is possible to use larger values of $\theta$, for example $\theta=\sum_{i\neq j}\gamma_{ij}=\sum_{i=1}^{d}(-\gamma_{ii})$ which is like a matrix norm. There may be other biological reasons to consider a particular $(\theta,P)$ pair and regard $\gamma$ as being determined by that pair. Explicit analytic solutions for the stationary distribution, Wright (1969), and the transition density in the diffusion, Griffiths (1979); Ethier and Griffiths (1993); Griffiths and Spanó (2010), are only known for parent independent mutation when $P$ has identical rows, in which case the stationary density is a Dirichlet distribution.

The coalescent duality with the Wright-Fisher diffusion is well known in population genetics. Briefly the argument is the following. Denote the multinomial distribution by

[TABLE]

where $\bm{n}=(n_{1},\ldots,n_{d})$, $\sum_{j=1}^{d}n_{j}=n$. Let

[TABLE]

be the probability distribution of a configuration $\bm{n}$ in a sample of $n$ individuals taken at time $t$ when the initial frequencies in the diffusion are $\bm{x}$. Note that

[TABLE]

Simplifying Eq. (3) and using basic theory

[TABLE]

Arguing probabilistically Eqs. (4) for $p(\bm{n};t,\bm{x})$ are the same as those for the sampling distribution in the coalescent with general mutations governed by $\theta$ and $P$ on the edges. The overall coalescence rate when $n$ individuals is $\frac{1}{2}n(n-1)$ and the overall mutation rate is $\frac{1}{2}n\theta$. If the first event back in time was a coalescence event then from the $n-1$ individuals after coalescence the configuration must be $\bm{n}-\bm{e}_{i}$ and a type $i$ individual chosen to give birth with probability $(n_{i}-1)/(n-1)$ to obtain a configuration of $\bm{n}$ for $i\in[d]$. If the first event back in time was a mutation, then from the $n$ individuals after the event to obtain a configuration of $\bm{n}$ a type $j$ individual mutates to a type $i$ individual from a configuration $\bm{n}-\bm{e}_{i}+\bm{e}_{j}$ with probability $(n_{j}+1-\delta_{ij})/n$ for $i,j\in[d]$.

Burden and Griffiths (2019) find that for small mutation rates, $\theta\rightarrow 0$, the probability of a configuration of $\bm{n}$ in a sample of $n$ individuals taken in the stationary distribution of the diffusion with generator in Eq. (1), $p(\bm{n};\theta)=\lim_{t\rightarrow\infty}p(\bm{n};t,\bm{x})$, is

[TABLE]

where $\bm{\pi}$ is the stationary distribution of $P$.

The main result of this paper, Theorem 4.1 in Section 4, extends the sampling formula to finite times $0\leq t<\infty$. Suppose a sample of genes taken at time $t$ is composed of families of genes from $l\geq 1$ founder lineages at time 0. Fig. 1 illustrates two possibilities contributing to the probability $p(\bm{n};t,\bm{x})$, which need to be considered. If there are $l>1$ founder lineages, then to $\mathcal{O}(\theta)$ there is at most one mutation event in the coalescent tree, and hence at most one mutant family arising from a mutation. However, if there is $l=1$ founder lineage, then it turns out that to $\mathcal{O}(\theta)$ there is at most one mutation event in that part of the coalescent tree between the most recent common ancestor of the sample (MRCA) and time $t$, but any number of coalescent events along the single founder edge from time zero to the MRCA. The reason for this is related to the asymptotic transition to the stationary distribution Eq. (5) as $t\rightarrow\infty$, as will become apparent in the next section. The proof of Theorem 4.1 uses this description and works through the combinatorial details of the type and size of families in a sample. Unlike in the stationary distribution there can be more than two types present when there are greater than two families of founder lineages.

The number of edges in a coalescent tree of $n$ individuals is a pure death process beginning at $n$ with death rates $\mu_{k}={k\choose 2}$, $k=n,\ldots,2$. The times between events are exponential random variables $T_{k}$ with rates $\mu_{k}$, $k=n,\ldots,2$. The number of edges in a coalescent tree at time $t$ back from the present time, $A_{n}(t)$, has a probability distribution

[TABLE]

where the notation is that

[TABLE]

and $\rho_{j}(t)=\exp\{-{j\choose 2}t\}$, (Tavaré, 1984; Griffiths, 1980). Denote $q_{nl}(t)=\mathbb{P}(A_{n}(t)=l)$ and the probability density of $T_{n}+\cdots+T_{l}$ as $f_{nl}(t)$. Then

[TABLE]

from which follows the identity

[TABLE]

As $n\rightarrow\infty$, $q_{\infty l}(t)$ and $f_{\infty l}(t)$ are proper distributions where $n_{[j]}/n_{(j)}$ is replaced by 1 in Eq. (6).

Combinatorics of the edge configuration in a coalescent tree are related to a Polya urn model. The probability $r_{\bm{n}\mid\bm{l}}$ of a configuration of edges $\bm{l}$ of $d$ types giving rise to a configuration $\bm{n}$ ($n>l$) is equivalent to considering a Polya urn beginning with $\bm{l}>\bm{0}$ coloured balls in which $\bm{n}-\bm{l}$ draws of the $d$ balls are made. If a ball of colour $j\in[d]$ is chosen on a draw from the urn, then it is replaced together with an additional ball of colour $j$. From classical theory, for which see Griffiths and Tavaré (2003),

[TABLE]

Here ${\cal DM}$ is the Dirichlet-multinomial distribution,

[TABLE]

defined in terms of the Dirichlet distribution,

[TABLE]

where $\bm{x}\in\Delta_{d}=\{0<x_{i}<1,i\in[d],|\bm{x}|=1\}$. That is, ${\cal DM}_{\bm{\theta},n}(\bm{n})$ is probability of drawing a configuration $\bm{n}$ in a sample of $n$ individuals apportioned with random probabilities $\bm{X}\sim\text{Dirichlet}(\bm{\theta})$.

It is useful to extend the definition of $r_{\bm{n}\mid\bm{l}}$ to allow some entries in $\bm{l}$ to be zero. Then $n_{i}=0$ if $l_{i}=0$ and

[TABLE]

In Theorem 4.1 the notation $\bm{n}\backslash ij$, $\bm{l}\backslash ij$ means the vectors with their $i$ and $j$th entries deleted. $r_{\bm{n}\backslash ij\mid\bm{l}\backslash ij}$ is equal to Eq. (9) with $l_{i}=0$, $l_{j}=0$ and $n\backslash ij\equiv n-n_{i}-n_{j}$. A calculation shows that $r_{\bm{n}\backslash ij\mid\bm{l}\backslash ij}$ is equal to the conditional probability $r_{\bm{n}\mid\bm{l}}/\mathbb{P}(n_{i},n_{j}\mid l_{i},l_{j})$ where the marginal distribution of $n_{i},n_{j}$ is

[TABLE]

An easy way to see this is to consider a Polya urn in which all colours other than $i$ and $j$ are considered identical. The probability that a particular edge, while $k$ edges in a coalescent tree, subtends $c$ edges in a sample of $n$ is a special case of Eq. (9)

[TABLE]

See also Griffiths and Tavaré (1988).

3 Sampling Distribution for Sampling Size $n=2$

To help understand Theorem 4.1, we first show how to calculate the probability of possible configurations when $n=2$. Looking back in time from $t$ to time zero there are $l=2$ founder lineages with probability $e^{-t}$ or $l=1$ founder lineages when the two initial lineages coalesce before $t$ with probability $1-e^{-t}$.

Consider first the contribution to the probability of a sample configuration $n_{a}=2$, at time $t$, when there are $l=2$ founder lineages. If we assume that $t=o(1/\theta)$, indicated by $t<<1/\theta$, then the probability that $l=2$ and that greater or equal to two mutations occur is e^{-t}\big{(}1-e^{-\theta t}-\theta te^{-\theta t}\big{)}=\mathcal{O}(\theta^{2}). For larger times, $t=\mathcal{O}(1/\theta)$, the probability $e^{-t}$ is $o(\theta)$, since $\lim_{\theta\to 0}e^{-1/\theta}/\theta=0$, and the $l=2$ contribution can safely be ignored to $\mathcal{O}(\theta)$. The probability that $l=2$ and $n_{a}=2$ is therefore

[TABLE]

In (12) there are two founder lineages at time $t$ with probability $e^{-t}$; the first term corresponds to the probability that both founder lineages are of type $a$ and no mutation occurs in $(0,t)$; the second term corresponds to either both lineages being of type $a$ and a mutation occurring from $a$ to $a$, not changing the type of the lineage, or one founder lineage is of type $a$ and the other is of type $j\neq a$ and a mutation occurs changing the lineage from type $j$ to type $a$. The probability of one mutation in $(0,t)$ on a single edge is $\frac{\theta}{2}te^{-\theta t}+o(\theta)$.

If $l=1$, the probability $1-e^{-t}$ that the two lineages coalesce is not small for large $t$, and care must be taken to include multiple mutations in the coalescent tree. Suppose coalescence occurs at time $w<t$ back. The two edges of the coalescent tree from the coalescent event at $t-w$ to the current time $t$ contribute a factor $e^{-w}$, which is $\mathcal{O}(1)$ for $w<<1/\theta$ and $o(\theta)$ for $w=\mathcal{O}(1/\theta)$. Thus it is sufficient to consider at most one mutation in this part of the tree. On the other hand, there is no restriction on the number of mutations on the single edge from time 0 to the coalescent event if the calculation is to be accurate to $\mathcal{O}(\theta)$ for large $t-w$, in particular, for $(t-w)\theta=\mathcal{O}(1)$. Writing the mutation rate matrix as $\gamma=\frac{\theta}{2}(P-I)$, the probability that $l=1$ and $n_{a}=2$ for all $t>0$ is therefore

[TABLE]

Probabilistically equation (LABEL:twol:1) is easiest to read from the first line, rather than from the simplified last line. In the first line the single founder is chosen to be type $j$ with probability $x_{j}$; the MRCA of the sample lineages is type $a$ from mutations along the single lineage; there could be no mutation along the two sample lineages, or one mutation from type $a$ to type $a$, not changing the type of the lineage. The integral accounts for the probability of coalescence of the two sample lineages at $w\in(0,t)$.

The total probability of a configuration $n_{a}=2$ is the sum of Eqs. (12) and (LABEL:twol:1). The stationary distribution is recovered by noting that $\bm{x}e^{t\gamma}\rightarrow\bm{\pi}$ as $t\rightarrow\infty$, giving

[TABLE]

If $t<<1/\theta$, then $e^{t\gamma}=I+\tfrac{1}{2}\theta t(P-I)+o(\theta)$, and Eqs. (12) and (LABEL:twol:1) give

[TABLE]

The other possible configuration of two types is $n_{a}=1$, $n_{b}=1$, $a\neq b$. If there are $l=2$ founder lineages, then by an argument similar to that for the $n_{a}=2$ case it is sufficient to consider at most one mutation in the tree. Then either the two founder lineages are of types $a$ and $b$ and there is no mutation along the edges; or one of the founder lineages is of type $a$ ($b$) one is of type $j$ and a mutation occurs from $j$ to $b$ ($a$) along the edge. The probability that $l=2$ and $n_{a}=1$, $n_{b}=1$ is therefore

[TABLE]

If there is $l=1$ founder lineage then by an argument similar to that for the $n_{a}=2$ case it is necessary to consider any number of mutations on the single edge of the coalescent tree from time 0 to the coalescence event, and sufficient to consider at most one mutation event on the two edges from the coalescence event to the present time $t$. The probability that $l=1$ and $n_{a}=1$, $n_{b}=1$ is therefore

[TABLE]

where $a\leftrightarrow b$ indicates the previous term with $a$ and $b$ interchanged. The total probability of a configuration $n_{a}=1$, $n_{b}=1$ is the sum of Eqs. (15) and (16). Taking the limit $t\rightarrow\infty$ we recover the stationary distribution

[TABLE]

If $t<<1/\theta$, then Eqs.  (15) and (16) give

[TABLE]

The above sampling distributions can be compared with an exact numerical calculation of the discrete-time, finite population, multi-allelic neutral Wright-Fisher model. Figure 2 shows the evolution of a $d=3$ allele neutral Wright-Fisher model with a population size $N=15$ and diffusion-limit instantaneous rate matrix

[TABLE]

with corresponding $\theta=2\max_{1\leq i\leq d}(-\gamma_{ii})=0.00720$. The conventions used for defining the mapping between the discrete Wright-Fisher model and its diffusion limit are set out in Burden and Tang (2016), Section 2. The initial distribution corresponding to an initial allele frequency $\mathbf{x}=(2,3,10)/15$ was evolved over discrete time steps $\delta t=1/N$ using the full $(N+1)(N+2)/2=136$ dimensional transition matrix defined in Eq. (1) of Burden and Tang (2016).

As is generally observed in population genetics models, the diffusion limit is approached rapidly for moderate population sizes. One observes that the ${\cal O}(\theta)$ approximation to the sampling distribution agrees well with the discrete Wright-Fisher model over both short and long time scales for all six possible allele combinations with $n=\sum_{a}n_{a}=2$. On the other hand the approximation in Eqs. (LABEL:twol:smalltheta) and (17) is only valid over short time scales as expected. This illustrates the importance of allowing for multiple mutations when there is a single founder lineage over long time scales in the small $\theta$ limit.

4 Sampling Distribution for Arbitrary Sampling Size

In a full $n$-coalescent tree the edge lengths $T_{n},\cdots,T_{2}$ are independent exponential random variables with rates ${n\choose 2},\ldots,{2\choose 2}$, however in a non-stationary partial coalescent tree back from time $t$ the edge lengths are constrained by the origin. If there are $l$ founder lineages then it is sometimes necessary to use q_{nl}(k;t)=\mathbb{E}\big{[}T_{k}\mathbb{I}\big{\{}A_{n}(t)=l\big{\}}\big{]}. Calculations of $q_{nl}(k;t)$ are made in Lemma 1.

Lemma 1

Explicit forms for expected edge lengths $T_{k}$ in a Kingman coalescent tree conditioned on $l$ founder lineages at time 0 are, for $1\leq l\leq k\leq n$ and $k\geq 2$,

[TABLE]

where

[TABLE]

and

[TABLE]

Proof

Calculations evaluate the expectation \mathbb{E}\big{[}T_{k}\mathbb{I}\big{\{}A_{n}(t)=l\big{\}}\big{]} using the identity that

[TABLE]

For $2\leq l<k\leq n$

[TABLE]

Eq. (20) is a convolution of $n-l+1$ functions with a Laplace transform, for $\phi\geq 0$,

[TABLE]

Therefore the integral is

[TABLE]

A modified argument is needed to show that the joint probability that $A_{n}(t)=l$ and mean length of edges while $l$ edges t-\big{(}T_{n}+\cdots+T_{l+1}\big{)} is actually the same as Eq. (21) when $k=l$. This expectation is

[TABLE]

Therefore Eq. (21) holds for $2\leq l\leq k\leq n$. $q_{n1}(1;t)$ follows from definition. ∎

The probability of a sample configuration for $n\geq 2$ is now calculated in Theorem 4.1. Combinatorics are more involved because it is necessary to consider the sample configuration from $l\leq n$ founder genes.

Theorem 4.1

The probability of a sample configuration $\bm{n}$ of $n$ individuals taken at time $t$ in the Wright-Fisher diffusion with generator in Eq. (1) as $\theta\to 0$ is

[TABLE]

where $r_{\bm{n}\mid\bm{l}},$ expressed in Eq. (8), is the probability of a configuration $\bm{n}$ in a partial coalescent tree from a configuration $\bm{l}\leq\bm{n}$;

[TABLE]

with q_{nl}(t):=\mathbb{P}\big{(}A_{n}(t)=l\big{)}; for $1\leq l\leq k\leq n$ and $k\geq 2$ (see Eq (6)); and q_{nl}(k;t):=\mathbb{E}\big{[}T_{k}\mathbb{I}\big{\{}A_{n}(t)=l\big{\}}\big{]} with explicit forms in Lemma 1. The Markov transition functions for mutation rate change are in the matrix $e^{\gamma t}$, where $\gamma=\frac{\theta}{2}(P-I)$; $m(\bm{l},l,\bm{x})$ is the multinomial distribution Eq. (2); and $p_{nk}(c)$ is defined in Eq. (11).

Proof

If there are no mutations occurring in the partial coalescent tree back to time 0 and the $l\geq 1$ founders have a configuration $\bm{l}$ then from a Polya urn argument considering $n-l$ branchings the probability of a configuration $\bm{n}$ at $t$ is $r_{\bm{n}\mid\bm{l}}$. The probability of a configuration of $\bm{n}$ when $l\geq 2$, and there is no mutation is

[TABLE]

The probability of a configuration $\bm{n}$, beginning from a configuration $\bm{l}$, when there is one mutation from $j$ to $i\neq j$ while $k$ edges is

[TABLE]

The argument used to obtain Eq. (24) is that the configuration from the $l$ founder edges to when $k$ edges; to after a mutation; then to when $n$ edges follows a Markov pattern $\bm{l}\to\bm{k}\to\bm{k}-\bm{e}_{j}+\bm{e}_{i}\to\bm{n}$. The sum in Eq. (24) can be simplified by summing over $\bm{k}$ keeping $k$, $k_{i}$ and $k_{j}$ fixed. Now by using Eq. (10)

[TABLE]

Summing over $\bm{k}\backslash ij$ with variables indexed by $i$ and $j$ held fixed, using Eq. (25),

[TABLE]

Finally the case when $l=1$ where coalescence to a single ancestor occurs before $t$ is evaluated. There is a need to consider mutations which may occur along the sample edges from the MRCA of the sample to the origin. If the MRCA occurs at time $w$ back from the current time then the type of the MRCA is determined by $\bm{x}e^{\gamma(t-w)}$. The probability of a configuration $\bm{n}=n\bm{e}_{a}$ when $l=1$ is therefore

[TABLE]

arguing that the MRCA occurs at time $w$ back, the identity of the MRCA is $a$ and there are either no mutations, or one from $a$ to $a$, to the MRCA.

The probability of a configuration $\bm{n}=n_{a}\bm{e}_{a}+n_{b}\bm{e}_{b}$ similarly needs a consideration of the type of the MRCA being either $a$ or $b$, and a mutation from the MRCA type to the other type while $k$ lineages. This probability is

[TABLE]

The total probability of a configuration of $\bm{n}$, Eq. (LABEL:config:3a) in the statement of the theorem, is found from Eq. (23), the simplified form of Eq. (24), and Eqs. (26) and (27). ∎

Corollary 1

(Burden and Griffiths, 2019)

The stationary sampling distribution is

[TABLE]

Proof

As $t\to\infty$ the MRCA is reached before $t$ with probability $1$, so only the last two terms in Eq. (LABEL:config:3a) need to be considered. The time to the MRCA $w$ is finite and $(\bm{x}e^{\gamma(t-w)})_{a}\to\pi_{a}$. The limit of the integrals in the term when $\bm{n}=n\bm{e}_{a}$ is

[TABLE]

and

[TABLE]

giving the first term in Eq. (LABEL:corr:11). The term in Eq. (LABEL:config:3a) where $\bm{n}=n_{a}\bm{e}_{a}+n_{b}\bm{e}_{b}$ converges to

[TABLE]

Simplification of the sum in Eq. (29) is done in Burden and Griffiths (2019). ∎

The combinatorial terms in Theorem 4.1 are much simpler if the initial frequency is $\bm{x}=\bm{e}_{a}$.

Corollary 2

Suppose the initial ancestor types are all $a$, so $\bm{x}=\bm{e}_{a}$. As $\theta\to 0$ the probability of a sample configuration where all individuals are type $a$ is

[TABLE]

the probability of a sample configuration with two types $a$ and $b$, $a\neq b$ is

[TABLE]

the probability of a sample configuration with two types $b$ and $c$, $a\neq b\neq c\neq a$, is

[TABLE]

and the probability of a sample configuration with greater than two types is of $o(\theta)$.

Proof

The probability of a configuration $\bm{n}=n\bm{e}_{a}$ is

[TABLE]

which is equal to Eq. (30) after expanding $e^{-\frac{\theta}{2}\sum_{k=2}^{n}kT_{k}}$ to $\mathcal{O}(\theta)$. The probability of a configuration $\bm{n}=n_{a}\bm{e}_{a}+n_{b}\bm{e}_{b}$ is

[TABLE]

which is simplifies to Eq. (31). A configuration $\bm{n}=n_{b}\bm{e}_{b}+n_{c}\bm{e}_{c}$, $b\neq c$, can only occur if the MRCA of the sample occurs before $t$ and the identity of the MRCA is $b$ or $c$. The probability of this sample configuration is

[TABLE]

which simplifies to Eq. (32). ∎

One may be interested only in time scales $t<<1/\theta$, in which case the following Corollaries 3 and 4 are relevant. If $t$ is fixed and $\theta\to 0$ then the probability of greater than one mutation along the edge from the MRCA to to the origin is $\mathcal{O}(\theta^{2})$ and the last two terms in Eq. (LABEL:config:3a) can be simplified.

Lemma 2

The following identities relating to the case $l=2$ will be of use in Corollaries 3 and 4:

[TABLE]

Proof

From Eq. (7) with $l=2$, $\int_{0}^{t}q_{n2}(x)dw=\int_{0}^{t}f_{n2}(x)dw=\mathbb{P}(T_{n}+\cdots+T_{2}<t)=\mathbb{P}(A_{n1}(t)=1)=q_{n1}(t)$. Alternatively one can differentiate Eq. (6) with $l=1$, and set the constant of integration by noting that $q_{n1}(0)=0$ for $n\geq 2$. For the second identity q_{n1}(k;t)=\mathbb{E}\big{[}T_{k}I\big{\{}A_{n}(t)=1\big{\}}\big{]}=\int_{0}^{t}\mathbb{E}\big{[}T_{k}I\big{\{}T_{n}+\cdots+T_{2}\in(w,w+dw)\big{\}}\big{]}=\int_{0}^{t}\mathbb{E}\big{[}T_{k}I\big{\{}A_{n}(w)=2\big{\}}\big{]}dw=\int_{0}^{t}q_{n2}(k;w)dw. The last identity is proved by noting that $\int_{0}^{t}q_{n2}(w)dw=q_{n1}(t)$, so $q_{n2}(t)=q^{\prime}_{n1}(t)$ and \int_{0}^{t}(t-w)q_{n2}(w)dw=\int_{0}^{t}(t-w)q^{\prime}_{n1}(w)dw=(t-w)q_{n1}(w)\Big{]}_{0}^{t}+\int_{0}^{t}q_{n1}(w)dw=0+q_{n1}(1;t). ∎

Corollary 3

If $t<<1/\theta$ as $\theta\to 0$ the last two terms in Eq. (LABEL:config:3a) become

[TABLE]

where expressions in Eq. (33) are defined previously.

Proof

In Eqs. (26) and (27) expand \bm{x}e^{\gamma(t-w)}=\bm{x}\big{(}I+\gamma(t-w)\Big{)}+\mathcal{O}(\theta^{2}). Then Eq. (26) becomes

[TABLE]

which simplifies to Eq. (33) after using the identities in Lemma 2. ∎

Corollary 4

Suppose the initial ancestor types are all $a$, so $\bm{x}=\bm{e}_{a}$. If $t<<1/\theta$ the probability of a sample configuration where all individuals are type $a$ is

[TABLE]

the probability of a sample configuration with two types $a$ and $b$ is

[TABLE]

and the probability of a sample configuration with greater than two types is of $\mathcal{O}(\theta^{2})$.

Proof

The last line of Eq. (30) is equal to

[TABLE]

Now add the first line of Eq. (30) to obtain Eq. (34). The probability of a configuration $\bm{n}=n_{a}\bm{e}_{a}+n_{b}\bm{e}_{b}$, Eq. (35), is found by noting that $(e^{\gamma(t-w)})_{ab}=\delta_{ab}+\mathcal{O}(\theta)$ so the second last line of Eq. (31) is

[TABLE]

and the last line is $\mathcal{O}(\theta^{2})$. Then Eq. (35) follows by adding Eq. (36) and the first line of Eq. (31). ∎

5 The infinitely-many-sites model and gene trees

In the infinitely-many-sites model genes are sequences. Mutations occur at positions never before mutant in the population of sequences. The configuration of mutations at segregating sites in a sample is equivalent to a unique gene tree, called a perfect phylogeny in biology. A gene tree of $n$ sequences of $g\leq d$ different types is constructed by labeling mutations in the coalescent tree of the $n$ genes, then describing each gene by its mutation path $\bm{z}$ from leaves to the root of the tree with multiplicities of the types $\bm{n}=(n_{1},\ldots,n_{g})$. The same ideas used in Theorem 1 lead to an approximate probability for a gene tree when $\theta$ is small. It is natural to consider gene trees because they are closely related to coalescent trees. The gene tree does not depend on the actual values of the labels, only their pattern. This is a very brief account; the reader is referred to Ethier and Griffiths (1997); Griffiths (1999); Griffiths and Tavaré (1999) for more details. An illustration of a coalescent tree showing these concepts is in Fig. 3.

An analog of the result in Theorem 1 for gene trees is

[TABLE]

The initial population of sequences has frequencies $\bm{x}$ which may have an infinite number of entries. In the labeling of $\bm{l}$ sample sequences only the non-zero number of sequences are included and relabeled.

The first line of (37) corresponds to the case where no mutation occurs in $(0,t)$. In this case the gene trees at time zero and time $t$, specified as above in terms of paths to the root, are equal, that is $T(\bm{l})=T$. However the trees have different multiplicities $\bm{l}$ and $\bm{n}$, and the argument leading to Eq.(23) applies. In the second line one mutation occurs adding a new singleton sequence to the tree. $T_{j}(\bm{l})$ is the gene tree $T(\bm{l})$ modified by mutation to one of the $j$th sequences having an additional mutation label added. In the sum over $\bf{k}$, there are two cases to consider; if $k_{j}>1$ there are still representatives of the $j$th sequence; and if $k_{j}=1$ the original sequence is removed. The new sequence is a singleton. Its position, denoted by $*$ is either at the end of the sequences, or in the position of the prior $j$th sequence which was removed in the second case.

The third and fourth lines correspond to cases where the coalescent tree has a MRCA in $(0,t)$. Then either the gene tree is a single sequence \big{(}(z_{0})\big{)} of $n$ copies at $t$ if there is no further mutation, or there is a mutation and then the gene tree at $t$ is T=\big{(}(z_{1}),(z_{0},z_{1})\big{)} with multiplicities $(n_{1},n_{2})$. In the third line T=\big{(}(z_{0})\big{)} is a tree which is the same type as a single gene chosen from the initial sequences $\bm{x}$. In the fourth line T=\big{(}(z_{1}),(z_{0},z_{1})\big{)} does not distinguish whether or not $z_{1}$ is from the initial population, or a new mutation. If there is no mutation along the single lineage it will be from the initial population, or if there is mutation it will denote a new mutation. To split the terms up into these two cases replace $q_{n1}(k;t)$ by $\int_{0}^{t}q_{n2}(k;w)e^{-\frac{\theta}{2}(t-w)}dw$ for no mutation and \int_{0}^{t}q_{n2}(k;w)\big{(}1-e^{-\frac{\theta}{2}(t-w)}\big{)}dw if mutation occurs.

The infinitely-many-alleles model is a model of gene frequencies where mutation always produces a new allele type. The distribution of the configuration of types in a sample of $n$ genes is the Ewens’ sampling formula (Ewens, 1972). The infinitely-many-alleles model is contained in the infinitely-many-sites model by considering the first co-ordinates of the sequences in the gene trees. The two models can be thought of at a macroscopic and microscopic level of detail in DNA sequences. There is a similar equation to (37) for the infinitely-many-alleles model, keeping the types labeled for convenience.

6 An approximate transition density for large sample sizes

Denote by $f(\bm{x},\bm{y};t,\theta)$ the transition density of the Wright-Fisher diffusion process $\{\bm{X}(t)\}_{t\geq 0}$ with initial frequencies $\bm{X}(0)=\bm{x}$ and final frequencies $\bm{X}(t)=\bm{y}$. The transition density $f(\bm{x},\bm{y},t)\equiv f(\bm{x},\bm{y};t,0)$ in a model with no mutation. For functions $g(\bm{X}(t))$ with finite expectation with $\theta\geq 0$

[TABLE]

where $\mathbb{E}^{\theta}_{\bm{x}}$ denotes expectation with parameter $\theta$ and initial frequencies $\bm{x}$.

An expression for the transition density as a mixture of Dirichlet distributions is known from Ethier and Griffiths (1993) and discussed in the review of Griffiths and Spanó (2010). A description of the structure of the population at $t$ and this transition density is now given. The population at time $t$ coalesces back to a random number of lineages $l$ at time 0, with probability $q_{\infty l}(t)$. The population is composed of families of individuals from these $l$ founders, with family sizes having an $l$ dimensional Dirichlet distribution with parameters $(1,\cdots,1)$. Family types are determined by the type of the founder lineages, so there is a $d$-dimensional $m(\bm{l},l,\bm{x})$ distribution of the accumulated number of families with types in $[d]$. The probability that in a subset $A\subseteq[d]$, $X_{j}(t)>0$ for $j\in A$ and $X_{j}(t)=0$ for $j\not\in A$ is the probability the set of founder lineage types is $A$. The transition density is then

[TABLE]

with the convention in the Dirichlet distribution that $y_{j}=0$ if $l_{j}=0$. That is

[TABLE]

Although the first summation is formally written beginning at 1, because of the convention that all the terms with $l<|A|$ are zero the summation really begins with $l=|A|$. See Ethier and Griffiths (1993) Eq. (1.27) for a more general expansion which includes Eq. (38) when their $\theta=0$. A first example: if $A=\{1,2\}$, so $y_{1}>0,y_{2}>0$ and $y_{1}+y_{2}=1$, the transition density is

[TABLE]

which is the joint density of $X_{1}(t),X_{2}(t)$ and the probability that $X_{j}(t)=0,j>2$. A second example: if $\bm{y}=\bm{e}_{1}$ the probability the population is fixed at time $t$ with type 1 individuals is the probability that all founder lineages are of type 1,

[TABLE]

It would be ideal to find an expansion of $f(\bm{x},\bm{y};t,\theta)$ in terms of $f(\bm{x},\bm{y};t,0)$ to the first order in $\theta$ similar to the sample transition distribution, however for fixed $\theta$ the probability of greater than one mutation in an infinite-leaf coalescent tree is one, not of order $\theta^{2}$. The error for a probability configuration in a sample of $n$ by omitting the probability of greater than one mutation (denoted by $Q_{1}$) is bounded by the probability of greater than one mutation in an unconstrained coalescent tree.

This probability is exactly

[TABLE]

Both terms in (40) are of order $\theta^{2}$. To see this for the first term note that

[TABLE]

Alternatively use the properties of the gamma function that as $\theta\to 0$

[TABLE]

where $\gamma$ is Euler’s constant and $\psi(\cdot)$ is the diagamma function with

[TABLE]

to obtain the result.

As $n\to\infty$ and $\theta\to 0$, there are three possible limits for $Q_{1}$, seen from the second line of (40);

[TABLE]

If $n\to\infty$ and $\theta\to 0$ such that $\theta\log n\to 0$ then the probability of greater than one mutation in an $n$-coalescent tree is $o(\theta)$; we take this approach and find an approximate transition density $f_{n}(\bm{x},\bm{y};t,\theta)$.

Theorem 6.1

An approximate transition density in the general mutation Wright-Fisher diffusion from the sample transition distribution as $n\to\infty$ and $\theta\to 0$ such that $\theta\log n\to 0$, and $\bm{n}/n\to\bm{y}$, with initial relative frequencies $\bm{x}$ is given by

[TABLE]

Proof

The proof follows immediately by taking $n\to\infty$ such that $\bm{n}/n\to\bm{y}$ in $p(\bm{n};t,\bm{x})$ of Theorem 4.1. The transition density for $\bm{X}(t)$ can include corners of the space where some variables are zero. Inversion of the sampling distributions in the expansion Eq. (LABEL:config:3a) follows from;

[TABLE]

The generalised Dirichlet distributions ${\cal D}_{\bm{l}}(\bm{y})$, ${\cal D}_{\bm{l}\backslash ij}(\bm{y}\backslash ij)$ have the property that if variables $\bm{y}$ or $\bm{y}\backslash ij$ are zero, the resulting Dirichlet distribution is of lower dimension than $d$. The transition density can include corners of the space where some of the variables are zero if some entries in $\bm{x}$ are zero.

The condition that $\theta\log n\to 0$ is needed to ensure that as $\theta\to 0$, $n\to\infty$ then $\frac{\theta}{2}\sum_{k=l}^{n}kq_{nl}(k;t)\to 0$ as well as other similar terms in (42). This is so because

[TABLE]

In Burden and Tang (2016); Burden and Griffiths (2019) the way a stationary distribution is developed is by finding that if $\theta$ is small, then the population is approximately fixed or contains just two allele types $a$, $b$ with relative frequencies on a line $y_{a}+y_{b}=1$. This is a boundary approximation. The transition density (42) when there are two types is

[TABLE]

The sum in (44) converges replacing $n$ by $\infty$ from (42) because for $0<y_{b}<1$ (considering the first term)

[TABLE]

The limit stationary distribution in (44) as $t\to\infty$ is

[TABLE]

on the line $y_{a}+y_{b}=1$, as in Burden and Tang (2016); Burden and Griffiths (2019), because $(e^{\gamma(t-w)})_{a}\to\pi_{a}$ and \int_{0}^{\infty}q_{\infty 2}(k;w)dw=\mathbb{E}\big{[}T_{k}\big{]}={k\choose 2}^{-1}.

7 Discussion

There is a natural curiosity in extending the stationary sampling distribution with small mutation rates (5), derived in Burden and Griffiths (2019), to the sampling distribution taken at time $t$ using coalescent methods. An important point is that because of the duality of the coalescent process with the Wright-Fisher diffusion this is the same as the sampling distribution at time $t$ in a Wright-Fisher diffusion. In this paper we find expressions for the sampling distribution in a sample of $n$ genes taken at time $t$ from a Wright-Fisher population which follows a diffusion model. The overall mutation rate $\theta$ is taken to be small and mutation rates between types are general. An extension to gene trees is described. The idea in obtaining formulae is simple and is based on the result that to order $\theta$ it is only necessary to consider at most one mutation in sample lineages up to coalescence. A sample then consists of families of genes from founder lineages which may be of different types and at most one family of types descendent from a mutation in the coalescent tree. An illustration of these gene families is in Fig. 1. There is a combinatorial complexity in Theorem 1 because of the different possible types at time zero. If the population at time zero consists of just one type, expressions for the sampling probability are much simpler as shown in Corollary 3. In a sample size of two the concepts are much easier, so that case is worked through first. A simulation for sample size two shows the approximation to be very accurate for small $\theta$.

Many authors have observed that samples at Single Nucleotide Polymorphism (SNP) sites with small mutation rates are either fixed, or contain bi-allelic states. If a site is fixed initially then explicit calculations in Corollary 2 and Corollary 4 model this feature of bi-allelic states. An application of the theory in this paper is to SNP positions in DNA sequences where mutation rates are small. Because the rate matrix $\gamma$ is constrained only by the requirement that its off-diagonal elements be ${\mathcal{O}}(\theta)$, the theory can accommodate any model class of nucleotide substitution rate matrices (see, for instance Yang, 2006, Table 1.1).

The frequency spectrum in DNA sequence data is often of interest. As an example we consider a frequency spectrum in which the ancestral bases are fixed at each individual SNP site, with relative proportions $z_{A},z_{T},z_{C},z_{G}$ along the sequence sites. If $t<<1/\theta$ then Corollary 4 applies at each SNP site and Eq. (35) shows that the probability of seeing $n_{A}$ entries $A$ at a site and $n_{T}$ entries $T$ when the ancestor type is $A$ is, to $\mathcal{O}(\theta^{2})$,

[TABLE]

with similar expressions for other pairs. The expected proportion of sites which are SNPs for which there are $n_{A}$ entries $A$ and $n_{T}$ entries $T$ at the site when the ancestor base is not known, denoted by $f(n_{A},n_{T})$, is then

[TABLE]

In a Bayesian approach a prior distribution could be chosen for $z_{A},z_{T},z_{C},z_{G}$. The derivation of $f(n_{A},n_{T})$ does not need independence between sites, which may be needed for further inference. If independence of SNP sites is assumed, then the likelihood of observing SNP site counts of $r(n_{A},n_{T})$ with $n_{A}$ entries $A$ and $n_{T}$ entries $T$ is $f(n_{A},n_{T})^{r(n_{A},n_{T})}$. Using a similar notation for the other base pairs the likelihood from segregating sites in the full data is, up to a multinomial coefficient,

[TABLE]

There is information in the non-segregating sites as well, which could be included in the likelihood. This is a sketch of a likelihood application. Full details are beyond the scope of this paper. Rate matrix estimation from site frequency data assuming stationarity is described in some detail in Burden and Tang (2017).

There are several reasons to consider large sample sizes by taking $n\to\infty$. The first is that current sample sizes can be large. The second is that the approximate sampling distributions are obtained in probability to $o(\theta)$ for fixed sample size $n$, and the error can possibly be large when $n\to\infty$. We show from (40) that when $n\to\infty$ this order of approximation can still be small if $n\to\infty$ and $\theta\to 0$ such that $n\theta\to\alpha$, where $\alpha<<1$. The third reason is that ideally we would like to use the coalescent approach to find an approximate transition density in the Wright-Fisher diffusion for small mutation rates. The approach taken in a sample cannot be used directly in an infinite-leaf coalescent tree because as $\theta\to 0$ the probability of greater than one mutation in the coalescent tree is 1, not $\mathcal{O}(\theta^{2})$. Nevertheless we can think of a large sample size model as approximating a discrete Wright-Fisher model with an effective population size of $N$. Then it is appropriate to consider sampling formulae when $n\to\infty$ and $\theta\to 0$ such that $n\theta\to\alpha<<1$ by thinking of $N=n$. $\alpha<<1$ is likely to be satisfied in applications, thinking of this as a population limit. Mutation rates at neutral genomic sites are typically less than $10^{-7}$ per generation per base, effective population sizes $N\sim\mathcal{O}(10^{4})$, so $\theta$ is typically $\mathcal{O}(10^{-3})$, and then $\alpha$ is typically less than $0.01$, which gives the asymptotic probability of more than one mutation in an infinite coalescent tree limit, Eq. (41), as $Q_{1}\sim 1-e^{-\alpha}(1+\alpha)<0.00005$.

8 Acknowledgement

We thank two referees for a careful reading of the manuscript and their comments and suggestions which have improved the paper.