PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Ann Appl Probab. Author manuscript; available in PMC 2010 July 28.
Published in final edited form as:
Ann Appl Probab. 2010 June; 20(3): 1005–1028.
doi:  10.1214/09-AAP646
PMCID: PMC2910927
NIHMSID: NIHMS152909

AN ASYMPTOTIC SAMPLING FORMULA FOR THE COALESCENT WITH RECOMBINATION

Abstract

Ewens sampling formula (ESF) is a one-parameter family of probability distributions with a number of intriguing combinatorial connections. This elegant closed-form formula first arose in biology as the stationary probability distribution of a sample configuration at one locus under the infinite-alleles model of mutation. Since its discovery in the early 1970s, the ESF has been used in various biological applications, and has sparked several interesting mathematical generalizations. In the population genetics community, extending the underlying random-mating model to include recombination has received much attention in the past, but no general closed-form sampling formula is currently known even for the simplest extension, that is, a model with two loci. In this paper, we show that it is possible to obtain useful closed-form results in the case the population-scaled recombination rate ρ is large but not necessarily infinite. Specifically, we consider an asymptotic expansion of the two-locus sampling formula in inverse powers of ρ and obtain closed-form expressions for the first few terms in the expansion. Our asymptotic sampling formula applies to arbitrary sample sizes and configurations.

Keywords and phrases: Ewens sampling formula, coalescent process, recombination, two-locus model, infinite alleles, Golding’s recursion

1. Introduction

The probability of a sample configuration provides a useful ground for analyzing genetic data. Popular applications include obtaining maximum likelihood estimates of model parameters and performing ancestral inference [see Stephens (2001)]. In principle, model-based full-likelihood analyses, such as that based on the coalescent [Kingman (1982a,b)], should be among the most powerful methods since they make full use of the data. However, in most cases it is intractable to obtain a closed-form formula for the probability of a given dataset. A well-known exception to this hurdle is the Ewens sampling formula (ESF), which describes the stationary probability distribution of a sample configuration under the one-locus infinite-alleles model in the diffusion limit [Ewens (1972)]. Notable biological applications of this closed-form formula include the test of selective neutrality [see Slatkin (1994); Watterson (1977)]. Hoppe (1984) provided a Pólya-like urn model interpretation of the formula, and recently Griffiths and Lessard (2005) provided a new combinatorial proof of the ESF and extended the framework to obtain new results for the case with a variable population size. We refer the reader to the latter paper for a nice summary of previous works related to the ESF. Note that the ESF also arises in several interesting contexts outside biology, including random partition structures and Bayesian statistics; see Arratia et al. (2003) for examples of intricate combinatorial connections. The ESF is a special case of the two-parameter sampling formula constructed by Pitman (1992, 1995) for exchangeable random partitions.

Golding (1984) considered generalizing the infinite-alleles model to include recombination and constructed a recursion relation satisfied by the two-locus sampling probability distribution at stationarity in the diffusion limit. Ethier and Griffiths (1990) later undertook a more mathematical analysis of the two-locus model and provided several interesting results. However, to date a general closed-form formula for the two-locus sampling distribution remains unknown. Indeed, it is widely recognized that recombination adds a formidably challenging layer of complexity to population genetics analysis. Because obtaining exact analytic results in the presence of recombination is difficult, recent research has focused on developing sophisticated and computationally-intensive Monte Carlo techniques. Examples of such techniques applied to the coalescent include Monte Carlo simulations [see Hudson (1985, 2001)], importance sampling [see De Iorio and Griffiths (2004a,b); Fearnhead and Donnelly (2001); Griffiths and Marjoram (1996); Griffiths et al. (2008); Stephens and Donnelly (2000)], and Markov Chain Monte Carlo methods [see Kuhner et al. (2000); Nielsen (2000); Wang and Rannala (2008)].

Being the simplest model with recombination, the two-locus case has been extensively studied in the past [Ethier and Griffiths (1990); Golding (1984); Griffiths (1981, 1991); Hudson (1985)], and a renewed wave of interest was recently sparked by Hudson (2001), who proposed a composite likelihood method which uses two-locus sampling probabilities as building blocks. LD-hat, a widely-used software package for estimating recombination rates, is based on this composite likelihood approach, and it has been used to produce a fine-scale map of recombination rate variation in the human genome [McVean et al. (2004); Myers et al. (2005)]. LDhat assumes a symmetric di-allelic recurrent mutation model at each locus and relies on the importance sampling scheme proposed by Fearnhead and Donnelly (2001) for the coalescent with recombination, to generate exhaustive lookup tables containing two-locus probabilities for all inequivalent sample configurations and a range of relevant parameter values. This process of generating exhaustive lookup tables is very computationally expensive. A fast and accurate method of estimating two-locus probabilities would be of practical value.

In this paper, we revisit the tantalizing open question of whether a closed-form sampling formula can be found for the coalescent with recombination. We show that, at least for the two-locus infinite-alleles model with the population-scaled recombination rate ρ large but not necessarily infinite, it is possible to obtain useful closed-form analytic results. Note that the aforementioned Monte Carlo methods generally become less efficient as ρ increases. Those methods involve sampling a large collection of genealogical histories consistent with the observed sample configuration, and, when ρ is large, the sampled genealogies tend to be very complicated; they typically contain many recombination events, and it may take a long time for every locus to reach a most recent common ancestor. However, contrary to this increased complexity in the standard coalescent, we actually expect the evolutionary dynamics to be easier to describe for large ρ, since the loci under consideration would then be less dependent. Hence, it seems reasonable to conjecture that there may exist a stochastic process simpler than the standard coalescent with recombination that describes the relevant degrees of freedom in the large ρ limit. We believe that our sampling formula may provide some hints as to what that dual process should be.

The work discussed here generalizes previous results [Ethier and Griffiths (1990); Golding (1984)] for ρ = , in which case the loci become independent and the two-locus sampling distribution is given by a product of one-locus ESFs. Our main results can be summarized as follows.

Main results

Consider the diffusion limit of the two-locus infinite-alleles model with population-scaled mutation rates θA and θB at the two loci. For a sample configuration n (defined later in the text), we use q(n|θA, θB, ρ) to denote the probability of observing n given the parameters θA, θB, and ρ. For an arbitrary n, our goal is to find an asymptotic expansion of q(n|θA, θB, ρ) in inverse powers of ρ, i.e., for large values of the recombination rate ρ, our goal is to find

q(nθA,θB,ρ)=q0(nθA,θB)+q1(nθA,θB)ρ+q2(nθA,θB)ρ2+O(1ρ3),

where q0, q1, and q2 are independent of ρ. As mentioned before, q0(n|θA, θB) is given by a product of one-locus ESFs. In this paper, we derive a closed-form formula for the first-order term q1(n|θA, θB). Further, we show that the second-order term q2(n|θA, θB) can be decomposed into two parts, one for which we obtain a closed-form formula and the other that satisfies a simple strict recursion. The latter can be easily evaluated using dynamic programming. Details of these results are described in Section 3. In a similar vein, in Section 4 we obtain a simple asymptotic formula for the joint probability distribution of the number of alleles observed at the two loci.

We remark that our work has practical value in genetic analysis. While this paper was under review, we applied the technique developed here to obtain analogous results for an arbitrary finite-alleles recurrent mutation model. See Jenkins and Song (2009) for details. In that paper, we performed an extensive assessment of the accuracy of our results for a particular finite-alleles model of mutation, and showed that they may be accurate even for moderate values of ρ, including a range that is of biological interest. The accuracy (not discussed here) of our results for the infinite-alleles model is very similar to that finite-alleles case.

2. Preliminaries

In this section, we review the ESF for the one-locus infinite-alleles model, as well as Golding’s (1984) recursion relation for the two-locus generalization. Our notational convention generally follows that of Ethier and Griffiths (1990).

Given a positive integer k, [k] denotes the k-set {1, …, k}. For a nonnegative real number x and a positive integer n, (x)n := x(x + 1) … (x + n − 1) denotes the nth ascending factorial of x. We use 0 to denote either a vector or a matrix of all zeroes; it will be clear from context which is intended. Throughout, we consider the diffusion limit of a neutral haploid exchangeable model of random mating with constant population size 2N. We refer to the haploid individuals in the population as gametes.

2.1. Ewens sampling formula for the one-locus model

In the one-locus model, a sample configuration is denoted by a vector of multiplicities n = (n1, …, nK), where ni denotes the number of gametes with allele i at the locus and K denotes the total number of distinct allelic types observed. We use n to denote i=1Kni, the total sample size. Under the infinite-alleles model, any two gametes can be compared to determine whether or not they have the same allele, but it is not possible to determine how the alleles are related when they are different. Therefore, allelic label is arbitrary. The probability of a mutation event at the locus per gamete per generation is denoted by u. In the diffusion limit, N → ∞ and u → 0 with the population-scaled mutation rate θ = 4Nu held fixed. Each mutation gives rise to a new allele that has never been seen before in the population. For the one-locus model just described, Ewens (1972) obtained the following result:

Proposition 2.1 (Ewens)

At stationarity in the diffusion limit of the one-locus infinite-alleles model with the scaled mutation parameter θ, the probability of an unordered sample configuration n = (n1, …, nK) is given by

p(nθ)=n!n1nK1α1!αn!θK(θ)n,
(2.1)

where αi denotes the number of allele types represented i times, i.e., αi := |{k | nk = i}|.

Let An external file that holds a picture, illustration, etc.
Object name is nihms152909ig1.jpg denote an ordered configuration of n sequentially sampled gametes such that the corresponding unordered configuration is given by n. By exchangeability, the probability of An external file that holds a picture, illustration, etc.
Object name is nihms152909ig1.jpg is invariant under all permutations of the sampling order. Hence, we can write this probability of an ordered sample as q(n) without ambiguity. It is given by

q(nθ)=p(nθ)[n!i=1Kni!1α1!αn!]1=[i=1K(ni1)!]θK(θ)n,
(2.2)

which follows from the fact that there are n!i=1Kni!1α1!αn! orderings corresponding to n [Hoppe (1984)]. To understand what we mean by ordered and unordered sampling configurations, it is helpful to relate the Ewens sampling formula to the theory of random partitions. If the gametes are labeled in order of appearance by 1, …, n, then the resulting sample configuration de-fines a random partition of [n], with gametes belonging to the same block if and only if they have the same allele. The quantity q(n | θ) is then the probability of a particular partition of [n] whose block sizes are given by the entries in n, while the quantity p(n | θ) is the probability of observing any partition of [n] with these block sizes. For example, if n = (2, 1, 1), then there are six partitions of [4] with these block sizes, and so p(n | θ) = 6q(n | θ). It is often more convenient to work with an ordered sample than with an unordered sample. In this paper, we will work with the former; i.e., we will work with q(n | θ) rather than p(n | θ).

In the coalescent process going backwards in time, at each event a lineage is lost either by coalescence or mutation. By consideration of the most recent event back in time, one can show that q(n | θ) satisfies

n(n1+θ)q(nθ)=i=1Kni(ni1)q(neiθ)+θi=1Kδni,1q(neiθ),
(2.3)

where δni,1 is the Kronecker delta and ei is a unit vector with the ith entry equal to one and all other entries equal to zero. The boundary condition is q(ei |θ) = 1 for all i [set membership] [K], and q(n |θ) is defined to be zero if n contains any negative component. It can be easily verified that the formula of q(n |θ) shown in (2.2) satisfies the recursion (2.3).

Ewens (1972) also obtained the following result regarding the number of allelic types:

Proposition 2.2 (Ewens)

Let Kn denote the number of distinct allelic types observed in a sample of size n. Then,

P(Kn=kθ)=s(n,k)θk(θ)n,
(2.4)

where s(n, k) are the unsigned Stirling numbers of the first kind. Note that (θ)n = s(n, 1)θ+ s(n, 2)θ2 + … + s(n, n)θn.

It follows from (2.1) and (2.4) that Kn is a sufficient statistic for θ.

2.2. Golding’s recursion for the two-locus case

Golding (1984) first generalized the one-locus recursion (2.3) to two loci, and Ethier and Griffiths (1990) later undertook a more mathematical study of the model. We denote the two loci by A and B, and use θA and θB to denote the respective population-scaled mutation rates. We use K and L to denote the number of distinct allelic types observed at locus A and locus B, respectively. The population-scaled recombination rate is denoted by ρ = 4Nr, where r is the probability of a recombination event between the two loci per gamete per generation. A key observation is that to obtain a closed system of equations, the type space must be extended to allow some gametes to be specified only at one of the two loci.

Definition 2.1 (Extended sample configuration for two loci)

The two-locus sample configuration is denoted by n = (a, b, c), where a = (a1, …, aK) with ai being the number of gametes with allele i at locus A and unspecified alleles at locus B, b = (b1, …, bL) with bj being the number of gametes with unspecified alleles at locus A and allele j at locus B, and c = (cij) is a K × L matrix with cij being the multiplicity of gametes with allele i at locus A and allele j at locus B. Further, we define

a=i=1Kai,ci·=j=1Lcij,c=i=1Kj=1Lcij,b=j=1Lbj,c·j=i=1Kcij,n=a+b+c.

We use q(a, b, c) to denote the sampling probability of an ordered sample with configuration (a, b, c). For ease of notation, we do not show the dependence on parameters. For 0 ≤ ρ < ∞, Golding’s (1984) recursion for q(a, b, c) takes the following form:

[n(n1)+θA(a+c)+θB(b+c)+ρc]q(a,b,c)=i=1Kai(ai1+2ci·)q(aei,b,c)+j=1Lbj(bj1+2c·j)q(a,bej,c)+i=1Kj=1L[cij(cij1)q(a,b,ceij)+2aibjq(aei,bej,c+eij)]+θAi=1K[j=1Lδai+ci·,1δcij,1q(a,b+ej,ceij)+δai,1δci·,0q(aei,b,c)]+θBj=1L[i=1Kδbj+c·j,1δcij,1q(a+ei,b,ceij)+δbj,1δc·j,0q(a,bej,c)]+ρi=1Kj=1Lcijq(a+ei,b+ej,ceij).
(2.5)

Relevant boundary conditions are q(ei, 0, 0) = q(0, ej, 0) = 1 for all i [set membership] [K] and j [set membership] [L]. For notational convenience, we deviate from Ethier and Griffiths (1990) and allow each summation to range over all allelic types. To be consistent, we define q(a, b, c) = 0 whenever any entry in a, b, or c is negative.

For ease of discussion, we define the following terms:

Definition 2.2 (Degree)

The degree of q(a, b, c) is defined to be a + b + 2c.

Definition 2.3 (Strictly recursive)

We say that a recursion relation is strictly recursive if it contains only a single term of the highest degree.

Except in the special case ρ = , a closed-form solution for q(a, b, c) is not known. Notice that the terms q(aei, bej, c + eij) and q(a + ei, b + ej, ceij) on the right-hand side of (2.5) have the same degree as q(a, b, c) on the left-hand side. Therefore, (2.5) is not strictly recursive. For each degree, we therefore need to solve a system of coupled equations, and this system grows very rapidly with n. For example, for a sample with a = 0, b = 0, and c = 40, computing q(0, 0, c) requires solving a system of more than 20, 000 coupled equations [Hudson (2001)]; this is around the limit of sample sizes that can be handled in a reasonable time. In the following section, we revisit the problem of obtaining a closed-form formula for q(a, b, c) and obtain an asymptotic expansion for large ρ.

3. An asymptotic sampling formula for the two-locus case

For large ρ, our objective is to find an asymptotic expansion of the form

q(a,b,c)=q0(a,b,c)+q1(a,b,c)ρ+q2(a,b,c)ρ2+O(1ρ3),
(3.1)

where q0, q1, and q2 are independent of ρ. Our closed-form formulas will be expressed using the following notation:

Definition 3.1

For a given multiplicity vector a = (a1, …, aK) with a=i=1Kai, we define

qA(a)=[i=1K(ai1)!]θAK(θA)a.
(3.2)

Similarly, for a given multiplicity vector b = (b1, …, bL) with b=i=1Lbi, we define

qB(b)=[j=1L(bj1)!]θBL(θB)b.
(3.3)

As discussed in Section 2.1, qA (respectively, qB) gives the probability of an ordered sample taken from locus A (respectively, B).

Definition 3.2 (Marginal configuration)

We use cA = (c)i[set membership][K] and cB = (c·j)j[set membership][L] to denote the marginal sample configurations of c restricted to locus A and locus B, respectively.

The leading order term q0(a, b, c) is equal to q(a, b, c) when ρ = , in which case the two loci are independent. Theorem 2.3 of Ethier and Griffiths (1990) states that q0(0, 0, c) = qA(cA)qB(cB). More generally, one can obtain the following result for the leading order contribution:

Proposition 3.1

In the asymptotic expansion (3.1) of the two-locus sampling formula, the zeroth order term q0(a, b, c) is given by

q0(a,b,c)=qA(a+cA)qB(b+cB).
(3.4)

Although this result is intuitively obvious, in Section 5.1 we provide a detailed new proof, since it well illustrates our general strategy. One of the main results of this paper is a closed-form formula for the next order term q1(a, b, c). The case with c = 0 admits a particularly simple solution:

Lemma 3.1

In the asymptotic expansion (3.1) of the two-locus sampling formula, the first order term satisfies

q1(a,b,0)=0

for arbitrary a and b.

That q1(a, b, 0) vanishes is not expected a priori. Below we shall see that q2(a, b, 0) ≠ 0 in general. For an arbitrary configuration matrix c of nonnegative integers, we obtain the following closed-form formula for q1(a, b, c):

Theorem 3.1

In the asymptotic expansion (3.1) of the two-locus sampling formula, the first order term q1(a, b, c) is given by

q1(a,b,c)=(c2)qA(a+cA)qB(b+cB)qB(b+cB)i=1K(ci·2)qA(a+cAei)qA(a+cA)j=1L(c·j2)qB(b+cBej)+i=1Kj=1L(cij2)qA(a+cAei)qB(b+cBej),
(3.5)

for arbitrary configurations a, b, c of non-negative integers.

Lemma 3.1 is used in proving Theorem 3.1. A proof of Theorem 3.1 is provided in Section 5.2, while a proof of Lemma 3.1 is given in Section 5.3. Note that the functional form of q0(a, b, c) and q1(a, b, c) in (3.4) and (3.5) has no explicit dependence on mutation; i.e., the dependence on mutation is completely absorbed into the marginal one-locus probabilities. It turns out that (3.4) and (3.5) are universal in that they also apply to an arbitrary finite-alleles model of mutation, with qA and qB replaced with appropriate marginal one-locus probabilities for the assumed mutation model. See Jenkins and Song (2009) for details.

In principle, similar arguments can be used to find the (j + 1)th order term given the jth, although a general expression does not seem to be easy to obtain. In Section 5.4, we provide a proof of the following result for q2(a, b, c):

Theorem 3.2

In the asymptotic expansion (3.1) of the two-locus sampling formula, the second order term q2(a, b, c) is of the form

q2(a,b,c)=q2(a+cA,b+cB,0)+σ(a,b,c),
(3.6)

where σ(a, b, c) is given by the analytic formula shown in Appendix A, and q2(a, b, 0) satisfies the following strict recursion:

[a(a+θA1)+b(b+θB1)]q2(a,b,0)=i=1Kai(ai1)q2(aei,b,0)+j=1Lbj(bj1)q2(a,bej,0)+θAi=1Kδai,1q2(aei,b,0)+θBj=1Lδbj,1q2(a,bej,0)+4[aθA(θA+a1)i=1Kδai,1][bθB(θB+b1)j=1Lδbj,1]qA(a)qB(b),
(3.7)

with boundary conditions q2(ei, 0, 0) = q2(0, ej, 0) = 0 for all i [set membership] [K] and j [set membership] [L].

In contrast to q1(a, b, 0) (c.f., Lemma 3.1), it turns out that q2(a, b, 0) does not vanish in general. We do not have an analytic solution for q2(a, b, 0), but note that (3.7) is strictly recursive and that it can be easily solved numerically using dynamic programming. Numerical study (not shown) suggests that the relative contribution of q2(a + cA, b + cB, 0) to q(a, b, c) is in most cases extremely small. (See Jenkins and Song (2009) for details.) Deriving an analytic expression for σ(a, b, c) in (3.6) is a laborious task, as the long equation in Appendix A suggests. We have written a computer program to verify numerically that our analytic result is correct.

4. Joint distribution of the number of alleles at the two loci in a sample

Following the same strategy as in the previous section, we can obtain the asymptotic behavior of the joint distribution of the number of alleles observed at the two loci in a sample. To make explicit the dependence of these numbers on the sample size, write the number of alleles at locus A as Ka,b,c and the number of alleles at locus B as La,b,c. Ethier and Griffiths (1990) proved that the probability p(a, b, c; k, l):= P(Ka,b,c = k, La,b,c = l) satisfies the recursion

[n(n1)+θA(a+c)+θB(b+c)+ρc]p(a,b,c;k,l)=a(a1+2c)p(a1,b,c;k,l)+b(b1+2c)p(a,b1,c;k,l)+c(c1)p(a,b,c1;k,l)+2abp(a1,b1,c+1;k,l)+θA[ap(a1,b,c;k1,l)+cp(a,b+1,c1;k1,l)]+θB[bp(a,b1,c;k,l1)+cp(a+1,b,c1;k,l1)]+ρcp(a+1,b+1,c1;k,l),
(4.1)

where p(a, b, c; k, l) = 0 if a < 0, b < 0, c < 0, k < 0, l < 0, a = b = c = 0, or k = l = 0. Equation (4.1) has a unique solution satisfying the initial conditions

p(1,0,0;k,l)=δk,1δl,0,p(0,1,0;k,l)=δk,0δl,1,

for k, l = 0, 1, …, n.

As with Golding’s recursion, equation (4.1) can be solved numerically, but quickly becomes computationally intractable with growing n. The only exception is the special case of ρ = ∞, for which the distribution is given by the product of (2.4) for each locus. In what follows, we use the following notation in writing an asymptotic series for p(a, b, c; k, l):

Definition 4.1

For loci A and B, respectively, we define the analogues of (2.4) as

pA(a;k)=s(a,k)θAk(θA)a,
(4.2)

and

pB(b;l)=s(b,l)θBl(θB)b,
(4.3)

where s(a, k) and s(b, l) are the Stirling numbers of the first kind.

We pose the expansion

p(a,b,c;k,l)=p0(a,b,c;k,l)+p1(a,b,c;k,l)ρ+O(1ρ2),
(4.4)

for large ρ. Then, in Section 5.5 we prove the following result for the zeroth order term:

Proposition 4.1

For an asymptotic expansion of the form (4.4) satisfying the recursion (4.1), p0(a, b, c; k, l) is given by

p0(a,b,c;k,l)=pA(a+c;k)pB(b+c;l).
(4.5)

Similar to Lemma 3.1, we obtain the following vanishing result for the first order term in the case of c = 0:

Lemma 4.1

For an asymptotic expansion of the form (4.4) satisfying the recursion (4.1), we have

p1(a,b,0;k,l)=0.

Using this lemma, it is then possible to obtain the following result for an arbitrary c:

Proposition 4.2

For an asymptotic expansion of the form (4.4) satisfying the recursion (4.1), p1(a, b, c; k, l) is given by

p1(a,b,c;k,l)=c(c1)2[pA(a+c;k)pA(a+c1;k)]×[pB(b+c;l)pB(b+c1;l)].
(4.6)

Proofs of Proposition 4.2 and Lemma 4.1 are provided in Section 5.6 and Section 5.7, respectively.

5. Proofs of main results

In what follows, we provide proofs of the results mentioned in the previous two sections.

5.1. Proof of Proposition 3.1

First assume c > 0. Substitute the expansion (3.1) into Golding’s recursion (2.5), divide by ρc and let ρ → ∞. We are then left with

q0(a,b,c)=i=1Kj=1Lcijcq0(a+ei,b+ej,ceij).
(5.1)

Now, applying (5.1) repeatedly gives

q0(a,b,c)=orderingsΠ(i,j)[K]×[L]cij!c!q0(a+cA,b+cB,0),

where the summation is over all distinct orderings of the c gametes with multiplicity c = (cij). There are c!Π(i,j)cij! such orderings and since the summand is independent of the ordering, we conclude

q0(a,b,c)=q0(a+cA,b+cB,0).
(5.2)

Clearly, (5.2) also holds for c = 0. From a coalescent perspective, this equation tells us that any gamete with specified alleles (i.e., “carrying ancestral material”) at both loci must undergo recombination instantaneously backwards in time.

Now, by substituting the asymptotic expansion (3.1) with c = 0 into Golding’s recursion (2.5) and letting ρ → ∞, we obtain

[n(n1)+θAa+θBb]q0(a,b,0)=i=1Kai(ai1)q0(aei,b,0)+j=1Lbj(bj1)q0(a,bej,0)+2i=1Kj=1Laibjq0(aei,bej,eij)+θAi=1Kδai,1q0(aei,b,0)+θBj=1Lδbj,1q0(a,bej,0).
(5.3)

Equation (5.2) implies q0(aei, bej, eij) = q0(a, b, 0), so with a bit of rearranging we are left with

[a(a+θA1)+b(b+θB1)]q0(a,b,0)=i=1Kai(ai1)q0(aei,b,0)+j=1Lbj(bj1)q0(a,bej,0)+θAi=1Kδai,1q0(aei,b,0)+θBj=1Lδbj,1q0(a,bej,0).
(5.4)

with boundary conditions q0(ei, 0, 0) = q0(0, ej, 0) = 1 for all i [set membership] [K] and j [set membership] [L]. Noting that (5.4) is the sum of two independent recursions of the form (2.3), one for each locus and each with appropriate boundary condition, we conclude that q0(a, b, 0) is given by

q0(a,b,0)=qA(a)qB(b),
(5.5)

a product of two (ordered) ESFs. It is straightforward to verify that (5.5) satisfies (5.4). Finally, using (5.2) and (5.5), we arrive at (3.4).

5.2. Proof of Theorem 3.1

First assume c > 0. Substitute the asymptotic expansion (3.1) into Golding’s recursion (2.5), eliminate terms of order ρ by applying (5.1), and let ρ. After applying (5.2) to the remaining terms and invoking (5.4), with some rearrangement we obtain

cq1(a,b,c)i=1Kj=1Lcijq1(a+ei,b+ej,ceij)=c(c1)q0(a+cA,b+cB,0)i=1Kci·(ci·1)q0(a+cAei,b+cB,0)j=1Lc·j(c·j1)q0(a+cA,b+cBej,0)+i=1Kj=1Lcij(cij1)q0(a+cAei,b+cBej,0).

Now, by utilizing (5.5), this can be written in the form

q1(a,b,c)=f(a,b,c)+i=1Kj=1Lcijcq1(a+ei,b+ej,ceij),
(5.6)

where

f(a,b,c):=(c1)qA(a+cA)qB(b+cB)qB(b+cB)i=1Kci·(ci·1)cqA(a+cAei)qA(a+cA)j=1Lc·j(c·j1)cqB(b+cBej)+i=1Kj=1Lcij(cij1)cqA(a+cAei)qB(b+cBej).
(5.7)

Above we assumed c > 0. We define f(a, b, c) = 0 if c = 0. Iterating the recursion (5.6), we may write q1(a, b, c) as

q1(a,b,c)=f(a,b,c)+i=1Kj=1Lcijc[f(a+ei,b+ej,ceij)+i=1Kj=1Lcijδiiδjjc1q1(a+ei+ei,b+ej+ej,ceijeij)].

Similarly, repeatedly iterating (5.6) yields

q1(a,b,c)=q1(a+cA,b+cB,0)+f(a,b,c)+i1j1ci1j1cf(a+ei1,b+ej1,cei1j1)+i1j1,i2j2ci1j1cci2j2δi1j1,i2j2c1×f(a+ei1+ei2,b+ej1+ej2,cei1j1ei2j2)++i1j1,,icjcΠijcij!c!f(a+cA,b+cB,0).
(5.8)

The key observation is that the right-hand side of (5.8) has a nice probabilistic interpretation which allows us to obtain a closed-form formula. To be more precise, consider the first summation

i1j1ci1j1cf(a+ei1,b+ej1,cei1j1).

For a fixed sample configuration c, this can be interpreted as the sum over all possible ways of throwing away a gamete at random and calculating f based on the remaining subsample, which we will denote c(c−1). Equivalently, it is the expected value of f with respect to subsampling without replacement c − 1 of the gametes in c. Write this as

E[f(A(c1),B(c1),C(c1))],

where C(c−1) is the random subsample obtained by sampling without replacement c − 1 gametes from c, and A(c1):=a+cACA(c1),B(c1):=b+cBCB(c1). Note that once the subsample c(c−1) is obtained, then a(c−1) and b(c−1) are fully specified. More generally, consider the (cm)th sum in (5.8). A particular term in the summation corresponds to an ordering of cm gametes in c, which, when removed leave a subsample c(m). With respect to this subsample, the summand is

i=1Kj=1Lcij!cij(m)!m!c!f(a(m),b(m),c(m)),

and for each such subsample c(m) there are (cmcc(m)) distinct orderings of the remaining types in c, with each ordering contributing the same amount to the sum. Here, (cmcc(m)) denotes the multinomial coefficient:

(cmcc(m))=(cm)!i=1Kj=1L(cijcij(m))!.

Gathering identical terms, the (cm)th sum in (5.8) can therefore be written over all distinct subsamples of c of size m:

c(m)(cmcc(m))i=1Kj=1Lcij!cij(m)!m!c!f(a(m),b(m),c(m))=c(m)1(cm)i=1Kj=1L(cijcij(m))f(a(m),b(m),c(m))=E[f(A(m),B(m),C(m))],

where, for a fixed m, C(m)=(Cij(m)) is a multivariate hypergeometric(c, c, m) random variable; that is,

P((i,j)[K]×[L][Cij(m)=cij(m)])=1(cm)(i,j)[K]×[L](cijcij(m)).

Furthermore, marginally we have

Cij(m)hypergeometric(c,cij,m),Ci·(m)hypergeometric(c,ci·,m),C·j(m)hypergeometric(c,c·j,m).

In summary, (5.8) can be written as

q1(a,b,c)=q1(a+cA,b+cB,0)+m=1cE[f(A(m),B(m),C(m))].
(5.9)

According to Lemma 3.1, the first term q1(a + cA, b + cB, 0) vanishes, so we are left with

q1(a,b,c)=m=1cE[f(A(m),B(m),C(m))].
(5.10)

Finally, since A(m)+CA(m)=a+cA and B(m)+CB(m)=b+cB, (5.7) and (5.10) together imply

q1(a,b,c)=m=1c[(m1)qA(a+cA)qB(b+cB)qB(b+cB)1mi=1KE[Ci·(m)(Ci·(m)1)]qA(a+cAei)qA(a+cA)1mj=1LE[C·j(m)(C·j(m)1)]qB(b+cBej)+1mi=1Kj=1LE[Cij(m)(Cij(m)1)]qA(a+cAei)qB(b+cBej)].

The moments in this equation are easy to compute and one can sum them over m to obtain the desired result (3.5).

5.3. Proof of Lemma 3.1

First note that for any sample (a, b, c) and any subsample of the form (a(1), b(1), c(1)), we have f(a(1), b(1), c(1)) = 0, since every term on right-hand side of (5.7) has a vanishing coefficient. So, equation (5.9) implies

q1(aei,bej,eij)=q1(a,b,0),
(5.11)

for any (i, j) [set membership] [K] × [L]. Now, substitute the asymptotic expansion (3.1) with c = 0 into Golding’s recursion (2.5). Note that terms of order ρ are absent since c = 0. Eliminate terms with coefficients independent of ρ by applying (5.3), multiply both sides of the recursion by ρ, and let ρ → ∞ to obtain the following:

[n(n1)+θAa+θBb]q1(a,b,0)=i=1Kai(ai1)q1(aei,b,0)+j=1Lbj(bj1)q1(a,bej,0)+2i=1Kj=1Laibjq1(aei,bej,eij)+θAi=1Kδai,1q1(aei,b,0)+θBj=1Lδbj,1q1(a,bej,0),

with boundary conditions q1(ei, 0, 0) = q1(0, ej, 0) = 0 for all i [set membership] [K] and j [set membership] [L]. This equation can be made strictly recursive by applying (5.11) to q1(aei, bej, eij). It therefore follows from the boundary conditions (for example by induction) that q1(a, b, 0) = 0.

5.4. Proof of Theorem 3.2

Here, we provide only an outline of a proof; details are similar to the proof of Theorem 3.1. Substitute the asymptotic expansion (3.1) into Golding’s recursion (2.5), eliminate terms with coefficients proportional to ρ or independent of ρ. Then, multiply both sides of the recursion by ρ and let ρ → ∞ to obtain

cq2(a,b,c)i=1Kj=1Lcijq2(a+ei,b+ej,ceij)=i=1Kai(ai1+2ci·)q1(aei,b,c)+j=1Lbj(bj1+2c·j)q1(a,bej,c)+i=1Kj=1L[cij(cij1)q1(a,b,ceij)+2aibjq1(aei,bej,c+eij)]+θAi=1K[j=1Lδai+ci·,1δcij,1q1(a,b+ej,ceij)+δai,1δci·,0q1(aei,b,c)]+θBj=1L[i=1Kδbj+c·j,1δcij,1q1(a+ei,b,ceij)+δbj,1δc·j,0q1(a,bej,c)][n(n1)+θA(a+c)+θB(b+c)]q1(a,b,c)
(5.12)

By substituting our expression (3.5) for q1(a, b, c), the right-hand side can be expressed as a function g(a, b, c) which is completely known but rather cumbersome to write down. As in the proof of Theorem 3.1, the same ‘unwrapping’ maneuver can be applied to rearrange (5.12) into the form (3.6), where

σ(a,b,c)=m=1cE[g(A(m),B(m),C(m))].

This time An external file that holds a picture, illustration, etc.
Object name is nihms152909ig2.jpg[g(A(m), B(m), C(m))] is a function of fourth-order moments of the multivariate hypergeometric distribution. The formula shown in Appendix A is obtained by evaluating the expectations and summing over m.

We now show that q2(a, b, 0) satisfies the recursion shown in (3.7). We will use the fact that for a sample (aei, bej, eij), we have

g(aei,bej,eij)=2(a1)(b1)qA(a)qB(b)2(b1)(ai1)qA(aei)qB(b)2(a1)(bj1)qA(a)qB(bej)+2(ai1)(bj1)qA(aei)qB(bej).
(5.13)

For an arbitrary c, g(a, b, c) is much more complicated.

Now, one can adopt the approach used in the proof of Lemma 3.1 to obtain a strict recursion for q2(a + cA, b + cB, 0). First, note that (3.6) and (5.13) imply

q2(aei,bej,eij)=q2(a,b,0)+E[g((Aei)(1),(Bej)(1),eij(1))]=q2(a,b,0)+g(aei,bej,eij)=q2(a,b,0)+2(a1)(b1)qA(a)qB(b)2(b1)(ai1)qA(aei)qB(b)2(a1)(bj1)qA(a)qB(bej)+2(ai1)(bj1)qA(aei)qB(bej).
(5.14)

As before, substitute the asymptotic expansion (3.1) for c = 0 into Golding’s recursion (2.5), eliminate terms with coefficients independent of ρ or proportional to ρ−1, and let ρ → ∞ to obtain

[n(n1)+θAa+θBb]q2(a,b,0)=i=1Kai(ai1)q2(aei,b,0)+j=1Lbj(bj1)q2(a,bej,0)+2i=1Kj=1Laibjq2(aei,bej,eij)+θAi=1Kδai,1q2(aei,b,0)+θBj=1Lδbj,1q2(a,bej,0),

with boundary conditions q2(ei, 0, 0) = q2(0, ej, 0) = 0 for all i [set membership] [K] and j [set membership] [L]. This equation can be made strictly recursive by applying (5.14) to q2(aei, bej, eij). After some simplification, this leads to the recursion (3.7).

5.5. Proof of Proposition 4.1

The proof is similar to the proof of Proposition 3.1, working with the system (4.1) rather than Golding’s recursion (2.5). First assume c > 0. Substitute the expansion (4.4) into the recursion (4.1), divide by ρc and let ρ. We are left with

p0(a,b,c;k,l)=p0(a+1,b+1,c1;k,l),

which implies

p0(a,b,c;k,l)=p0(a+c,b+c,0;k,l).
(5.15)

Clearly, (5.15) also holds for c = 0.

Now, by substituting the asymptotic expansion (4.4) with c = 0 into (4.1) and letting ρ, we obtain

[n(n1)+θAa+θBb]p0(a,b,0;k,l)=a(a1)p0(a1,b,0;k,l)+b(b1)p0(a,b1,c;k,l)+2abp0(a1,b1,1;k,l)+θAap0(a1,b,0;k1,l)+θBbp0(a,b1,0;k,l1).
(5.16)

After invoking (5.15) on p0(a − 1, b − 1, 1; k, l) and rearranging, we are left with

[a(a+θA1)+b(b+θB1)]p0(a,b,0;k,l)=a(a1)p0(a1,b,0;k,l)+b(b1)p0(a,b1,0;k,l)+θAap0(a1,b,0;k1,l)+θBbp0(a,b1,0;k,l1),
(5.17)

with boundary conditions p0(1, 0, 0; k, l) = δk,1δl,0, and p0(0, 1, 0; k, l) = δk,0δl,1. Equation (5.17) can be expressed as a linear sum of two independent recursions:

(a1+θA)p0A(a;k)=(a1)p0A(a1;k)+θAp0A(a1;k1),(b1+θB)p0B(b;l)=(b1)p0B(b1;l)+θBp0B(b1;l1),

with respective boundary conditions p0A(1;k)=δk,1 and p0B(1;l)=δl,1. These recursions are precisely those considered by Ewens (1972, eq. 21), with respective solutions (4.2) and (4.3). Hence, p0A(a;k)=pA(a;k) and p0B(b;l)=pB(b;l), and it is straightforward to verify that pA(a; k)pB(b; l) satisfies (5.17). Substituting this solution into (5.15), we arrive at (4.5), as required.

5.6. Proof of Proposition 4.2

First assume c > 0. Substitute the asymptotic expansion (4.4) into the recursion (4.1), eliminate terms with coefficients linear in ρ by applying (5.15), and let ρ. After applying (5.15) to the remaining terms and invoking (5.17), with some rearrangement we obtain

p1(a,b,c;k,l)p1(a+1,b+1,c1;k,l)=(c1)[p0(a+c,b+c,0;k,l)p0(a+c1,b+c,0;k,l)p0(a+c,b+c1,0;k,l)+p0(a+c1,b+c1,0;k,l)].
(5.18)

Applying the recursion repeatedly, this becomes

p1(a,b,c;k,l)=p1(a+c,b+c,0;k,l)+[p0(a+c,b+c,0;k,l)p0(a+c1,b+c,0;k,l)p0(a+c,b+c1,0;k,l)+p0(a+c1,b+c1,0;k,l)]m=0c1(c1m).
(5.19)

According to Lemma 4.1, the first term p1(a + c, b + c, 0; k, l) vanishes. Hence, since p0(a, b, c; k, l) is given by (4.5), the right-hand side of (5.19) is fully known. With some rearrangement we are left with (4.6).

5.7. Proof of Lemma 4.1

First note that (5.18) implies

p1(a1,b1,1;k,l)=pl(a,b,0;k,l).
(5.20)

Now, substitute the asymptotic expansion (4.4) with c = 0 into (4.1), eliminate leading order terms by applying (5.16), and let ρ. The result is made strictly recursive by invoking (5.20), and we obtain

[a(a+θA1)+b(b+θB1)]p1(a,b,0;k,l)=a(a1)p1(a1,b,0;k,l)+b(b1)p1(a,b1,0;k,l)+θAap1(a1,b,0;k1,l)+θBbp1(a,b1,0;k,l1),

with boundary conditions p1(1, 0, 0; k, l) = p1(0, 1, 0; k, l) = 0, for k, l = 0, …, n. It therefore follows (for example by induction) that p1(a, b, 0; k, l) = 0.

Acknowledgments

We thank Robert C. Griffiths, Charles H. Langley, and Joshua Paul for useful discussions.

APPENDIX A: EXPRESSION FOR σ(a, b, c)

We use QA to denote qA(a + cA), QiA to denote qA(a + cAei), QikA to denote qA(a + cAeiek), and so on. Then,

σ(a,b,c)=c3[(c1)(c+1)(3c2)8+(c1)(3a+3b+2c1)+6ab]QAQBθA(c1)2i=1Kδai,0δci·,1QiAQBθB(c1)2j=1Lδbj,0δc·j,1QAQjB+i=1K[θAc(c3)+2a+4b44ci·(ci·1)(2b+c1)ci·(ai+ci·1)]QiAQB+12i=1K[θA2δci·,2+56ai4ci·6]ci·(ci·1)QiiAQB+j=1L[θBc(c3)+2b+4a44c·j(c·j1)(2a+c1)c·j(bj+c·j1)]QAQjB+12j=1L[θB2δc·j,2+56bj4c·j6]c·j(c·j1)QAQjjB+i,k=1Kci·(ci·1)ck·(ck·1)8QikAQB+j,l=1Lc·j(c·j1)c·l(c·l1)8QAQjlBθA+θBc(c5)+2a+2b44i=1Kj=1Lcij(cij1)QiAQjB+i=1Kj=1L[ci·(ci·1)c·j(c·j1)4+cij(cij+12ci·+2ci·c·j2c·j)2+cijbj(ci·1)+cijai(c·j1)+2aibjcij+θB2δbj,0δc·j,1δcij,1(ci·1)+θA2δai,0δci·,1δcij,1(c·j1)]QiAQjB+12i=1K[ai+ci·1θA2δci·,2]j=1Lcij(cij1)QiiAQjB+12j=1L[bj+c·j1θB2δc·j,2]i=1Kcij(cij1)QiAQjjB14i=1Kj=1Lk=1Kcij(cij1)ck·(ck·1)QikAQjB14i=1Kj=1Ll=1Lcij(cij1)c·l(c·l1)QiAQjlB+18i=1Kj=1Lk=1Kl=1Lcij(cij1)ckl(ckl1)QikAQjlB112i=1Kj=1Lcij(cij1)(2cij1)QiiAQjjB.

To check the correctness of the above expression, we also solved the recursion (5.12) numerically for all sample configurations of sizes n = 10, 20, and 30 (with K, L ≤ 2), and confirmed that the above analytic expression agreed in all cases. We also implemented a Mathematica program to solve q(a, b, c) exactly in the special case K = L = 1. The program can return series expansions in terms of ρ−1 as ρ → ∞ which are symbolic in θA and θB. We could then compare the first three terms against q0, q1, and q2, for various sample configurations (a, b, c).

Footnotes

AMS 2000 subject classifications: Primary 92D15; secondary 65C50,92D10

Contributor Information

Paul A. Jenkins, Computer Science Division, University of California, Berkeley, Berkeley, CA 94720, USA.

Yun S. Song, Department of Statistics and Computer Science Division, University of California, Berkeley, Berkeley, CA 94720, USA.

References

  • Arratia A, Barbour AD, Tavaré S. Logarithmic Combinatorial Structures: A Probabilistic Approach. European Mathematical Society Publishing House; Switzerland: 2003.
  • De Iorio M, Griffiths RC. Importance sampling on coalescent histories I. Adv Appl Prob. 2004a;36:417–433.
  • De Iorio M, Griffiths RC. Importance sampling on coalescent histories II. Adv Appl Prob. 2004b;36:434–454.
  • Ethier SN, Griffiths RC. On the two-locus sampling distribution. J Math Biol. 1990;29:131–159.
  • Ewens WJ. The sampling theory of selectively neutral alleles. Theor Popul Biol. 1972;3:87–112. [PubMed]
  • Fearnhead P, Donnelly P. Estimating recombination rates from population genetic data. Genetics. 2001;159:1299–1318. [PubMed]
  • Golding GB. The sampling distribution of linkage disequilibrium. Genetics. 1984;108:257–274. [PubMed]
  • Griffiths RC. Neutral two-locus multiple allele models with recombination. Theor Popul Biol. 1981;19:169–186.
  • Griffiths RC. The two-locus ancestral graph. In: Basawa IV, Taylor RL, editors. Selected proceedings of the Shefield symposium on applied probability: 18. IMS Lecture Notes—Monograph series. Vol. 18. 1991. pp. 100–117.
  • Griffiths RC, Lessard S. Ewens’ sampling formula and related formulae: combinatorial proofs, extensions to variable population size and applications to ages of alleles. Theor Popul Biol. 2005;68:167–77. [PubMed]
  • Griffiths RC, Marjoram P. Ancestral inference from samples of DNA sequences with recombination. J Comput Biol. 1996;3:479–502. [PubMed]
  • Griffiths RC, Jenkins PA, Song YS. Importance sampling and the two-locus model with subdivided population structure. Adv Appl Prob. 2008;40:473–500. [PMC free article] [PubMed]
  • Hoppe F. Pólya-like urns and the Ewens’ sampling formula. J Math Biol. 1984;20:91–94.
  • Hudson RR. The sampling distribution of linkage disequilibrium under an infinite allele model without selection. Genetics. 1985;109:611–631. [PubMed]
  • Hudson RR. Two-locus sampling distributions and their application. Genetics. 2001;159:1805–1817. [PubMed]
  • Jenkins PA, Song YS. Closed-form two-locus sampling distributions: accuracy and universality. Genetics. 2009 In press. [PubMed]
  • Kingman JFC. On the genealogy of large populations. J Appl Prob. 1982a;19:27–43.
  • Kingman JFC. The coalescent. Stoch Process Appl. 1982b;13:235–248.
  • Kuhner MK, Yamato J, Felsenstein J. Maximum likelihood estimation of recombination rates from population data. Genetics. 2000;156:1393–1401. [PubMed]
  • McVean GAT, Myers S, Hunt S, Deloukas P, Bentley DR, Donnelly P. The fine-scale structure of recombination rate variation in the human genome. Science. 2004;304:581–584. [PubMed]
  • Myers S, Bottolo L, Freeman C, McVean G, Donnelly P. A fine-scale map of recombination rates and hotspots across the human genome. Science. 2005;310:321–324. [PubMed]
  • Nielsen R. Estimation of population parameters and recombination rates from single nucleotide polymorphisms. Genetics. 2000;154:931–942. [PubMed]
  • Pitman J. Technical Report 345. Dept. Statistics; U.C. Berkeley: 1992. The two-parameter generalization of Ewens’ random partition structure.
  • Pitman J. Exchangeable and partially exchangeable random partitions. Probab Th Rel Fields. 1995;102:145–158.
  • Slatkin M. An exact test for neutrality based on the Ewens sampling distribution. Genet Res. 1994;64:71–4. [PubMed]
  • Stephens M. Inference under the coalescent. In: Balding D, Bishop M, Cannings C, editors. Handbook of Statistical Genetics. chapter 8. Wiley; Chichester, UK: 2001.
  • Stephens M, Donnelly P. Inference in molecular population genetics. J R Statist Soc B. 2000;62:605–655.
  • Wang Y, Rannala B. Bayesian inference of fine-scale recombination rates using population genomic data. Phil Trans R Soc B. 2008;363:3921–3930. [PMC free article] [PubMed]
  • Watterson GA. Heterosis or neutrality? Genetics. 1977;85:789–814. [PubMed]