記住我
Allele frequencies vary across populations and loci, even in the presence of migration. While most differences may be due to genetic drift, divergent selection will further increase differentiation at some loci. Identifying those is key in studying local adaptation, but remains statistically challenging. A particularly elegant way to describe allele frequency differences among populations connected by migration is the F-model, which measures differences in allele frequencies by population specific FST coefficients. This model readily accounts for multiple evolutionary forces by partitioning FST coefficients into locus- and population-specific components reflecting selection and drift, respectively. Here we present an extension of this model to linked loci by means of a hidden Markov model (HMM), which characterizes the effect of selection on linked markers through correlations in the locus specific component along the genome. Using extensive simulations, we show that the statistical power of our method is up to twofold higher than that of previous implementations that assume sites to be independent. We finally evidence selection in the human genome by applying our method to data from the Human Genome Diversity Project (HGDP).
MIGRATION is a major evolutionary force homogenizing evolutionary trajectories of populations by promoting the exchange of genetic material. At some loci, however, the influx of new genetic material may be modulated by selection. In case of strong local adaptation, for instance, migrants may carry maladapted alleles that are selected against. Identifying loci that contribute to local adaptation is of major interests in evolutionary biology because these loci are thought to constitute the first step toward ecological speciation (e.g., Wu 2001; Feder et al. 2012) and allow us to understand the role of selection in shaping phenotypic differences between populations and species (e.g., Bonin et al. 2006; Fournier-Level et al. 2011).
A simple, yet flexible and useful, approach to identify loci contributing to local adaptation is to scan the genome using statistics that quantify divergence between populations. One frequently used statistic is FST, which measures population differentiation, and loci with greatly elevated FST have been reported for many population comparisons (e.g., Jones et al. 2012; Andrew and Rieseberg 2013; Stölting et al. 2013). While other statistics measuring absolute divergence (Cruickshank and Hahn 2014) or incongruence between a population tree and locus-specific genealogies (Durand et al. 2011; Peter 2016) may be more suited in some situations, genome scans suffer from two inherent limitations. First, multiple evolutionary scenarios may explain the deviations in those statistics, making interpretation difficult (e.g., Cruickshank and Hahn 2014; Eriksson and Manica 2012). Second, the definition of outliers is arbitrary, allowing for the detection of candidate loci only. Indeed, loci also vary in their divergence between populations that were never subjected to selection, but outlier approaches would still identify outliers.
Multiple methods have thus been developed that explicitly incorporate the stochastic effects of genetic drift. A first important step to improve the reliability of outlier scans was the proposal to compare observed values of such statistics against the distribution expected under a null model. Among the first, Beaumont and Nichols (1996) proposed to obtain the distribution of FST through simulations performed under an island model. While the idea to evidence selection by comparing FST to its expectations is far from new (e.g., Lewontin and Krakauer 1973), the difficulty of properly parameterizing the null model was quickly realized (e.g., Nei and Maruyama 1975). The success of the method by Beaumont and Nichols (1996) relies on tailoring the parameters of the underlying island model to match the observed heterozygosity at each locus—an approach that is also easily extended to structured island models (Excoffier et al. 2009).
A more formal approach is given by means of the F-model (Rannala and Hartigan 1996; Balding 2003; Falush et al. 2003; Gaggiotti and Foll 2010), under which allele frequencies are measured by locus and population specific 
coefficients that reflect the amount of drift that occurred in population j at locus l since its divergence from a common ancestral population. In the case of biallelic loci, the current frequencies 
are then given by a beta distribution (Beaumont and Balding 2004)
(1)where pl are the frequencies in the ancestral population and θlj is given by
It is straightforward to extend this model to account for different evolutionary forces that affect the degree of genetic differentiation. For instance, Beaumont and Balding (2004) proposed to partition the effects of genetic drift and selection into locus-specific and population-specific components αl and βj, as well as a locus-by-population specific error term γij:
(2)Loci with αl ≠ 0 are interpreted to be affected by either balancing (αl < 0) or divergent (αl > 0) selection, either because they are targets of selection or through hitch-hiking (Beaumont and Balding 2004). Such loci may be identified by contrasting models with αl = 0 or αl ≠ 0 for each locus l, either through Bayesian variable selection (Riebler et al. 2008) or via reversible-jump MCMC, as is done in the popular software BayeScan (Foll and Gaggiotti 2008).
A common problem of this, and many other, genome-scan methods is the assumption of independence among loci, which is easily violated when working with genomic data. By evaluating information from multiple linked loci jointly, however, the statistical power to detect outlier regions is likely increased considerably. Indeed, even a weak signal of divergence may become detectable if it is shared among multiple loci. Similarly, false positives may be avoided as their signal is unlikely to be shared with linked loci.
Unfortunately, fully accounting for linkage is often statistically challenging as well as computationally very costly. One solution is to split the problem by first inferring haplotypes for each sample, and then performing selection scans on the haplotype structure. The extended haplotype homozygosity (EHH) and its derived statistics (Sabeti et al. 2002; Voight et al. 2006; Sabeti et al. 2007; Tang et al. 2007), for instance, identify shared haplotypes of exceptional length. More recently, Fariello et al. (2013) introduced methods that identify haplotype clusters with particularly large frequency differences between populations and showed that using haplotypes rather than single markers increases power substantially.
An alternative solution is to model linkage through the autocorrelation of hierarchical parameters along the genome, which does not require knowledge of the underlying haplotype structure. Boitard et al. (2009) and Kern and Haussler (2010), for instance, proposed a genome-scan method in which each locus was classified as selected or neutral, and then used a Hidden Markov Model (HMM) to account for the fact that linked loci likely belonged to the same class, while ignoring autocorrelation in the genetic data itself.
Here, we build on this idea to develop a genome-scan method based on the F-model. While an HMM implementation of the F-model was previously proposed to deal with linked sites when inferring admixture proportions (Falush et al. 2003), we use it here to characterize autocorrelations in the strength of selection αl among linked markers. As we show using both simulations and an application to human data, aggregating information across loci results in an increase in power of up to twofold at the same false-discovery rate (FDR).
Materials and MethodsA model for genetic differentiation and observationsWe assume the classic F-model, in which J populations diverged from a common ancestral population. Since divergence, each population experienced genetic drift at a different rate. We quantify this drift of population 
at locus 
by θjl. We further assume each locus to be biallelic with ancestral frequencies pl, in which case the current frequencies 
are given by a beta distribution (Beaumont and Balding 2004), as shown in (2). We thus have
(3)where ql = 1 − pl, 
, 
and 
is the gamma function.
Let njl denote the allele counts in a sample of Njl haplotypes from population j at locus l, which is given by a binomial distribution
and hence
(4)Equations (3) and (4) combine to a beta-binomial distribution
(5)Model of selectionWe decompose θjl into a population-specific component βj shared by all loci, and a locus-specific component αl shared by all populations:
Here, the locus-specific component αl quantifies an excess or dearth of differentiation, which is attributed to the effect of either divergent or balancing selection, respectively (Beaumont and Balding 2004). Note that we adopt here the formulation of Foll and Gaggiotti (2008) and omit the error term γij of the original model of Beaumont and Balding (2004) shown in (2), as there is generally not enough information to estimate these parameters from the data (Beaumont and Balding 2004).
To account for autocorrelation among the locus-specific component, we propose to discretize αl = α (sl), where 
are the states of a ladder-type Markov model with m = 2smax + 1 states such that
(6)for some positive parameters αmax. The transition matrix of this Markov model shall be a finite-state birth-and-death process
(7)with elements 
denoting the probabilities to go from state i at locus l − 1 to state j at locus l given the strength of autocorrelation measured by the positive scaling parameter κ and the known distance dl between these loci, either in physical or in recombination space. Here, L is the m × m generating matrix
where the middle row at position smax + 1 reflects neutrality, and is given by the element
As exemplified in Figure 1, the two parameters μ and ν control the distribution of sites affected by selection (i.e., having αl ≠ 0) in the genome, with ν affecting the number of selected regions and μ their extent and selection strength, with higher values leading to more sites affected selection. It is important to note that we do not assume all sites with αl ≠ 0 to be targets of selection. Instead, many will be linked to a target of selection and experience αl ≠ 0 due to hitch-hiking.
Figure 1 (A) Expected proportion of neutral sites as a function of rates μ and ν. (B, C) Example paths of αl along 1000 loci simulated at a distance of dl = 100 with smax = 10 positive and negative states up to αmax = 3.0. Autocorrelation among loci was simulated with log(κ) = −3.0, v = 0.02, and μ = 0.91 (B, square) or μ = 0.74 (C, circle), respectively. The two cases correspond to an expected proportion of 20% and 10% of the genome under selection, as marked in (A).
The stationary distribution of this Markov chain is given by
with
Note that, as κ → ∞, our model approaches that of Foll and Gaggiotti (2008) implemented in BayeScan, but with discretized αl.
Hierarchical island models, first introduced by Slatkin and Voelm (1991), address the fact that divergence might vary among groups of populations. They were previously used to infer divergent selection, both using a simulation approach (Excoffier et al. 2009) as well as in the case of F-models (Foll et al. 2014). Here, we describe how our model is readily extended to additional hierarchies.
Consider G groups each subdivided into Jg populations with population-specific allele frequencies 
that derive from group-specific frequencies pgl as described above with group-specific parameters μg, vg, and κg. Analogously, we now assume group-specific frequencies to have diverged from a global ancestral frequency Pl according to locus-specific and group-specific parameters Θgl. Specifically,
such that
(8)where Ql = 1 − Pl and qgl = 1 − pgl. The parameter Θgl is given by
(9)As above, Bg quantifies group specific drift, 
are the states of a Markov model with m states and transition matrix 
with parameters μ and ν, a positive scaling parameter κ, and A (Sl) and Amax defined as in (6). Hence, we assume independent HMM models of the exact same structure at both levels of the hierarchy, as outlined in Figure 2. Additional levels could be added analogously.
Figure 2 A directed acyclic graph (DAG) of the proposed model with two hierarchical levels.
InferenceWe developed a Bayesian inference scheme for the parameters of the proposed model using a Markov chain Monte Carlo (MCMC) approach with Metropolis–Hastings updates, as detailed in the Supplementary Material. As priors, we used
Following Beaumont and Balding (2004), we used μb = −2 and 
throughout. We further set ap = bp = 1.
To identify candidate regions under selection, we used MCMC samples to determine the FDRs
for divergent and balancing selection, respectively, where 
and 
denote the full data.
We implemented the proposed Bayesian inference scheme in the easy-to-use C++ program Flink.
Given the heavy computational burden of the proposed model, we introduce several approximations. Most importantly, we group the distances dl into E + 1 ensembles such that el = log2dl, 
and use the same transition matrix Q(2e) for all loci in ensemble e. We then calculate Q(1) for the first ensemble using the computationally cheap yet accurate approximation
with 
where D = 2smax + 1 is the dimensionality of the transition matrix (Ferrer-Admetlla et al. 2016). The transition matrices of all other ensembles is obtained through the recursion 
(See Supplementary Information for other details regarding the implementation).
To quantify the benefits of accounting for autocorrelation in the locus specific components αl among linked loci, we first compared our method implemented in Flink against the method implemented in BayeScan (Foll and Gaggiotti 2008) on simulated data. All simulations were conducted under the model laid out above for a single group, using routines available in Flink and with parameter settings similar to those used in (Foll and Gaggiotti 2008). Specifically, we focused on a reference simulation in which we sampled N = 50 haplotypes from J = 10 populations with βj chosen such that 
in the neutral case (αl = 0). Following Foll and Gaggiotti (2008), we simulated all 
and about 20% of sites affected by selection (i.e., with αl ≠ 0) by setting μ = 0.91 and v = 0.02. We further set smax = 10 (resulting in m = 21 states) and 
and simulated 103 loci for each of 10 chromosomes, with a distance of 100 positions between adjacent sites and strength of autocorrelation log(κ) = −3. We then varied the number of populations J, the sample size N, 
or the strength of autocorrelation κ individually, while keeping all other parameters constant (Table 1). We further added a case without linkage (i.e., κ → ∞) by simulating each locus on its own chromosome.
To infer parameters with Flink, we set smax and αmax to the true values and ran the MCMC for 7·105 iterations, of which we discarded the first 20·105 as burn-in. During the chain, we recorded parameter values every 100 iterations as posterior samples. To infer parameters with Bayescan, we used version 2.1 and set the prior odds for the neutral model to 50, which we found to result in the same power as Flink in the reference simulation (see below) and in the absence of linkage (κ → ∞). We identified loci under selection at an FDR threshold of 5% for both methods.
Power of inference:We first evaluated the power of Flink in inferring the hierarchical parameters βj, ν, μ, and κ. As shown through the distributions of posterior means across all simulations, these estimates were very accurate and unbiased, regardless of the parameter values used in the simulations (Figure 3). This suggests that the power to identify selected loci is not limited by the number of loci we used to infer hierarchical parameters.
Figure 3 Boxplot of the parameters β1 (left), ν and μ (center), and log(κ) (right). The values are obtained from the mean of the posterior distributions obtained using Flink on the 10 simulations run for each of the set of parameters reported in Table 1. The red dotted lines show the true values of the respective parameters.
We next studied the impact of the sample size and the strength of population differentiation on the power (the true positive rate) to identify loci affected by selection (i.e., loci with αl ≠ 0). In line with findings reported by Foll and Gaggiotti (2008), power generally increased with 
, the number of sampled haplotypes, and the number of sampled populations (Figure 4, A–C). Larger sample sizes or stronger differentiation was particularly relevant for detecting loci under balancing selection, for which the power was generally lower and virtually zero at low differentiation 
or if only few populations were sampled (J = 2). Importantly, the FDR was below the chosen 5% threshold in 100% and 98.6% of all simulations conducted for loci identified as affected by divergent and balancing selection, respectively (Supplemental Material, Figure S1). The false positive rates (FPR) for these classes was < 0.1% in 98.6% and 97.1% of all simulations (see Figure 4 for neutral sites).
Figure 4 The true positive rate (power) in classifying loci as neutral (black) or under divergent (orange) or balancing selection (blue) as a function of the FST between populations (A), the number of haplotypes N (B), the number of populations J (C), and the strength of autocorrelation κ (D). Lines indicate the mean and range of true positive rates obtained with Flink (solid) and BayeScan (dashed) across 10 replicate simulations. Filled dots and the vertical gray line indicate the reference simulation shown in each plot.
Compared to BayeScan run on the same set of simulations, Flink had a higher power at the same FDR across all simulations, and often considerably so, unless if very many populations were sampled (Figure 4). If J = 10 populations were sampled, for instance, the power of Flink was about 0.2 higher for loci under divergent selection, and even up to 0.4 higher for those under balancing selection (Figure 4, A and B). Importantly, this increase in power described here is fully explained by Flink accounting for autocorrelation among the αl values as we chose the prior odds in BayeScan to result in the same power if the strength of autocorrelation vanishes (i.e., κ → ∞). Exploiting information from linked sites to identify divergent or balancing selection can thus strongly increase power, certainly if linkage extends to many loci. This is maybe best illustrated by the much higher power of Flink to identify loci under balancing selection at low differentiation (
, Figure 4A), in which case even many neutral loci are expected to show virtually no difference in allele frequencies and only an aggregation of loci with a subtle reduction in 
can be interpreted as a reliable signal for a target of selection in the region (Foll and Gaggiotti 2008).
Thanks to careful optimization, there is little to no overhead of our implementation compared to that of BayeScan. On the reference simulation of 104 loci from 10 populations, for instance, Flink took on average 130 min on a modern computer if calculations were spread over four CPU cores. On the same data, BayeScan took 361 min. However, we note that comparing the two implementations is difficult due to many settings that strongly impact run times, such as the number of iterations or the use of pilot runs in BayeScan. Without pilot runs, the run time of BayeScan reduced to 182 min on average for the default number of iterations (105 including burn-in). In the same time, Flink runs for close to 106 iterations, but also requires more to converge.
But since computation times scale linearly with the number of loci, they remain prohibitively slow for whole genome applications in a single run. However, computations are easily spread across many computers by analyzing the genome in independent chunks such as for each chromosome or chromosome arm independently. This is justified because (1) linkage does not persist across chromosome boundaries and is usually weak across the centromere, and (2) because our simulations indicate that 104 polymorphic loci were sufficient to estimate the hierarchical parameters accurately.
Effect of model misspecificationThe F-model makes the explicit assumption that the allele frequencies in a structured population can be characterized by a multinomial Dirichlet distribution. This distribution is appropriate for a wide range of demographic models, but not if some pairs of populations share a more recent ancestry than others (Beaumont and Balding 2004; Excoffier et al. 2009). Unsurprisingly, several previous studies found high FPRs when challenging BayeScan with models of isolation-by-distance (IBD), recent range expansions, recent admixture, or asymmetric divergence (e.g., Lotterhos and Whitlock 2014; Luu et al. 2017). These high FPRs are partially mitigated by choosing higher prior odds (e.g., 50 as used here, Lotterhos and Whitlock 2014) or when using the hierarchical version of BayeScan (Foll et al. 2014), particularly in case of asymmetric divergence. In the case of a recent range expansion or recent admixture, however, the F-model is unlikely to be appropriate and other methods have been shown t
留言 (0)