Bayesian Inference of Selection in the Wright-Fisher Diffusion Model

11/07/2017 ∙ by Jeffrey J. Gory, et al. ∙ 0

The increasing availability of population-level allele frequency data across one or more related populations necessitates the development of methods that can efficiently estimate population genetics parameters, such as the strength of selection acting on the population(s), from such data. Existing methods for this problem in the setting of the Wright-Fisher diffusion model are primarily likelihood-based, and rely on numerical approximation for likelihood computation and on bootstrapping for assessment of variability in the resulting estimates, requiring extensive computation. Recent work (Jenkins and Spano, 2015) has provided a method for obtaining exact samples from general Wright-Fisher diffusion processes, enabling the development of methods for Bayesian estimation in this setting. We develop and implement a Bayesian method for estimating the strength of selection based on the Wright-Fisher diffusion for data sampled at a single time point. The method utilizes the work of Jenkins and Spano (2015) to develop a Markov chain Monte Carlo algorithm to draw samples from the joint posterior distribution of the selection coefficient and the allele frequencies. We demonstrate that when assumptions about the initial allele frequencies are accurate the method performs well for both simulated data and for an empirical data set on hypoxia in flies (Zhou et al. 2011), where we find evidence for strong positive selection in a region of chromosome 2L previously identified by Ronen et al. (2013). We discuss possible extensions of our method to the more general settings commonly encountered in practice, highlighting the advantages of Bayesian approaches to inference in this setting.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The Wright-Fisher model (Fisher, 1930; Wright, 1931) is the most widely used model for the evolution of genetic variation at a locus over time. In its most basic form, the Wright-Fisher model assumes a population of individuals, each with two possible genetic variants (called alleles) at the locus of interest, for a total of alleles in the population. Of those alleles, are assumed to be of type , and the remaining alleles are assumed to be of type . The Wright-Fisher model traces the proportion of alleles across generations under the assumption of discrete, non-overlapping generations of constant size, in which the next generation of alleles is produced by sampling from the current generation at random with replacement. Under this model, it is easy to see that the number of alleles at time , denoted

, is a discrete-time Markov process with transition probabilities given by the Binomial distribution with

trials and probability of success. In this model, states and are absorbing states.

Application of the basic Wright-Fisher model is limited, however, in several ways. For example, it deals only with variation at a single locus in a single population, and it neglects other important aspects of the evolutionary process, such as selection for advantageous or against deleterious mutations, changes in population size over time, and the generation of new allelic variants through spontaneous mutation. While these processes may be straight-forwardly modeled, inference in the discrete-time framework in these more complex situations is generally not computationally tractable.

For this reason, two common continuous-time approximations are typically employed. The first is generally referred to as Kingman’s coalescent process (Kingman, 1982a, b, c), and arises by examining the limiting distribution of the time at which a collection of sampled alleles last shared a common ancestor as the population size becomes large relative to the sample size. The coalescent approach has proven extremely useful for species-level phylogenetic inference for moderate numbers of species or populations (generally, three to 20 species/populations) under the basic Wright-Fisher model in which the effects of selection are negligible and in which no gene flow between distinct species occurs (Edwards, 2009; Knowles and Kubatko, 2010). However, inference in the presence of selection and/or migration between species or populations is much more difficult using Kingman’s coalescent (Wakeley and Hey, 1997; Nielsen and Wakeley, 2001; Hey and Nielsen, 2007; Wakeley, 2009), and both simulation and empirical studies indicate that the inferences obtained under models that ignore these processes can be biased in important ways (Eckert and Carstens, 2008; Leache et al., 2014; Burbrink and Guiher, 2015).

The second continuous-time approximation to the Wright-Fisher model is based on modeling , the proportion of alleles in a population of size at time , as a diffusion process. This process is defined through the following stochastic differential equation (SDE)


which will be referred to as the Wright-Fisher (WF) diffusion. In (1), the initial condition is , is a standard Brownian motion process and the coefficient is assumed to satisfy sufficient smoothness conditions such that (1) admits a unique weak solution. Such conditions can be found in Section 5.3 of Karatzas and Shreve (1998). In applications, the drift function is selected to incorporate complex evolutionary processes such as selection and mutation. Because such processes are easily specified in this approach, the WF diffusion has been increasingly applied to settings in which the number of populations is small (typically on the order of one to three populations) and primary interest is on inferring the contribution of these evolutionary processes to the current genetic configuration of the population(s).

Practical use of (1) is complicated by the fact that the transition density for the WF diffusion does not have a closed form. Current approaches using this model are therefore based on approximations of this density. A typical solution is to approximate the transition density using a numerical scheme for the SDE (see Kloeden and Platen, 1992, for a classical overview). Other recent approaches to estimating the transition density of an SDE include methods based on importance sampling (see Pedersen, 1995; Elerian et al., 2001; Brandt and Santa-Clara, 2002; Durham and Gallant, 2002; Lin et al., 2010, and the references therein for a review), closed-form approximations based on a Hermite polynomials expansion (Aït-Sahalia, 2002, 2008), a strategy based on approximating the solution of the Kolmogorov forward equation (Lo, 1988), and exact sampling (Beskos et al., 2006b). Another approach, specific to approximating the transition density of the WF diffusion, involves truncation of a spectral representation of the target density (Song and Steinrücken, 2012; Steinrücken et al., 2013).

Given the wealth of genetic data currently available as a result of ever-improving sequencing technologies, an additional challenge for the diffusion-based inferential framework is the desire to use data from many loci simultaneously for the inference of evolutionary process parameters. This is typically done by making the assumption that the loci included in the study are unlinked, and thus independent, as would typically be the case for single nucleotide polymorphism (SNP) data collected across the genome. Under this assumption, the likelihood is computed by multiplying single-locus likelihoods computed under the model across loci. When this assumption is violated and loci are actually linked, the likelihood is viewed as a composite likelihood, and inference proceeds as usual.

This framework has been considered in a variety of contexts. Williamson et al. (2005) considered specification of the drift coefficient to include the effect of selection, and fit models that allowed for changes in population size throughout time. To carry out the computations required by the WF diffusion, the Crank-Nicolson approximation was used to solve the relevant SDE in order to compute the likelihood. Estimation of the selection coefficient and population size parameters was then carried out in a likelihood framework. These authors later extended the method to include multiple populations with the possibility of gene flow between populations, again using a finite differencing method to provide a numerical solution to the SDE (Gutenkunst et al., 2009)

. As before, inference was carried out in a likelihood framework, with variance estimates and significance values for hypothesis testing obtained using bootstrapping. These methods were implemented in the program

for up to three populations (Gutenkunst et al., 2010), and have been applied to a variety of organisms, including humans, cattle, rice, and bees (Robinson et al., 2014). In addition, recent attention has been given to methodology for efficient computation of estimates of uncertainty in parameter estimates, with applications to model selection (Coffman et al., 2015; Robinson et al., 2014).

Despite a great deal of previous work in this area, inference in a Bayesian framework has been limited. Schraiber et al. (2016) and Ferrer-Admetlla et al. (2016) applied Bayesian techniques to infer selection in allele frequency time series data. Their methods require approximation of the transition density of the WF diffusion. Jenkins and Spanó (2015) recently developed an algorithm for obtaining exact draws from the WF diffusion with general drift function, using draws from the neutral WF diffusion with mutation. The work of Jenkins and Spanó (2015) thus enables the development of Markov chain Monte Carlo (MCMC) approaches to inference of evolutionary process parameters from genome-wide SNP data based on the composite likelihood, but without the need for numerical approximation of the solution to the SDE.

Here, we implement one such method for a single-population model in which interest is focused on estimation of the selection parameter. Unlike Schraiber et al. (2016), Ferrer-Admetlla et al. (2016), Foll et al. (2015), and Malaspinas et al. (2012), who considered time series data, we assume that allele frequency data are only available at a single time point. We note that previous work in this single time point setting has focused primarily on likelihood- or summary statistic-based inference, rather than on Bayesian approaches. In particular, Nielsen et al. (2005) proposed a composite likelihood statistic for the detection of selective sweeps that allowed estimation of the magnitude of the selection coefficient as well as the location of the sweep. This method, called SweepFinder, was later made more computationally efficient in an implementation called SweeD (Pavlidis et al., 2013). Similarly, the linkage disequilibrium-based method of Kim and Nielsen (2004) has recently been implemented with improved computational efficiency in the software OmegaPlus (Alachiotis et al., 2012). Finally, the methodology proposed in Živković et al. (2015) could also be applied to data sampled at a single time point. All of these approaches, however, involve computation of statistics (e.g., composite likelihoods) at various positions in the genome, with no naturally associated measure of variance, aside from those obtained via bootstrap procedures, as in Živković et al. (2015). We thus propose a Bayesian approach to the problem of data sampled at a single time point.

In the next section, we describe the details of our proposed model, as well as the computational approach we take to implement inference in an MCMC-based Bayesian framework. We apply the method to simulated data to study its performance, which provides insight into the role of the prior distribution in inference under the model. We then apply the method to an empirical data set (Zhou et al., 2011) consisting of SNP data on flies that were subjected to either a hypoxic environment or a control (e.g., normal oxygen) environment for 200 generations, with the goal of estimating the extent of selection for the flies that were oxygen-deprived. We compare our results with those of Ronen et al. (2013), who analyzed the same data. Finally, we discuss potential extensions of our results to more general settings, as well as the advantages of inference in the Bayesian setting for these problems.

2 Methods

Suppose that a total of individuals are sampled, and that the type of allele ( or ) is recorded for each individual at each of loci. At locus , the number of individuals with allele is denoted . Let

be the complete vector of data across all

loci. For each locus, is assumed to be an observation from a Binomial distribution with trials and probability of success. Note that represents the proportion of alleles in the population at locus at time – the time at which the sample is collected – however, we will suppress the dependence on throughout the manuscript, in order to simplify notation. Conditional on the vector , the likelihood is given by


Note that this assumes conditional independence of the data across loci; if the loci are not independent given the allele frequencies, then the likelihood in (2) can be viewed as a pseudo-likelihood. This is the same likelihood function used in other applications (e.g., Williamson et al., 2005).

The allele frequencies , are assumed to be realizations (at time ) from WF diffusion processes with starting values . Often, is a parameter of interest in the model, however, in this case we will assume that is fixed and known. The drift coefficient is specified to incorporate the evolutionary processes of interest and associated model parameters. The goal of the inference procedure is to estimate the evolutionary model parameters given the data and carry out standard inference (e.g., hypothesis testing).

Here we consider a model in which selection and mutation operate within a single population. Selection refers to the evolutionary mechanism in which one allele is preferred over other potential alleles to survive to the next generation, while mutation refers to the origination of new types of alleles within the population. Over time, both selection and mutation can act to alter the frequency of allele at a particular locus in the population beyond what is expected under the basic Wright-Fisher model that does not include these processes. The WF diffusion corresponding to the model with mutation and selection is given by


where is the population-scaled selection parameter and is the population-scaled mutation parameter (see Etheridge, 2011, Chapter 5). For each locus , define the density of the allele frequency at time under this WF diffusion by , where is the initial allele frequency at locus .

Assuming that , , and are known, we focus on estimating the extent of selection via inference of the parameter . We take a Bayesian viewpoint and base our statistical inference on the posterior density function


where we use the notation

and is a specified prior distribution on the selection parameter . Subsequently, we use an MCMC algorithm to explore the distribution given by (4).

To that end, note that a typical Metropolis within Gibbs (MwG) sampler would alternate between updating and , conditional on everything else, using some proposal-acceptance mechanism. However, a direct implementation of this procedure suffers from several issues. For example, during the “update given ” step, under the assumption that the proposal distribution for is symmetric, the MwG accept-reject mechanism requires evaluating the acceptance ratio

where and are the current and proposed states respectively. However, since the transition density for the WF diffusion does not have a closed form, the ratio above cannot be evaluated exactly in practice.

Most of the methods described in Section 1 for estimating come with a significant computational cost. Such methods are typically designed to work in a maximum likelihood estimation setting, where one would have a relatively small number of transition density evaluations. In a Bayesian setting however, simulating a Markov chain would require a very large number of transition density evaluations and thus many of these methods become impracticable.

In this manuscript we suggest approaches based on exact sampling. Our goal is to use the transition density in an MCMC setting and thus we require computational efficiency. We propose two sampling strategies: (i) an exact sampler which avoids evaluating the transition density entirely but suffers from some computational overhead and (ii) an efficient approximate sampler based on exact simulations of the diffusion process (3). In both cases, our method relies on the ability to obtain exact draws from the density . We explain how this is done in detail in the Appendix. Henceforth, we assume that we can simulate variates

Although such simulation is possible, it can be extremely time-consuming.

Exact Sampler. Our exact MCMC sampler uses a joint proposal mechanism as follows. Let be the current state of the Markov chain exploring the distribution (4) and consider simulating a new proposed state from the (joint) prior distribution as follows:

1:Draw ;
2:Given , draw for ;
3:Set ;
4:Return .
Algorithm 1 Joint simulation of from the prior model

In this case one can derive the acceptance ratio for to be

which does not depend on and can therefore be easily evaluated. Exact draws from are used to generate proposals, but it is not necessary to evaluate to compute the acceptance ratio. Thus, unlike existing methods, this strategy does not require one to approximate the transition density.

There are, however, a number of drawbacks that one should consider. Given that this algorithm is essentially an independent Metropolis-Hastings sampler, it is well known that it may mix very poorly. Further, the second sampling step requires one to simulate the SDE (3). As mentioned above and discussed in greater detail in Section 3, such simulation can be extremely time-consuming. Since each update in the MCMC procedure requires this sampling step, simulation of the SDE becomes a bottleneck in the algorithm. Due to the combination of poor mixing and a slow sampling step, we found this algorithm to be impracticable computationally, especially when studying a large number of loci simultaneously. Nonetheless, given sufficient computing resources, we believe this procedure would perform adequately.

Approximate MwG. A natural solution to the issues described above is to find an analytical expression for the transition density , which would allow one to implement a standard MwG algorithm. As discussed earlier, this is a difficult task. To take advantage of our ability to sample exactly from

, we suggest approximating the transition density using a kernel density estimate



is the density of the standard Gaussian distribution,

is a bandwidth parameter, and are independent and identically distributed draws from the density . We choose according to Scott’s Rule (Scott, 1992) and in our simulations we use , see Section 3. It is known that the mean integrated squared error in this case is

where the expected value is taken with respect to the Monte Carlo draws. Once the transition density has been approximated, the posterior distribution (4) can be replaced by




The distribution (6) can then be explored directly using a standard MwG algorithm.

This approximate MwG algorithm requires a kernel density estimate of for each potential set of proposed values for , , , and . For this to be feasible we must restrict ourselves to a finite set of plausible values for these parameters. We choose to treat both and as fixed and known, and to assume that and are contained in the finite sets and , respectively. The initial allele frequencies are generally not known in practice, so it is not reasonable to treat as fixed and known. Instead, we view each as an independent draw from a discrete distribution that puts positive mass on . This allows us to account for uncertainty in , but presents a challenge in that we cannot estimate and simultaneously. We therefore “integrate out” to leave as the only parameter to be estimated.

Ultimately, for each and , we use (5) to produce a kernel density estimate of the transition density and then sum over to compute an alternative estimate that does not depend on . Specifically, we compute


The density estimate replaces in (7) and the algorithm proceeds as if were known.

The upfront computational cost of this approach is substantial as one must simulate many draws from many transition densities in order to generate a sufficient number of accurate kernel density estimates with which to effectively explore the parameter space. However, once these density estimates have been produced, the MCMC algorithm is able to run efficiently. This is in contrast to the exact sampler, which does not have the same upfront costs but performs each update inefficiently. Further, by updating the parameters individually, conditional on all other parameters in the model, the MwG sampler mixes better than the exact sampler and therefore requires far fewer updates to obtain a sample of approximately independent draws from the target posterior.

3 Results

Simulated Data. In order to assess the viability of our approach we first applied it to simulated data with known selection parameter. Specifically, for each of three different selection parameters (, , and ) we simulated 500 data sets with individuals and loci. This was accomplished by using the algorithm described in the Appendix to simulate each from the SDE given in (3) and drawing each corresponding from the appropriate Binomial distribution. Model parameters were given values matching those used with the empirical fly data in the example described below. Specifically, for each data set we set and , and let follow the ending distribution of allele frequencies for the neutral fly population in the example.

We conducted two sets of simulations. In both cases we defined the set of plausible values for as and assumed a distribution for the initial allele frequencies that put mass on the set . As such, we computed kernel density estimates with one corresponding to each combination of and . Each of these density estimates was based on 10,000 draws from the true transition density, produced via exact simulation of the SDE in (3). We then summed over as in (8) to obtain an estimate for each of the values of in .

Generating the samples required to compute these density estimates took nearly 1,600 hours on a NVIDIA Tesla C2050 GPU computing card. However, a disproportionate amount of this time was required to generate samples for the largest values of . It took about 607 hours for alone, about 384 hours for , and about 215 hours for . In contrast, generating the required samples for all values of in the set took a total of just 32 minutes. Large values of require more time because the proposals in the acceptance-rejection algorithm are generated from a process with . Thus, as increases in magnitude there is a bigger discrepancy between the process used for proposals and the target process, so far fewer proposals are accepted.

For the first set of simulations we defined as a discretized version of the distribution of that was used to generate the data whereas in the second set of simulations we defined

as a discrete Uniform distribution over the set

. Altering the assumed distribution for allows us to assess the importance of this choice and its impact on inference for .

For all simulations we placed a discrete Uniform prior on so that for each . Our proposal density for was a discrete Uniform centered on the current value of that put positive mass on the five grid values on either side of the current value ( for and otherwise). Proposals for each came from a distribution where was the current value for and was the proposal variance. Since these proposal densities are symmetric they do not play a role in the acceptance ratios of the MwG algorithm. For each simulated data set, a Markov chain was run for 10,000 steps. The first 1,000 steps were discarded as burn-in and every tenth step thereafter was retained to form a sample of size 900 from the posterior distribution. Trace plots and autocorrelation plots for a random subsample of the replicates (not shown) suggest that the Markov chains converge to a stationary distribution within 1,000 steps and that retaining only every tenth step eliminates virtually all of the autocorrelation among values in the posterior sample.

For each of the three simulated values of (0.0, 5.5, and 11.0) and each choice of we generated a posterior sample for each of the 500 simulated data sets. The posterior mean of was then computed for each of these samples, yielding 500 estimates of for each simulated value of and choice of . Histograms of these posterior estimates when was chosen to approximate the initial allele frequency used to generate the data can be seen in Figure 1. For all three values of the empirical distribution of the estimator covers the true value of . However, the empirical distribution of the estimator is not always centered at the true value of

. Thus, this approach appears to be able to successfully distinguish among weak, strong, and moderate selection, but it may not yield an unbiased estimator of


Additionally, a 99% credible interval for

was computed from the posterior sample of each replicate. When , 97.2% of these intervals included the true value of , whereas for and the coverage was 99.2% and 84.4%, respectively. In all three cases the average length of the credible intervals was about three units. Thus, the lower coverage for is presumably due not to narrower intervals, but rather to the fact that those intervals tend to lie above . Ultimately, our approach may lead to slightly biased estimates, but we feel that we effectively quantify the uncertainty of the estimates.

Figure 1: Histograms showing estimates of the selection parameter from 500 simulations of data with true values of (white), (light gray), and (dark gray) when is chosen to reflect the distribution of initial allele frequencies used to generate the data. The vertical lines indicate the true values of .

Histograms of the posterior estimates for the cases of and when was chosen to be a discrete Uniform can be seen in Figure 2. The plot for is excluded because all of the Markov chains became trapped in a state with , which is the lower boundary of our set of plausible values for . This suggests that our posterior estimate for would have been less than if the algorithm had been able to search for lesser values of . Further, none of the 500 replicates for any of the simulated values of yields a single 99% credible interval that captures the true value of . It is evident from these results that our approach relies on an accurate assumption for the distribution of the initial allele frequencies. Assuming a discrete Uniform distribution for puts more mass on higher values of than the distribution used to generate the data did. Consequently, we estimate values for that are too small when we assume such a distribution for . We address the dependence of our method on an accurate choice for in Section 4.

Figure 2: Histograms showing estimates of the selection parameter from 500 simulations of data with true values of (light gray) and (dark gray) when is chosen to be Uniform over . The vertical lines indicate the true values of .

Example Data: Hypoxia in Flies. An example to which this methodology can be applied is a study conducted by Zhou et al. (2011) in which flies (Drosophila melanogaster) were subjected to a hypoxic environment or a normal oxygen environment for generations. Flies in the hypoxic environment were under selective pressure to adapt to the lack of oxygen in their environment. The data consist of allele frequencies at the end of the study among approximately 200 surviving flies in each population. These frequencies were measured at 75,999 loci along the 2L chromosome for the population under selective pressure and 71,335 loci along the 2L chromosome for the neutral population. A subsequent analysis of these data conducted by Ronen et al. (2013) identified a region of the 2L chromosome showing strong evidence of selection. Within this region we have measurements at 459 loci for the population under selective pressure and measurements at 324 loci for the neutral population.

We applied our proposed methodology to the region of the 2L chromosome of the fly data of Zhou et al. (2011) that Ronen et al. (2013) identified as being subject to selection. We used the same settings as in the first set of simulations including choosing to be a discretized version of the distribution of allele frequencies at the end of the study for the neutral population. Thus, we assumed that both the neutral population and the population under selective pressure began the study with a distribution of allele frequencies similar to that of the neutral population at the end of the study. As such, we should estimate a selection parameter near zero for the neutral population. A Markov chain was run to obtain a posterior sample of size 900 for each of the two populations of flies.

Figure 3 shows a histogram of the marginal posterior sample of for each population of flies. These posterior samples correspond to the region of the 2L chromosome identified by Ronen et al. (2013) as being subject to selection. The white bars show the posterior sample for the neutral population whereas the gray bars show the posterior sample for the population subjected to the hypoxic environment. For the neutral population the posterior mean for is and a corresponding 99% credible interval contains zero. Thus, as expected, there is no evidence of selection. For the population subjected to the hypoxic environment the posterior mean for is and the entire sample from the posterior density of lies above zero. This indicates strong positive selection for this population.

Figure 3: Histograms of the posterior sample of from fitting our model to the neutral fly population (white) and the fly population subjected to a hypoxic environment (gray).

Finally, to determine if any other regions of the 2L chromosome display evidence of selection we applied our proposed methodolgy to sliding windows of loci with each window shifted 250 loci from the preceding window. This was done for the population subjected to the hypoxic environment. Note that this procedure assumes that all loci within a sliding window have identical values for the selection coefficient, a condition that is unlikely to be strictly true. However, we view this as a reasonable approximation given that the relatively small window size for this fairly dense SNP sample means that SNPs within a single window are likely to be linked. In this situation, the selection coefficient estimated in our model can be viewed as an “average” value for the strength of selection across the loci included in the window. Thus, our method can be used to identify regions of the genome that are likely to be under selection, but does not necessarily pinpoint individual SNPs. For each window, a Markov chain was run for 10,000 steps and every tenth step after the first 1,000 was retained to form a posterior sample of size 900. These posterior samples were used to produce 99% credible intervals for .

Figure 4 displays 99% credible intervals for for the sliding windows along the 2L chromosome for the population subjected to the hypoxic environment. Intervals displayed in gray contain zero whereas those shaded black do not. We find evidence of selection in more regions than Ronen et al. (2013) did, including several regions under negative selection. However, we detect the strongest signal in the same region that Ronen et al. (2013) identified. In fact, the credible interval that lies furthest above zero falls directly in the relatively narrow region (highlighted by a gray vertical bar in Figure 4) that Ronen et al. (2013) identified using the XP-SFselect method. To further compare our method to existing approaches, we used the SweeD software (Pavlidis et al., 2013) with grid size 30,000 to compute the composite likelihood ratio (CLR) at various points along the chromosomal region (shown in medium-gray at the bottom of Figure 4 with scaling given by the right-handed y-axis). In general, the magnitude of the CLR statistic corresponds to the signal indicated by our credible intervals.

Figure 4: Results of our analysis of selection for the 2L chromosome of the fly population subjected to the hypoxic environment. Location (x-axis) is displayed in units of . The left-handed y-axis gives units for 99% credible intervals for for sliding windows along the chromosome. Intervals containing zero are gray whereas those not containing zero are black. The right-handed y-axis gives units for the composite likelihood ratio statistic computed by the SweeD software (Pavlidis et al., 2013), shown in medium-gray at the bottom of the plot. The vertical gray bar indicates a region previously identified by Ronen et al. (2013) as being under positive selection, which corresponds to the strongest signal of selection detected by our analysis.

4 Discussion

Here we propose a Bayesian method for inference of population genetics parameters under the Wright-Fisher diffusion model. In particular, we consider estimation of the selection coefficient from allele frequency data sampled at a single time point by constructing an MCMC algorithm that allows one to draw from the joint posterior distribution of the selection coefficient and the allele frequencies. We show that when assumptions about the initial allele frequencies are accurate our method performs well in simulations as well as in the analysis of an empirical data set on hypoxia in flies.

Existing methods for inference in similar settings have used a maximum likelihood framework and have relied on numerical approximation of the solution of the SDE corresponding to the WF diffusion. In addition to requiring significant computation, a drawback of such methods is that they do not provide a natural way to assess variability in the estimates. Recent work in this area has involved the development of bootstrap procedures for estimation of error variances (Coffman et al., 2015; Robinson et al., 2014), requiring significant additional computation time. An advantage of a Bayesian approach is that MCMC can be used to approximate the entire posterior distribution, providing an automatic assessment of uncertainty.

Bayesian approaches to this problem require an efficient method for sampling from the WF diffusion or for accurately approximating its transition density. The important work of Jenkins and Spanó (2015) provides a method to draw exact samples from WF diffusions that incorporate features of the evolutionary process such as mutation and selection. Their work will enable much more general developments in the application of Bayesian methods to population-level diffusion models. Our example on hypoxia in flies demonstrates the potential utility of these approaches; for example, we compute 99% credible intervals along a chromosomal region using a sliding window approach, and recover evidence for strong positive selection in a region identified in previous studies.

Although the recent work of Jenkins and Spanó (2015) enables MCMC-based Bayesian inference in this setting, the method we have developed based on this work is also fairly computationally intensive. This is because of the difficulty of obtaining exact draws from the general WF diffusion model. To deal with this, we chose to use kernel density estimates of the transition density, and most of our computational effort is devoted to computing these kernel density estimates. The computational effort needed for this step will vary depending on the true underlying parameters that generated the data. For our empirical data set for which , the method required a great deal of effort because is large compared to the neutral Wright-Fisher model, and sampling from the WF diffusion is inefficient. When the true value of is closer to zero, transition density estimates for larger values of will not be needed, and the overall computation time will be reduced. Also, although we have exploited GPU computing in order to make this approach computationally feasible, it is possible that more efficient algorithms could be developed to speed computation further. It may also be possible to use alternative approaches to the kernel density estimates we have used that might provide similar accuracy with reduced computational expense.

Our simulation studies clearly demonstrate that the choice of prior distribution for the initial allele frequencies across loci will have a strong effect on the posterior distribution. This is intuitive in the case of selection: a large initial allele frequency and strong negative selection can yield the same ending allele frequency distribution as a small initial allele frequency with strong positive selection. Thus, two markedly different processes can result in very similar data, and it is not possible to disentangle the effects of these two parameters. This means that a requirement for effective use of this method is either prior knowledge of the approximate allele frequency distribution at neutrality, or a thorough exploration of the sensitivity of the results to the choice of prior distribution. Fortunately, many studies in this area do devote effort to the assessment of neutral processes in the populations under consideration and thus there will often be prior information about the distribution of initial allele frequencies available to guide choice of the prior distribution (see, e.g., Li and Stephan, 2005; Williamson et al., 2005).

The method presented here addresses perhaps the simplest scenario one might consider in this setting: a single population and a single parameter of interest (the selection coefficient). In practice, however, data may be collected for multiple populations or species at multiple time points, and inference for all relevant parameters (e.g., divergence times, migration rates, and selection coefficients) may be desired. While diffusion models such as this are now commonly used to model up to three populations and to infer multiple parameters simultaneously, they are generally computationally demanding and, as mentioned above, the current likelihood-based approaches require significant additional computation in order to carry out complete statistical inference that includes an assessment of uncertainty. As genome-scale data become more widely available and sampling of individuals within populations allows more precise estimates of population allele frequencies throughout the genome, methods based on allele frequencies provide a unique opportunity for efficient population-level evolutionary inference. Methods that utilize a Bayesian framework for these problems have the potential to contribute in important ways to these developments, and our Bayesian framework provides a useful starting point for these future developments.

5 Acknowledgement

This work was supported by National Science Foundation grants DMS-1209142 and DMS-1407604 (R.H.) and DMS-1106706 (L.S.K.).


Exact simulation of the neutral Wright-Fisher diffusion

Consider a neutral WF diffusion process


where the drift function is , for and (note the absence of the selection parameter). The necessary details for simulating a variate from the model (9) for a fixed , given that are given in Jenkins and Spanó (2015) and we briefly review the main ideas. It is well-known (Griffiths, 1979; Tavaré, 1984; Ethier and Griffiths, 1993; Griffiths and Spanò, 2013) that, conditional on , the variate

has a probability density function given by


where , and and are the Binomial and Beta density functions, given by

The mixture weights in (10) have an alternating series expansion



with the convention that for and . Sampling from the density (10) is achieved via Agorithm 2 described below (see Jenkins and Spanó, 2015).

1:Draw from the discrete distribution ;
2:Draw ;
3:Draw ;
4:Return .
Algorithm 2 Exact Sampling for the neutral WF diffusion:

The computational burden in Algorithm 2 is in Step 1, since the weights do not have a closed form expression and (11) suggests an infinite amount of computation. However, given the alternating series representation in (11), one can use ideas described in Chapter 5 of Devroye (1986) to devise an exact procedure for simulating the variate , which only requires finite computing time. Details of such a procedure are provided in Section 3.2 of Jenkins and Spanó (2015).

Exact simulation of general WF-diffusion models

As we mention above, our approach relies on being able to obtain exact draws from a general Wright-Fisher diffusion (1). We briefly review the general setup for sampling an SDE. Assume that is a stochastic process described via


where the functions and are assumed to be smooth enough such that (12) has a unique weak solution. Such conditions can be found in Karatzas and Shreve (1998), for example. Our goal is to simulate exact draws from the distribution of , for some fixed , conditional on . In general, the transition density for the process is not available in closed form, not even in an infinite series representation as is the case with the neutral Wright-Fisher diffusion. Recently, in a series of papers (Beskos and Roberts, 2005; Beskos et al., 2006a, 2008), the authors have developed a rejection sampling approach, which yields an exact (distribution-wise) skeleton of the full path for a very large class of diffusions. Briefly, let and let be a typical element of . Let be the probability measure induced by onto and be another probability measure on which is user-selected in an appropriate way. Under regularity conditions (see Beskos and Roberts, 2005) one can use the Girsanov formula to derive and expression for the Radon-Nikodym derivative

where, in principle, one can arrange that . Thus, a rejection sampling strategy will be appropriate:

1:Sample ;
2:Accept as a draw from with probability .
Algorithm 3 Exact sampling for general diffusions

The proposal distribution has to be selected in a way that Step 1 of Algorithm 3 can be done efficiently, i.e. a biased Brownian Bridge measure. We refer the reader to Beskos and Roberts (2005) for a detailed description of all of the conditions and the general setup.

General Wright-Fisher diffusions. If the diffusion coefficient in (12) takes the form

then the SDE (12) is a general WF diffusion. As Jenkins and Spanó (2015) suggest, in this case, it is efficient to select the proposal measure to be the law induced by the neutral WF diffusion (9).


  • Aït-Sahalia (2002) Aït-Sahalia, Y. (2002). Maximum likelihood estimation of discretely sampled diffusions: A closed-form approximation approach. Econometrica, 70:223–262.
  • Aït-Sahalia (2008) Aït-Sahalia, Y. (2008). Closed-form likelihood expansions for multivariate diffusions. The Annals of Statistics, 36:906–937.
  • Alachiotis et al. (2012) Alachiotis, N., Stamatakis, A., and Pavlidis, P. (2012). OmegaPlus: a scalable tool for rapid detection of selective sweeps in whole-genome datasets. Bioinformatics, 28:2274–2275.
  • Beskos et al. (2006a) Beskos, A., Papaspiliopoulos, O., and Roberts, G. O. (2006a). Retrospective exact simulation of diffusion sample paths with applications. Bernoulli, 12:1077–1098.
  • Beskos et al. (2008) Beskos, A., Papaspiliopoulos, O., and Roberts, G. O. (2008). A factorisation of diffusion measure and finite sample path constructions. Methodology and Computing in Applied Probability, 10:85–104.
  • Beskos et al. (2006b) Beskos, A., Papaspiliopoulos, O., Roberts, G. O., and Fearnhead, P. (2006b). Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68:333–382.
  • Beskos and Roberts (2005) Beskos, A. and Roberts, G. O. (2005). Exact simulation of diffusions. The Annals of Applied Probability, 15:2422–2444.
  • Brandt and Santa-Clara (2002) Brandt, M. W. and Santa-Clara, P. (2002). Simulated likelihood estimation of diffusions with an application to exchange rate dynamics in incomplete markets. Journal of Financial Economics, 63:161–210.
  • Burbrink and Guiher (2015) Burbrink, F. T. and Guiher, T. J. (2015). Considering gene flow when using coalescent methods to delimit lineages of North American pitvipers of the genus Agkistrodon. Zoological Journal of the Linnean Society, 173:505–526.
  • Coffman et al. (2015) Coffman, A. C., Hsieh, P. H., Gravel, S., and Gutenkunst, R. N. (2015). Computationally efficient composite likelihood statistics for demographic inference. Molecular Biology and Evolution, 33:591–593.
  • Devroye (1986) Devroye, L. (1986). Non-Uniform Random Variate Generation. Springer-Verlag.
  • Durham and Gallant (2002) Durham, G. B. and Gallant, A. R. (2002). Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes. Journal of Business & Economic Statistics, 20:297–316.
  • Eckert and Carstens (2008) Eckert, A. J. and Carstens, B. C. (2008). Does gene flow destroy phylogenetic signal? The performance of three methods for estimating species phylogenies in the presence of gene flow. Molecular Phylogenetics and Evolution, 49:832–842.
  • Edwards (2009) Edwards, S. V. (2009). Is a new and general theory of molecular systematics emerging? Evolution, 63:1–19.
  • Elerian et al. (2001) Elerian, O., Chib, S., and Shephard, N. (2001). Likelihood inference for discretely observed nonlinear diffusions. Econometrica, 69:959–993.
  • Etheridge (2011) Etheridge, A. (2011). Some Mathematical Models from Population Genetics: École D’Été de Probabilités de Saint-Flour XXXIX-2009, volume 39. Springer Science & Business Media.
  • Ethier and Griffiths (1993) Ethier, S. N. and Griffiths, R. C. (1993). The transition function of a Fleming-Viot process. Annals of Probability, 21:1571–1590.
  • Ferrer-Admetlla et al. (2016) Ferrer-Admetlla, A., Leuenberger, C., Jensen, J. D., and Wegmann, D. (2016).

    An approximate Markov model for the Wright-Fisher diffusion and its application to time series data.

    Genetics, 203:831–846.
  • Fisher (1930) Fisher, R. A. (1930). The Genetical Theory of Natural Selection. Clarendon Press.
  • Foll et al. (2015) Foll, M., Shim, H., and Jensen, J. D. (2015). WFABC: a Wright-Fisher ABC-based approach for inferring effective population sizes and selection coefficients from time-sampled data. Molecular Ecology Resources, 15:87–98.
  • Griffiths (1979) Griffiths, R. C. (1979). A transition density expansion for a multi-allele diffusion model. Advances in Applied Probability, 11:310–325.
  • Griffiths and Spanò (2013) Griffiths, R. C. and Spanò, D. (2013). Orthogonal polynomial kernels and canonical correlations for Dirichlet measures. Bernoulli, 19:548–598.
  • Gutenkunst et al. (2010) Gutenkunst, R., Hernandez, R., Williamson, S., and Bustamante, C. (2010). Diffusion approximations for demographic inference: DaDi. Nature Precedings ¡¿.
  • Gutenkunst et al. (2009) Gutenkunst, R. N., Hernandez, R. D., Williamson, S. H., and Bustamante, C. D. (2009). Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genetics, 5:e10000695.
  • Hey and Nielsen (2007) Hey, J. and Nielsen, R. (2007). Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proceedings of the National Academy of Sciences USA, 104:2785–2790.
  • Jenkins and Spanó (2015) Jenkins, P. A. and Spanó, D. (2015). Exact simulation of the Wright-Fisher diffusion. arXiv:1506.06998v1.
  • Karatzas and Shreve (1998) Karatzas, I. and Shreve, S. (1998). Brownian motion and stochastic calculus, volume 113. Springer Science & Business Media.
  • Kim and Nielsen (2004) Kim, Y. and Nielsen, R. (2004). Linkage disequilibrium as a signature of selective sweeps. Genetics, 167:1513?1524.
  • Kingman (1982a) Kingman, J. F. C. (1982a). The coalescent. Stochastic Processes and their Applications, 13:235–248.
  • Kingman (1982b) Kingman, J. F. C. (1982b). Exchangeability and the evolution of large populations. In Koch, G. and Spizzichino, F., editors, Exchangeability in Probability and Statistics, pages 97–112. North-Holland, Amsterdam.
  • Kingman (1982c) Kingman, J. F. C. (1982c). On the genealogy of large populations. Journal of Applied Probability, 19A:27–43.
  • Kloeden and Platen (1992) Kloeden, P. E. and Platen, E. (1992). Numerical Solution of Stochastic Differential Equations (Stochastic Modelling and Applied Probability). Springer-Verlag.
  • Knowles and Kubatko (2010) Knowles, L. L. and Kubatko, L. S. (2010). Estimating Species Trees: Practical and Theoretical Aspects. Wiley-Blackwell.
  • Leache et al. (2014) Leache, A. D., Harris, R. B., Rannala, B., and Yang, Z. (2014). The influence of gene flow on species tree estimation: A simulation study. Systematic Biology, 63:17–30.
  • Li and Stephan (2005) Li, H. and Stephan, W. (2005). Maximum-likelihood methods for detecting recent positive selection and localizing the selected site in the genome. Genetics, 171:377–384.
  • Lin et al. (2010) Lin, M., Chen, R., and Mykland, P. (2010). On generating Monte Carlo samples of continuous diffusion bridges. Journal of the American Statistical Association, 105:820–838.
  • Lo (1988) Lo, A. W. (1988). Maximum likelihood estimation of generalized Itô processes with discretely sampled data. Econometric Theory, 4:231–247.
  • Malaspinas et al. (2012) Malaspinas, A.-S., Malaspinas, O., Evans, S. N., and Slatkin, M. (2012). Estimating allele age and selection coefficient from time-serial data. Genetics, 192:599–607.
  • Nielsen and Wakeley (2001) Nielsen, R. and Wakeley, J. (2001). Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics, 158:885–896.
  • Nielsen et al. (2005) Nielsen, R., Williamson, S., Kim, Y., Hubisz, M. J., Clark, A. G., and Bustamante, C. (2005). Genomic scans for selective sweeps using snp data. Genome Research, 15:1566?1575.
  • Pavlidis et al. (2013) Pavlidis, P., Živković, D., Stamatakis, A., and Alachiotis, N. (2013). SweeD: Likelihood-based detection of selective sweeps in thousands of genomes. Molecular Biology and Evolution, 30:2224–2234.
  • Pedersen (1995) Pedersen, A. R. (1995). A new approach to maximum likelihood estimation for stochastic differential equations based on discrete observations. Scandinavian Journal of Statistics, 22:55–71.
  • Robinson et al. (2014) Robinson, J. D., Coffman, A. J., Hickerson, M. J., and Gutenkunst, R. N. (2014). Sampling strategies for frequency spectrum-based population genomic inference. BMC Evolutionary Biology, 14:254.
  • Ronen et al. (2013) Ronen, R., Upda, N., Halperin, E., and Bafna, V. (2013). Learning natural selection from the site frequency spectrum. Genetics, 195:181–193.
  • Schraiber et al. (2016) Schraiber, J. G., Evans, S. N., and Slatkin, M. (2016). Bayesian inference of natural selection from allele frequency time series. Genetics, 203:493–511.
  • Scott (1992) Scott, D. W. (1992). Multivariate Density Estimation: Theory, Practice, and Visualization. John Wiley.
  • Song and Steinrücken (2012) Song, Y. S. and Steinrücken, M. (2012). A simple method for finding explicit analytic transition densities of diffusion processes with general diploid selection. Genetics, 190:1117–1129.
  • Steinrücken et al. (2013) Steinrücken, M., Wang, Y. X. R., and Song, Y. S. (2013). An explicit transition density expansion for a multi-allelic Wright-Fisher diffusion with general diploid selection. Theoretical Population Biology, 83:1–14.
  • Tavaré (1984) Tavaré, S. (1984). Line-of-descent and genealogical processes, and their applications in population genetics models. Theoretical Population Biology, 26:119–164.
  • Živković et al. (2015) Živković, D., Steinrücken, M., Song, Y. S., and Stephan, W. (2015). Transition densities and sample frequency spectra of diffusion processes with selection and variable population size. Genetics, 200:601–617.
  • Wakeley (2009) Wakeley, J. (2009). Coalescent Theory: An Introduction. Roberts and Company.
  • Wakeley and Hey (1997) Wakeley, J. and Hey, J. (1997). Estimating ancestral population parameters. Genetics, 145:847–855.
  • Williamson et al. (2005) Williamson, S. H., Hernandez, R., Fledel-Alon, A., Zhu, L., Nielsen, R., and Bustamante, C. D. (2005). Simultaneous inference of selection and population growth from patterns of variation in the human genome. Proceedings of the National Academy of Sciences, 102:7882–7887.
  • Wright (1931) Wright, S. (1931). Evolution in Mendelian populations. Genetics, 16:97–159.
  • Zhou et al. (2011) Zhou, D., Upda, N., Gersten, M., Visk, D. W., Bashir, A., and et al. (2011). Experimental selection of hypoxia-tolerant Drosophila melanogaster. Proceedings of the National Academy of Sciences, 108:2349–2354.