On the Recombination Rate Estimation in the Presence of Population Substructure.

As recombination events are not uniformly distributed along the human genome, the estimation of fine-scale recombination maps, e.g. HapMap Project, has been one of the major research endeavors over the last couple of years. For simulation studies, these estimates provide realistic reference scenarios to design future study and to develop novel methodology. To achieve a feasible framework for the estimation of such recombination maps, existing methodology uses sample probabilities for a two-locus model with recombination, with recent advances allowing for computationally fast implementations. In this work, we extend the existing theoretical framework for the recombination rate estimation to the presence of population substructure. We show under which assumptions the existing methodology can still be applied. We illustrate our extension of the methodology by an extensive simulation study.

Introduction

The discovery that recombination events in human genome are not uniformly distributed, but concentrated in specific genomic regions, which are typically referred to as recombination hotspots [1], was one of the driving forces behind the HapMap Project [2]. The characterization and understanding of the local linkage disequilibrium structure between genetic variants is fundamental for the discovery of disease susceptibility loci (DSLs). The peak recombination rate within recombination hotspots can be multiple times higher than the recombination rate outside such genomic areas [3]. To create recombination hotspots in simulated data, commonly used software tools based on coalescent simulations, e.g. cosi, incorporate recombination rates that vary along the chromosome [4]. Such software tools, in order to generate realistic data, require an accurate fine-scale recombination rate map as input parameter. McVean et al. developed the LDhat software which estimates the recombination rate along the genome according to a assumed step-wise constant model [5]. The approach is computationally fast and can be applied at a genome-wide level. In 2007, Auton and McVean extended this approach and incorporated recombination hotspots into the assumed form of the recombination rate-model [6]. Both approaches are based on the so-called composite-likelihood approximation [7] that calculates the data likelihood as the product of pairwise two-locus probabilities. In order to obtain these sample probabilities, an Importance Sampling (IS) scheme by Fearnhead and Donnelly [8] was applied which requires exhaustive lookup tables of probabilities for a wide range of recombination rates. These methods were used to estimate the recombination rate for the HapMap samples CEU, YRI and ASN [2] (resp. later the 1000Genomes [9] samples), separately. As explained in [10], the overall genetic map of recombination rates was produced by comparison of the total map lengths to that estimated by the pedigree method in [11] and averaging these maps. In 2009, Jenkins and Song proposed approaches to calculate the sample probability for a given configuration with an analytic asymptotic formula of order two in the reciprocal recombination rate. Their approach was initially intended for an infinite-allele model [12] and later extended to a finite-allele model [13]. In particular, they showed how to transform the symmetric diallelic mutation model in Fearnhead and Donnelly [8] resp. LDhat to a so-called Parent-Independent-Mutation (PIM) model and were therefore able to present a different approach to calculate the required two-locus probabilities for LDhat. The key aspect they showed was that the calculations were independent from the specific value of the recombination rate. Thus, they were able to evaluate the sample probability instantly for every recombination rate. This substantially reduces the computational burden. In addition, for large recombination rate values, Monte Carlo based methods as the IS scheme by Fearnhead and Donnelly are less efficient, since the number of recombination events increases and the sampled genealogies become very complicated, whereas the asymptotic sampling formula by Jenkins and Song becomes very precise for large values. In 2012, Jenkins and Song derived an arbitrary order asymptotic expansion for the finite-allele model, using the generator of the corresponding Wright-Fisher diffusion [14]. Together with the application of Pade approximations, this provided a way to evaluate two-locus probabilities efficiently for a wide range of recombination rates. The mentioned diffusion generator used in [14] describes the diffusion limit of the Markov Chain explained in [15]. In this communication, we extend this underlying Markov Chain model to the presence of population substructure with a finite number of subpopulations. Under the assumption of strong migration between the subpopulations, we derive the corresponding diffusion limit of a specific weighted mean of subpopulation frequencies. The key result of our work is, that differences between subpopulation frequencies disappear and the diffusion limit has the form as in the panmictic case with rescaled effective population size. This implies the possibility to combine subpopulation samples and evaluate the corresponding sample probability with the existing methodology without any additional computation effort. The advantage for the estimation of recombination rates lies in the increased sample size and the more realistic underlying model.

Methods

We start with a summary of the results in [13]. We denote the sample, i.e. the genetic data at locus A and B for n study subjects, by s. The parameter ρ = 4N r is the population-scaled recombination rate. The parameters θ and θ refer to the mutation rates at locus A and B. Let r be the number of possible alleles at the first locus A and r be the number of possible alleles at the other locus B. Denote the different alleles by A 1, …, A and B 1, …, B . For notational convenience, we introduce [k]: = {1, …, k}. First, we summarize the work in [14] to emphasize the differences to our model later. Write the sample configuration s for a single population by with where c gives the number of gametes with allele A at locus A and B at locus B. Similar, where a gives the number of gametes with allele A at locus A and unspecified allele at locus B. b is defined analogously. Further, write , , and n = a + b + c. This notation is also used in [13]. Jenkins and Song utilize the fact that, under a stationary distribution, the expectation of a suitable function, which is applied to the generator (described in [14]), is equal to zero. If q(s; ρ) describes the sample probability with reference to ρ, they derived that the sample probabilities can be calculated from the following linear system with boundary conditions where π and π denote the stationary distributions of the mutation models for locus A and B. Since this linear system has to be solved for every recombination rate ρ separately and grows rapidly in n, they propose an asymptotic expansion in

Now, the corresponding recursions are solved only once for a fixed mutation model and the sample probabilities, with reference to several recombination rates, can be evaluated by plugging in the recombination parameter ρ. The estimation procedures in [5] and [6] are based on two-locus probabilities from the model in [8] with respect to some mutation parameter θ , recombination rate ρ and a symmetric diallelic mutation model

As shown in [13], both models are in line if we set r = r = 2, θ = 2θ and

As explained in the introduction, the generator in [14] corresponds to the diffusion limit of the discrete Markov Chain in [15]. In the next section, we introduce the population substructured extension of this discrete Markov Chain, derive the diffusion limit of weighted mean frequencies and show that we can utilize the same idea as in [14] to improve the recombination rate estimation framework.

Model setting

As explained, we consider a two-locus model. As in [14], let r be the number of possible alleles at the first locus A and r be the number of possible alleles at the other locus B. The generations are supposed to be non-overlapping and the population is monoecious. We suppose the underlying model for the discrete Wright-Fisher model as it is described in [15], which leads to a diffusion approximation, corresponding to the generator in [14]. But in addition, we assume that the population is subdivided into Γ < ∞ subpopulations. Denote by , α = 1, …, Γ, the fractions of the corresponding subpopulation sizes and let N be the overall population size. Thus, and . Set for α = 1, …, Γ. We model the migration procedure as in the discrete model in [16]. In this context, let m denote the probability that an individual in subpopulation α was in subpopulation β one generation before. These probabilities are given by the backward migration matrix M. The following diagram summarizes the lifecycle of the model (compare with [16]).

The Markov Chain Z , which describes the relative frequencies (P ), is completely defined in S1 Appendix. P denotes the relative frequency of type A B in subpopulation α in the model Z , we suppress the dependency on N for convenience. Note that the subpopulation sizes are assumed to be effectively infinite except shorlty before and after population regulation [16].

Parameter assumptions

In order to obtain a diffusion limit, we make the following assumptions on the model parameters.

Biological parameters

We use the well established assumptions on the parameter scaling for the biological mechanisms. These conditions depend on the effective population size (section 3.7 in [17]). The effective population size for our model is derived below. We write and where u and u are the mutation rates per individual per generation and P as well as P are the transition matrices for mutation for location A resp. B with and as in [14]. See Eq (A.1) and (A.2) in S1 Appendix for the definition of the mutation model. For convenience, we omit technical conditions of the truncation N, since we are only interested in the diffusion limit resp. large N. We assume resp.

For the recombination fractions in Eq (A.3) resp. (A.4) in S1 Appendix, we suppose

Migration

The backward migration matrix M is assumed to be a stochastic matrix, which is irreducible and aperiodic. Migration is independent of time and N (see below) resp. N. This is the strong migration assumption. Denote the stationary distribution of M by ξ. This means that migration dominates all other evolutionary forces. This assumption was used in several papers. We adopted it from [16]. Additionally, this suggestion was analyzed in [18], in other papers of [19-21] and [22-24].

Diffusion approximation

The next step is to analyze the diffusion limit of the Markov Chain Z . As in [15] and [16], this is essentially an application of the diffusion approximation Theorem in [25] in S1 Appendix, which is also stated as Theorem (A.3). We describe by simple calculation that the Theorem is still applicable for our extended Markov Chain S1 Appendix. Important for the derivation of the diffusion limit is the connection between N and N , the effective population size. Nagylaki observed and described this connection in the setting of a subdivided population with strong migration in a one-locus model in [16]. He stated where

Therefore, N ≤ N and in particular N is equal to N if and only if . This scenario is called conservative migration. These results are still valid in our two-locus scenario. Analogous to [16], we define where and where d = P − P . In addition,

The idea is to set the two timescales in Theorem (A.3) to ε = 2N and δ = 1. Recombination and mutation work on the slow timescale ε , but the strong migration assumption implies that migration works on the fast timescale δ . Then, it follows that the scaled X ([2N •]) converges to a diffusion process X with explicitly known generator (A.17). The process Y describes the derivation of the frequencies within the subpopulations from the weighted mean frequencies. The result is, that Y ([2N t]]) converges to 0, for every t > 0. The interpretation is that, since migration works faster than the biological mechanisms, the differences in relative frequencies within the subpopulations disappear and approach the mean frequencies. In S1 Appendix, it is shown that the process X [2N •] converges weakly to a process X with state space (A.16), which is associated to the generator (A.17) resp. (A.23) with state space (A.22). This generator has the same form as the generator which was used in [14]. The difference lies in the scaled effective population size. This observation is in agreement with the one-locus results in [16]. The assumptions about the mutation model (Eqs (2) and (3)) imply an important fact about the discrete Markov Chains.

Lemma 1: The Markov chain Z has a unique stationary distribution μ .

Proof See S1D Appendix.

Another important observation is that there exists a unique stationary distribution μ for the diffusion process X.

Lemma 2: The diffusion process X has a unique stationary distribution μ.

Proof This is a consequence of the restriction of the results in [26] to a finite number of possible alleles.

Sample configurations

As described in the introduction, the main interest lies in the efficient evaluation of approximate sample probabilities. Therefore, we need to extend the definition of the latter object to incorporate our scenario. A sample configuration for a subdivided population is defined by where Here, we have where c gives the number of gametes in subpopulation α with allele A at locus A and B at locus B. Similar, where a gives the number of gametes in subpopulation α with allele A at locus A and unspecified allele at locus B. b is defined analogously. This is the straightforward extension of the definition in [14].

Sample probabilities

Recall the stationary distributions μ of the discrete model and define with where and P analogously. The probability of a sample configuration for the model can be expressed as

If we have Γ = 1, this is in line with the approach in [14]. Based on the observations above, we can derive that we can approximate these sample probabilities by combined samples for the diffusion limit. More precise:

Lemma 3: where ∪ s is the combined sample for one population with , and for all i, j.

Proof. See S1E Appendix.

μ is the stationary distribution of the diffusion process X and (x ) in the state space K (A.22). We do not distinguish in the notation between these objects in different characterizations, since they can be identified with each other.

Main result

Thus, our main result is that we can combine the samples over all subpopulations and evaluate the corresponding sample probability for the diffusion process, which is stated above. The advantage is that the computational effort is still the same as in the panmictic case. The only adjustment is done by the rescaling of the effective population size. To clarify this connection, note that the generator for a single population in [14], has the same form and depends on the recombination rate ρ and mutation rates θ resp. θ . As mentioned above, Jenkins and Song evaluated the expression for the sample probability, given biological parameters, with a sophisticated recursion technique, which was derived by the fact that

Here denotes the associated generator of the diffusion process. For our scenario, we are dealing with a rescaled effective population size, conditioned on the migration model, by the factor δ in Eq (5) and can apply the same technique.

Results/Simulation study

In order to give a better impression about the results, we present an empirical example for the theory. We consider a model with two subpopulations. A realistic choice for N is about N = 10,000, according to [4]. The fractions of subpopulation sizes are denoted by and . We choose the backward migration matrix to be equal to

Note that this matrix satisfies the migration assumptions and implies . We consider a diallelic model for both loci, i.e. and assume that the mutation rates are the same for loci A and loci B. Furthermore, we choose as in [14]. As it is used in [14] (as explained in the Methods section, this corresponds to the mutation model used in LDhat), reasonable parameter values for θ are θ = θ = θ ∈ [0.001, 0.01] and ρ = 50. To test our results, we simulate the discrete model Z of the relative frequencies in both subpopulations with arbitrary initial distributions according to Eqs (A.4), (A.5), (A.6) and (A.7) in S1 Appendix. For this, we compute the corresponding mutation and recombination rates u 0 and r 0 satisfying and for θ = 0.01, ρ = 50 and N = 10,000, to get realistic parameter values for the Markov Chain. To obtain an estimate of the sample probability under the stationary distribution, we run the model with respect to u 0 and r 0 over 1012 generations and estimate the corresponding quantity for some example configurations every 105 generations. We do this for two different scenarios of subpopulation sizes. The first choice is (q 1, q 2) = ξ and the second . In the first case, the condition for conservative migration is satisfied, δ 1 = 1. Therefore, the theory predicts that the effects of the subdivision disappear. In the second scenario, we compute

Based on these observations, we compare the empirical results for the discrete model with the results from the diffusion approximation. As explained, we have to take care about the corresponding scaled mutation and recombination rates. In the first case, this implies θ = 0.01 and ρ = 50, since δ 1 = 1, and in the second case θ ≈ 0.004235 and ρ ≈ 21.17. The sample probabilities for the diffusion approximation can either be obtained by the software package asf [13] or can be calculated by solving the linear system Eq (1). Solving the linear system gives us the exact value for the sample probabilities, the software asf only calculates a second order approximation. We implement a solver for the linear system, similar to the goldings recursion program from Richard R. Hudson [7], and use these results (the results are consistent with the the sample probabilities from asf). The comparison of the values is given in the Tables 1 and 2. We only consider sample configurations with full gamete information, e.g. a 1 = a 2 = b 1 = b 2 = 0. As we can see in the Tables 1 and 2, the observations are in line with the theoretical predictions.

Table 1

Comparison of estimated and expected probabilities for the scenario (q 1, q 2) = (ξ 1, ξ 2).

sample configuration	discrete	combined sample	diffusion
((0, 0, 1, 0), (0, 0, 0, 1))	0.00123218	(0, 0, 1, 1)	0.00123137
((0, 0, 1, 1), (0, 0, 1, 0))	0.000615191	(0, 0, 2, 1)	0.00061418
((1, 0, 2, 0), (1, 0, 2, 0))	6.1255e-05	(2, 0, 4, 0)	6.1262e-05
((4, 0, 0, 0), (2, 0, 0, 0))	0.245318	(6, 0, 0, 0)	0.244383
((3, 0, 0, 0), (2, 1, 0, 0))	0.000243998	(5, 1, 0, 0)	0.000244123
((1, 1, 0, 0), (0, 0, 0, 1))	1.46358e-06	(1, 1, 0, 1)	1.50463e-06

6 different example sample configurations. The discrete value corresponds to the estimated probability from the discrete model, the diffusion value to exact solution of the linear system.

Table 2

Comparison of estimated and expected probabilities for the scenario .

sample configuration	discrete	combined sample	diffusion
((0, 0, 1, 0), (0, 0, 0, 1))	0.000522292	(0, 0, 1, 1)	0.000525972
((0, 0, 1, 1), (0, 0, 1, 0))	0.000260697	(0, 0, 2, 1)	0.00026271
((1, 0, 2, 0), (1, 0, 2, 0))	2.62979e-05	(2, 0, 4, 0)	2.62440e-05
((4, 0, 0, 0), (2, 0, 0, 0))	0.248091	(6, 0, 0, 0)	0.247599
((3, 0, 0, 0), (2, 1, 0, 0))	0.000104652	(5, 1, 0, 0)	0.000104806
((1, 1, 0, 0), (0, 0, 0, 1))	2.60423e-07	(1, 1, 0, 1)	2.67692e-07

Description: See Table 1.

Discussion

In their work, Jenkins and Song derived an approximate sampling formula of arbitrary order from the generator of the Wright-Fisher diffusion [14]. Relying on first results in [13] and the detailed study of accuracy in [14], we restate that the application of the approximate sampling formula has huge advantages in relation to the computational demanding Monte Carlo based methods, especially for large values of ρ. These theoretical results can be applied to make the existing recombination rate estimation procedures as implemented in LDhat and rhomap more efficient. In this communication, we go one step further and extend the underlying discrete Wright-Fisher model to the scenario in which the sample population consists of distinct subpopulations. We describe the required assumptions for the biological parameters and the migration model, in order to derive the corresponding diffusion limit of the weighted mean frequencies. The strong migration assumption implies that the form of the generator for the diffusion limit is the same as in the panmictic case. Under our mutation and recombination model, we can show the existence and uniqueness of stationary distributions for the discrete Markov Chains and the diffusion process. We describe the connection between these stationary distributions, which is essentially the motivation for the approach to compute approximate sample probabilities. An empirical example shows that our theoretical derivations are valid. The conclusion is that the probability of a sample containing subsamples from each subpopulation can be calculated as the probability of the combined sample with the original approach of [14]. The difference lies in the rescaled effective population size. The advantage is that one can handle 2-loci data for multiple subpopulations without any additional computational burden.

We believe that the suggested extension can help to achieve a more realistic setting for the recombination rate estimation and increased efficiency due to larger sample sizes.

Theoretical derivations and description.

In this file, we describe the underlying Markov Chain, the theoretical tools behind the diffusion limit and the proofs of Lemma 1 and 3.

(PDF)

Click here for additional data file.