PaperlayerView exact revision on arXiv

Source-linked AI summary

WFABC: a Wright-Fisher ABC-based approach for inferring effective population sizes and selection coefficients from time-sampled data

Matthieu Foll, Hyunjin Shim, Jeffrey D. Jensen

arXiv:1409.8529v1q-bio.PE

TL;DR

As time-sampled sequencing becomes more available, methods for estimating population genetic parameters from these data remain limited. This paper implements an ABC-based software approach and reports that a nearly fully recessive lethal model best explains medionigra variation in Panaxia dominula.

  • Problem

    The paper notes that relatively few time-sampled datasets were available because genetic markers were limited.

    Source: Introduction · p. 3

  • Method

    The authors implement and expand an ABC approach in software, using the Jorde and Ryman Fs′ statistic as a summary statistic for Ne.

    Source: Introduction · p. 4, Ne-‐based ABC method · p. 5

  • Results

    For the P. dominula medionigra data, a nearly fully recessive lethal model is the best explanation for the observed variation.

    Source: Discussion · p. 21

  • Takeaways & Limitations

    The approach enables estimation of genome-wide average Ne from time-sampled data, including a single allele trajectory.

    Source: Discussion · p. 20

    Likelihood-method approximations are limited to particular effective-population-size and sampling settings, with a stated boundary involving 5000.

    Source: Discussion · pp. 19–20

Abstract

from arXiv · show

With novel developments in sequencing technologies, time-sampled data are becoming more available and accessible. Naturally, there have been efforts in parallel to infer population genetic parameters from these datasets. Here, we compare and analyze four recent approaches based on the Wright-Fisher model for inferring selection coefficients (s) given effective population size (Ne), with simulated temporal datasets. Furthermore, we demonstrate the advantage of a recently proposed ABC-based method that is able to correctly infer genome-wide average Ne from time-serial data, which is then set as a prior for inferring per-site selection coefficients accurately and precisely. We implement this ABC method in a new software and apply it to a classical time-serial dataset of the medionigra genotype in the moth Panaxia dominula. We show that a recessive lethal model is the best explanation for the observed variation in allele frequency by implementing an estimator of the dominance ratio (h).

Introduction

Time-sampled sequencing is becoming more common and may improve inference of population genetic parameters, but methods for these data—and especially incorporating selection—remain limited. The authors introduce ABC software for inferring per-site effective population sizes and selection coefficients, and compare it with existing likelihood-based methods.

Introduction

  • Time-sampled sequencing has become more common with technological advances, and allele trajectories may improve inference compared with single-time-point analyses.Source: Introduction · p. 3
  • Few methods estimate population genetic parameters from time-sampled data, and efforts to incorporate selection into those procedures have only recently begun.Source: Introduction · p. 3
  • The expanded availability of temporal data offers an opportunity to distinguish genetic drift from selection with greater precision and power.Source: Introduction · p. 3
  • The authors introduce software implementing and expanding an ABC approach to infer per-site Ne and s from time-sampled data.Source: Introduction · p. 4
  • They compare the ABC approach with existing likelihood-based methods to guide users toward suitable methods for different datasets.Source: Introduction · p. 4

Materials and Methods

This section describes WFABC’s two-step ABC approach for estimating shared effective population size and locus-specific selection, alongside likelihood-based alternatives and simulated comparisons that vary ascertainment and sampling conditions.

Ne-‐based ABC method

  • WFABC estimates shared Ne from genome-wide Fs’ statistics, then simulates Wright–Fisher trajectories using that Ne distribution and retains those closest to observed locus-specific summaries to infer selection.Fs’ values are averaged across sites and times, with a Bayesian bootstrap used to obtain a distribution for Ne.Source: Ne-‐based ABC method · pp. 4–5

Likelihood-‐based methods

  • The three likelihood-based HMM methods differ in scope: Bollback et al. (2008) estimates s and Ne, Malaspinas et al. (2012) also estimates allele age, and Mathieson and McVean (2013) estimates s assuming known Ne.Their likelihood calculations use numerical approximations, an EM algorithm, or a one-step diffusion process, respectively.Source: Likelihood-‐based methods · pp. 6–7

Simulated datasets for testing

  • Simulations can readily implement ascertainment schemes that are difficult for likelihood methods, including conditioning loci to remain polymorphic today and observing a minimum allele frequency across the trajectory.The Malaspinas et al. (2012) conditional implementation removes absorbing states 0 and 1.Source: Simulated datasets for testing · pp. 8–9
  • The likelihood methods use different ascertainment processes that produce data with differing informativeness, preventing direct performance comparison; they were therefore compared separately with WFABC.The first ascertainment case, observing any trajectory, is available in datasets but not to likelihood-based methods.Source: Simulated datasets for testing · p. 9
  • Wright–Fisher simulations used 1,000 replicates at Ne = 200, 1,000, or 5,000, with s sets (0, 0.005, 0.01, 0.015, 0.02) and (0.1, 0.2, 0.3, 0.4); WFABC retained the best 1% of 500,000 simulations.Using Euclidean distance between observed and simulated Fsd’ and Fsi’ statistics, the rejection ABC estimate of s was the posterior-distribution mean.Source: Simulated datasets for testing · p. 9
  • The simulations assessed WFABC performance across different sampling time points and sample sizes, and examined ascertainment using an unconditional case with 10% initial minor-allele frequency and a case where a new mutation occurs in the first generation.Source: Simulated datasets for testing · p. 10
  • For the realistic last-time-point polymorphism comparison, simulations adjusted generation counts to allow that condition, introduced a mutation in generation one, and drew 100 binomial samples at 12 time points.WFABC was compared with Malaspinas et al. (2012) and Bollback et al. (2008), using separate s search ranges for small and large selection values.Source: Simulated datasets for testing · p. 10
  • WFABC uses multiple sites to estimate Ne, whereas Bollback et al. (2008) and Malaspinas et al. (2012) estimate parameters from a single allele-frequency trajectory, which contains limited information about Ne.For comparison, Ne was fixed to its true value in all scenarios to evaluate the methods.Source: Simulated datasets for testing · p. 11
  • A final multilocus scenario inferred both s and Ne: with Ne fixed at 1,000, 10,000 trajectories were simulated, 500 under selection; WFABC used all trajectories to estimate Ne before s at each locus.WFABC used a uniform s prior of [−0.2, 0.6].Source: Simulated datasets for testing · p. 11

Results

The Results evaluate WFABC’s performance across selection, population-size, ploidy, and multilocus settings, alongside its application to moth data. They also examine how method efficiency varies with whether candidate mutations and population sizes are known.

Performance of the examined estimation procedures

  • With only 2 sampling time points, WFABC estimated s accurately for moderate s<0.2; estimates improved with more sampling time points and larger sample sizes.Source: Performance of the examined estimation… · p. 12
  • When datasets were conditioned on a minimum frequency of 5% at one time point, WFABC performed better than Mathieson and McVean (2013) for both small and large s.Mathieson and McVean (2013) consistently overestimated small s values, while their underestimation for large s compensated for that bias.Source: Performance of the examined estimation… · p. 12
  • For s values larger than 0.1, WFABC performed significantly better than the Malaspinas et al. (2012) approach in accuracy and precision, producing estimates with 10-fold less bias and 2-fold less error.Source: Performance of the examined estimation… · p. 13
  • Bollback et al. (2008) consistently estimated s=0 because the global likelihood maximum was near zero, so the method was excluded from further analyses and the data application.A narrow search range around the true s=0.4 could recover a local maximum, but the full likelihood surface placed the global maximum near zero.Source: Performance of the examined estimation… · pp. 13–14
  • For Ne in [200,1000], the two methods estimated the selection coefficient with high accuracy and precision; at Ne=5000 and s>0.01, 1000 replicates were not complete.Source: Performance of the examined estimation… · pp. 14–15
  • 0.049 vs. 0.10 was the selection-coefficient RMSE across 500 selected trajectories for WFABC and Malaspinas et al. (2012), respectively; WFABC excelled in large-s, large-Ne diploid, haploid, and multilocus cases.The summary identifies WFABC as superior when both s and Ne are large, for any haploid cases, and when multiple loci are available.Source: Performance of the examined estimation… · p. 15

Comparison of computational efficiency

Data Application

  • To estimate h, ABC added persistence summaries: tl counts generations below 5% without loss and th above 95% without fixation; moth data yielded tl=54 and th=0.The number of generations that an allele persists at low frequency is informative about dominance.Source: Data Application · p. 18
  • At 2Ne=1000, the joint-posterior mode was s=-1 and h=0.043, supporting lethal bimacula and deleterious medionigra phenotypes (wAA=0, wAa=0.96).The inferred medionigra heterozygote fitness was wAa=0.96.Source: Data Application · p. 18
  • The joint posterior indicates medionigra is either very strongly selected against and almost completely recessive (h≈0), or co-dominant with weaker selection (s≈-0.2).Density is greater in the recessive-lethal region, but the two-dimensional 90% highest posterior region includes the alternative co-dominant state.Source: Data Application · p. 18
  • The joint posterior mode is s=-1 and h=0, but its two-dimensional 90% highest posterior region includes s=0, consistent with a small enough population size allowing genetic drift alone to explain the observed pattern.Source: Data Application · p. 18

Discussion

In simulations, estimating genome-wide Ne before per-site selection coefficients improved selection-estimate precision, while WFABC and likelihood methods had similar performance in most cases. WFABC also supported the moth data’s nearly fully recessive lethal explanation, but diffusion-based likelihoods and all methods’ demographic assumptions constrain their scope.

Discussion

  • The performance gap between the ABC and likelihood methods was very small in most cases, suggesting Fsd' and Fsi' may be close to statistically sufficient.Source: Discussion · p. 19
  • To approximate the Wright-Fisher model with diffusion, likelihood methods assume continuity in state space and time, requiring small s and large Ne; thus they are limited to large effective population sizes.Source: Discussion · p. 19
  • WFABC was slightly less precise in some small-s cases, but its implementation handles complex ascertainment and estimates genome-wide Ne from multi-locus data.Source: Discussion · p. 20
  • Accurate multi-locus Ne estimation increased selection-coefficient precision and reduced overall RMSE versus the Malaspinas et al. (2012) method.The approach estimates genome-wide Ne before per-site selection coefficients.Source: Discussion · p. 20
  • WFABC applied to Panaxia dominula supported a nearly fully recessive lethal model as the best explanation for medionigra allele-frequency variation.The authors report this as consistent with Mathieson and McVean (2013).Source: Discussion · p. 21
  • Using Fs' on post-1950 observations, when the allele frequency first fell below 0.001, yielded 2Ne=927, consistent with previous estimates.Source: Discussion · p. 21
  • WFABC can co-estimate h and accommodate s=-1, but distinguishing lethality from very strong negative selection (s<-0.5) was difficult, while h had very small variance.Source: Discussion · p. 21
  • The methods assume linkage equilibrium and generally omit demographic history; extending time-series inference to demographic parameters remains a challenge.The Mathieson and McVean (2013) method is the exception, integrating structured populations with a lattice model; ABC estimators are more readily extended to complex demographic models.Source: Discussion · pp. 21–22

Data Accessibility

The section notes that WFABC software is available online and identifies the sources of the P. dominula dataset.

Data Accessibility

  • WFABC software used to create simulated data is available from the software page.The passage gives the URL http://jensenlab.epfl.ch/.Source: Data Accessibility · p. 27
  • The P. dominula dataset is from Cook and Jones (1996) and Jones (2000).The passage also says the dataset is provided, but does not specify where.Source: Data Accessibility · p. 27

Tables

The tables compare selection-coefficient error and simulation runtime across methods and settings. For small s at Ne=1000, Malaspinas et al. (2012) has lower RMSE than WFABC; for large s, the relative RMSE advantage varies by s.

Tables

  • For small s at Ne=1000, Malaspinas et al. (2012) has lower RMSE than WFABC at all five tested values, including 0.011 vs 0.018 at s=0.02.The table reports bias as well as RMSE for each method and s value.Source: Table 1
  • At Ne=1000, WFABC RMSE is lower for s=0.3 (0.046 vs 0.084) and 0.4 (0.059 vs 0.13); Malaspinas et al. (2012) is lower at s=0, 0.1, and 0.2.These values are from the large-s Wright-Fisher diploid scenario.Source: Table 2 · p. 28
  • WFABC runtime is lower in every matched Ne and s setting, including 300 versus 80,000 seconds for small s at Ne=5000.The times are average seconds per simulation replicate.Source: Table 3 · p. 28
  • The Mathieson and McVean (2013) method constantly overestimates small s in Figure 1.Figure 1 presents boxplots of estimated small selection coefficients.Source: Figure 1
Loading 1409.8529v1…