A principled approach to deriving approximate conditional sampling distributions in population genetics models with recombination

Abstract

The multilocus conditional sampling distribution (CSD) describes the probability that an additionally sampled DNA sequence is of a certain type, given that a collection of sequences has already been observed. The CSD has a wide range of applications in both computational biology and population genomics analysis, including phasing genotype data into haplotype data, imputing missing data, estimating recombination rates, inferring local ancestry in admixed populations, and importance sampling of coalescent genealogies. Unfortunately, the true CSD under the coalescent with recombination is not known, so approximations, formulated as hidden Markov models, have been proposed in the past. These approximations have led to a number of useful statistical tools, but it is important to recognize that they were not derived from, though were certainly motivated by, principles underlying the coalescent process. The goal of this article is to develop a principled approach to derive improved CSDs directly from the underlying population genetics model. Our approach is based on the diffusion process approximation and the resulting mathematical expressions admit intuitive genealogical interpretations, which we utilize to introduce further approximations and make our method scalable in the number of loci. The general algorithm presented here applies to an arbitrary number of loci and an arbitrary finite-alleles recurrent mutation model. Empirical results are provided to demonstrate that our new CSDs are in general substantially more accurate than previously proposed approximations.

Figure 1.—

Illustrations of a genealogy and conditional genealogy for a two-locus (k = 2), two-allele model. The two loci of a haplotype are each represented by a circle, with the shading (light or dark) indicating the allelic type at that locus. Mutation events, along with the locus and resulting haplotype, are indicated by small arrows. Recombination events (always taking the left loci from the left side and the right locus from the right side), along with the resulting haplotype, are indicated by dotted circles. (a) A genealogy with n = 4. It is easy to verify that, starting with the MRCA and following the genealogy forward in time, the sample configuration n shown at the leaves is obtained. (b) An “observed” genealogy with n = 3 and a conditional genealogy with m = 1. Absorption events are indicated by dotted arrows into . Following the combined genealogy forward in time, it is easy to check that the conditional sample m shown at the leaf of is obtained.

Figure 2.—

Illustration of a conditional genealogy using the approximation . Absorption events are indicated by dotted arrows into the “trunk” ancestry . Comparing with Figure 1b, observe that is time invariant and extends infinitely into the past.

Figure 3.—

Relative error of CSDs for θ₀ = 1 and ρ₀ = 4. See (14) for definition of . With θ₀ = 1 and ρ₀ = 4, we used a coalescent simulator to generate 250 data sets, each with 25 haplotypes and 10 loci. Then, requisite k-locus, n-haplotype conditional configurations {C⁽ⁱ⁾}_{i=1, …, 250} were obtained using method M1 described in the text. (a) k ∈ {2, 3, 4, 5, 6, 8, 10}, n = 6, and ρ = ρ₀. (b) k = 4, n = 6, and ρ ∈ {0, 2, 4, 6, 8, 12, 16, 20}. (c) k = 4, n ∈ {2, 4, 6, 8, 10, 14, 20}, and ρ = ρ₀.

Figure 4.—

Relative error of CSDs for θ₀ = 0.01 and ρ₀ = 0.1. See (14) for definition of . With θ₀ = 0.01 and ρ₀ = 0.1, we used a coalescent simulator to generate 250 data sets, each with 25 haplotypes and 500 loci. Then, requisite k-locus n-haplotype conditional configurations {C⁽ⁱ⁾}_{i=1, …, 250} were obtained using method M2 described in the text. (a) k ∈ {2, 3, 4, 5, 6, 8, 10}, n = 6, and ρ = ρ₀. (b) k = 4, n = 6, and ρ ∈ {0, 4, 8, 12, 16, 20, 30, 40, 50} × 10⁻². (c) k = 4, n ∈ {2, 4, 6, 8, 10, 14, 20}, and ρ = ρ₀.

Figure 5.—

Relative error of PAC likelihoods for θ₀ = 1 and ρ₀ = 4. See (15) for definition of . With θ₀ = 1 and ρ₀ = 4, we used a coalescent simulator to generate 250 data sets, each with 25 haplotypes and 10 loci. Then, requisite k-locus n-haplotype configurations {n⁽ⁱ⁾}_{i=1, …, 250} were obtained using method M1 described in the text. (a) k = 3, n = 25, and ρ ∈ {0, 2, 4, 6, 8, 12, 16, 20}. (b) k = 5, n = 25, and ρ ∈ {0, 2, 4, 6, 8, 12, 16, 20}.

Figure 6.—

Relative error of PAC likelihoods for θ₀ = 0.01 and ρ₀ = 0.1. See (15) for definition of . With θ₀ = 0.01 and ρ₀ = 0.1, we used a coalescent simulator to generate 250 data sets, each with 25 haplotypes and 500 loci. Then, requisite k-locus n-haplotype configurations {n⁽ⁱ⁾}_{i=1, …, 250} were obtained using method M2 described in the text. (a) k = 3, n = 25, and ρ ∈ {0, 4, 8, 12, 16, 20, 30, 40, 50} × 10⁻². (b) k = 5, n = 25, and ρ ∈ {0, 4, 8, 12, 16, 20, 30, 40, 50} × 10⁻².

Figure 7.—

Approximate values of signed PACErr for θ₀ = 0.01 and ρ₀ = 0.1, corresponding to Figure 6b. The correspondence between the symbols and 's is the same as in previous figures.