(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



Isolating selective from non-selective forces using site frequency ratios [1]

['Jody Hey', 'Department Of Biology', 'Temple University', 'Philadelphia', 'Pennsylvania', 'United States Of America', 'Vitor A. C. Pavinato']

Date: 2025-05 

Abstract A new method is introduced for estimating the distribution of mutation fitness effects using site frequency spectra. Unlike previous methods, which make assumptions about non-selective factors, or that try to incorporate such factors into the underlying model, this new method mostly avoids non-selective effects by working with the ratios of counts of selected sites to neutral sites. An expression for the likelihood of a set of selected/neutral ratios is found by treating the ratio of two Poisson random variables as the ratio of two gaussian random variables. This approach also avoids the need to estimate the relative mutation rates of selected and neutral sites. Simulations over a wide range of demographic models, with linked selection effects show that the new SFRatios method performs well for statistical tests of selection, and it performs well for estimating the distribution of selection effects. Performance was better with weak selection models and for expansion and structured demographic models than for bottleneck models. Applications to two populations of Drosophila melanogaster reveal clear but very weak selection on synonymous sites. For nonsynonymous sites, selection was found to be consistent with previous estimates and stronger for an African population than for one from North Carolina.

Author summary A new statistical method is presented for estimating the distribution of strengths of natural selection acting on mutations in natural populations using the distribution of polymorphic site allele frequencies. In order to isolate the impact of selection, separately from other demographic and genomic factors that can shape allele frequencies, our method uses the ratio of the frequency of candidate selected variants to the ratio of the frequency of neutral variants of the same frequency. An expression for the overall likelihood across the range of frequency ratios is developed using a gaussian approximation. Testing of the method, called SFRatios, finds that it performs reasonably well across a range of strengths of selection and demographic histories. Applications to two Drosophila populations find estimates of the strength of selection on nonsynonymous coding variants consistent with previous estimates and estimates for synonymous variation quite close to selectively neutral.

Citation: Hey J, Pavinato VAC (2025) Isolating selective from non-selective forces using site frequency ratios. PLoS Genet 21(4): e1011427. https://doi.org/10.1371/journal.pgen.1011427 Editor: Kirk E. Lohmueller,, University of California Los Angeles, UNITED STATES OF AMERICA Received: September 12, 2024; Accepted: March 24, 2025; Published: April 21, 2025 Copyright: © 2025 Hey. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The SFRatios program, along with simulation scripts, as well as script for building the Drosophila data sets, and other scripts used in the analysis of the site frequency spectra from the Drosophila populations, are available at https://github.com/jodyhey/SF_Ratios. The pipeline and scripts for making the Drosophila site frequency spectra are available at https://github.com/jodyhey/SFRatios. Funding: This research was supported by National Institutes of Health (NIH.gov) research grant R01GM144468 to JH. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Population genomics is often a science of sifting signal from noise as investigators regularly seek to distill the signs of natural selection from the confusing patterns of variation that arises from other factors [1–4]. These other factors are quite diverse with some being especially noise-like (genetic drift, recombination and gene-conversion events, and the effects of linked random mutations) and others that have a directional (i.e., non-random) component, as occurs when the demographic history departs from the assumptions of the investigator’s model. Here we describe a new approach for estimating the distribution of selection coefficients acting on mutations, but that does so while largely sidestepping the confounding effects of all these other factors. For questions about selective effects investigators have often employed a classic body of theory on the distribution of allele frequencies [5,6], also known as the site frequency spectrum, or SFS. An important theoretical advance was the realization that the count of observed sites with an allele at a particular frequency could be modelled as a Poisson random variable, with an expected value that depends on a particular model of directional selection and population demography [7,8]. These Poisson random field (PRF) models provide accessible likelihood formulae, not only for the estimation of single selection coefficients, but also for the estimation of the distribution of selection coefficients [9–11]. The original PRF work was limited to constant size Wright-Fisher (WF) populations. To allow for departures from WF models, these methods have been adapted for joint estimates of selection and demography under models of population size change [12–15] and with gene flow between subpopulations [16]. Nevertheless, real populations can have histories that vary in many ways not accounted for with these methods. Nor do such methods account for the many other non-selective non-demographic factors that can shape the distribution of allele frequencies. For example, if rates of gene conversion and gene conversion bias are high enough, they will alter the site frequency spectrum, as will variation in mutation rates if some sites have a high enough mutation to result in some polymorphic sites being caused by multiple mutations. Variation in recombination rates can also affect allele frequencies by shaping the degree to which some sites, more than others, are affected by selection on linked sites. Finally, the structure of the sample across subpopulations can have a very large effect on the site frequency spectrum, one that may easily not be appreciated if there are unknown subdivisions within the sampled population(s). One way to improve upon methods that attempt to jointly estimate selection and other factors is to include in the analysis a set of neutral control variants that are thought to be affected only by non-selective factors. In particular, Eyre-Walker and colleagues developed a method that uses the joint likelihood of selected and neutral variants, and represents shared, non-selective factors by a series of nuisance parameters, one for each frequency bin [10,17]. This approach has been adopted in a number of studies and applications [18–21]. However the method of using both selective and neutral sites, while potentially solving one problem, introduces another complication, which is that the respective mutation rates for each of the two classes must be estimated. This can be done using a previously estimated mutation rate and by assuming a particular demographic history, or by jointly estimating that history. However in most applications it is handled by including in the likelihood function factors and , the number of sampled sites where a neutral or selected mutation could occur, respectively. If these values are known without appreciable error, then the counts of sites that are invariant with respect to a sister species can be obtained, and these can be used to estimate the two mutation rates. One important benefit of this approach is that the divergence measures can be used in turn to estimate the rate of adaptive substitution [9,17,20,21]. For methods that do not include neutral controls or divergence between sister species, and rely only upon the frequency distribution of polymorphisms, the actual mutation rate is not a parameter of much interest. In these cases, it is the changing relative height of the polymorphism count across frequency bins that informs on the effects of selection. However, the methods that use both selected and neutral variants all depend on knowing or successfully estimating the underlying mutation rate, which often depends on knowing the value of the number of sampled neutral positions, . This value will typically include a large number of invariant positions. However if a subset of these is not variable because they are actually under selective constraint, then value will be too large. Another issue that arises when using values as a means to include divergence in the analyses, is that non-selective factors may have changed over the course of the divergence process [10]. Here we take a new approach to using a neutral control set for isolating the effects of direct selection on a set of variants. But unlike other methods that depend on estimates of the overall mutation rates to selective and neutral variants, our method does not depend on estimates of the mutation rates for each class of variant. The method depends not at all on estimates of species divergence or estimates of the number of sampled selected and neutral positions.

Discussion The promise of the SFRatios method is to enable the estimation of selection intensity without having to consider either divergence between species, or the underlying mutation rates, or the other non-selective factors that will also have shaped polymorphism patterns. Given the simulation results, the method works well across a wide array of demographic scenarios, particularly when selection is not strong. For strong selection (e.g., mean most simulation results revealed an estimator bias towards weaker (less negative) selection. Although not examined here in depth, the approach of using ratios should also be effective for other non-selective factors, in addition to demography, so long as they are shared by the selected and neutral variants. For example, if the two sets of variants are sampled near each other (as in this study), then linked selection effects, including background selection and selective sweeps, will have affected both in similar ways. Similarly, the use of ratios should accommodate factors that affect mutation biases and gene conversion biases if they are shared by selected and neutral variants, as is partly the case in this study. However, if the selected and neutral sets differ because of factors that are not shared, such as differing levels of biased gene conversion due to base composition differences, as occurs in mammals [30], then the ratio-based estimate may suffer. If investigators have ready access to counts of invariant sites, that otherwise match the criteria for sampling selected and neutral polymorphic sites, they have other tools available, including fastDFE [19]. However counting invariant sites presents a different set of challenges than counting polymorphic sites. For example, it may require assuming some particular fraction of the genome is susceptible to the class of mutations being studied (as when working with synonymous and nonsynonymous variants). A larger difficulty arises when the sampling effort of invariant sites is unknown or uneven, such as when working with VCF files of polymorphic sites, or when working with pooled samples with differing or unknown sampling efforts. These issues are compounded in their difficulty if counts of invariant sites are to be divided into those that are fixed for an ancestral allele and those that are fixed for a derived allele, as required by many methods. The SFRatios method depends strongly on the quality of the neutral control set of SNPs under three main criteria. The first is selective neutrality, which can be inferred in relative terms using site frequency analyses and divergence patterns [26], but it is difficult to know with certainty. The second criterion is that the neutral set share as many of the non-selective factors as possible with the selected set of SNPs. All parts of the genome share a demographic history, but not all share equally in background selection or mutation and recombination related processes. This second criterion can often be met, at least approximately, by having the neutral SNPs be near to, and interspersed among, the selected SNPs. The third criterion is to avoid increasing the variance of ratios by being sure to record the frequencies of both types of SNPs, selected and control, on the same set of genomes. This criterion arises from the very large effect that unforeseen population structure can have on the SFS of a sample. Consider for example an SFS for 10 genomes sample from a population that, unbeknownst to the investigator, actually includes two divergent subpopulations. In this case the partition of the sample (i.e., the numbers of individuals in each subpopulation) will have a very large effect on the SFS. If the neutral and selected SNPs were drawn from different individuals, then there will be a strong chance that the two sets will have a different partition with respect to the subpopulations, with very large and differing effects on their respective SFSs. In comparison with fastDFE over a wide range of demographies and gamma DFEs, we observed that both SFRatios and fastDFE performed inconsistently, with each method performing better under some of the models. These results highlight the challenge of developing a general estimator that can perform well across a very wide range of DFEs and demographies. It is possible some of the challenge may be computational, particularly if discontinuities in the likelihood surface arise for some parts of the parameter space as a byproduct of approximations used in numerical integration. We did observe that SFRatios performed better for lognormal DFEs (compare Fig 4 with S4 and S5 Figs), which is noteworthy as all of the best-fitting models for the Drosophila data sets were lognormal DFEs (Table 2). Our study of the Zambia Drosophila population can be compared to estimates from other methods that use site frequency data, but that also relied upon mutation rate or divergence measures. Table 3 shows the best fitting model results of three studies, all of which compared lognormal and gamma densities for . These estimates bracket that obtained here that had a mean estimate of -1643.8. Huber et al., [31] estimated a mean of -738.26 in the Zambia sample, while Ragsdale et al., obtained an estimated mean of -2760 on the Zambian sample [32] and another of -7414 on a Rwandan sample [33]. All these other studies rely upon estimates of the number of invariant sites as well as fitting of a demographic model. In our method, by using the ratio of two SFSs, one being the SFS of putative neutral sites, we avoid both of these complications. PPT PowerPoint slide

PNG larger image

TIFF original image Download: Table 3. Previous DFE estimates for nonsynonymous mutations. https://doi.org/10.1371/journal.pgen.1011427.t003

Materials and methods Simulations With SLiM3 [34], we simulated a series of diverse demographies with linked selection effects under both a fixed and a continuous distribution of values, all without dominance. For the distribution, we used an inverted lognormal distribution with a maximum set to 1 and a minimum of -100000 and an inverted gamma distribution with a maximum of 0 and a minimum of -100000. For each combination of distribution and demographic model we ran 20 simulations each with a total genome length of 4 megabase pairs (Mbp) for each of five lognormal distributions that ranged from very weak selection (mean ) to fairly strong selection (mean ) (see Fig 4). To achieve faster running times, we split each genome into 400 fragments of 10 kilobase pairs (Kbp), allowing us to simulate these smaller fragments faster in independent sub-simulations. Each of these 10 Kbp genomic fragments was designed to model a typical Drosophila gene, including introns and flanking sequences. Each fragment included eight consecutive neutral/selected pairs of size 810 + 324 = 1134 bp, and a neutral segments of 928 bases. Selected fragments had only non-neutral mutations (i.e., , while neutral fragments carried only neutral mutations (i.e., ). Base diploid population size ( was set to 1000, with recombination and mutation rates per base pair of , giving population level rates ( ) similar to that seen in human populations (i.e., ). For each simulation, the site frequency spectra (SFS) for all neutral variants sampled from the sub-simulations were summed (as were those for the selected variants), then ratios calculated as the selected count for a frequency bin, divided by the neutral count for the corresponding frequency bin. All simulations began with a burn-in period of generations to allow the population to reach mutation-selection-drift equilibrium [35]. For simulations with changing population size, values were kept constant by rescaling selection coefficients at each generation in the simulation at which population size changed. We chose this approach, rather than the alternative of keeping constant, to be consistent with our model in which (or a distribution of ) is fixed. We considered the following demographic models. (1) Constant population size at . (2). Population expansion, with an initial population of 1000 jumping to 10000, followed by 100 generations before sampling. Population bottleneck, with an initial population of 1000 jumping to 100, followed by 100 generations before sampling. Population structure, with an initial population of 1000, splitting into two populations each of 1000, followed by 1000 generations before sampling equally from both populations. In addition we simulated the human African-Origin (AO) model as inferred by Gravel et al. [24]. This model has some of the features present in the described models, but also some more complex ones. This model includes multiple populations, a bottleneck at the founding of non-African populations, expansion of European and East Asians subpopulations, bottleneck, and population structure as well gene flow among populations. Because the OA simulations with the inferred population sizes are too computationally expensive, we ran a set of neutral simulations to identify the best scaling factor for the simulation parameters that could recover the site-frequency spectra of the original model without any scaling. We then scaled the simulation parameters: mutation and recombination rates, population sizes, migration rates, and growth rates of the exponential population growth phase with a factor of 10. We did not attempt to simulate whole genomes, but 400 sub-simulations of fragments of size 10 Kbp as was done for the other models. For all models described above, we sampled individuals at the end of the simulation. For the AO model, we sampled at the end by taking an equal proportion of the three populations: Africans, Europeans, and East Asians. For all demographic models and lognormal densities, selected and neutral SFSs were generated for samples sizes of both 200 chromosomes and 50 chromosomes to assess the effect of sample sizes. All SFSs were folded before analysis. For distributions of the population selection coefficient we considered a normal distribution, as well as inverted lognormal and gamma distributions, i.e., extending to rather than . Rather than have the upper limit at zero, which would not allow for the inclusion of strictly neutral mutations, we set the upper limit at 1, to include the possibility of weakly advantageous mutations as well as neutral mutations. These densities are then: (10) for an inverted lognormal density with maximum , expectation and standard deviation for the natural logarithm of , and (11) for an inverted gamma density with maximum , mean and shape parameter . We also considered a normal (gaussian) distribution. as well as a mixed distributions, that included a continuous distribution (lognormal, gamma or normal, as described) and a point mass at zero. For any continuous density , the corresponding mixture with a point mass at zero is , where is the Dirac delta function. fastDFE applications The fastDFE program [19] was used to analyze SLiM simulated data sets generated under several inverted gamma distributions with a maximum at zero. fastDFE can run on folded SFSs, but also requires the count of invariant sites, which were obtained from the SLiM runs. We used the GammaExpParametrization model in fastDFE with maximum set to zero and divided the estimated mean by 2, as fastDFE is based on a parameterization of selection strength, in contrast to as used by SFRatios. Drosophila applications We extracted synonymous, nonsynonymous, and short intron site-frequency spectrums from whole-genome sequencing data sets of two Drosophila melanogaster populations for the four autosomal chromosome arms. The North Carolina population has 200 inbred lines [36] while the Zambia collection is based on 197 haploid embryos [37]. Because both data sets have missing data in some lines at some positions, allele counts at all SNPs were down sampled to the expected SFS with a uniform sample size of 160. Because of uncertainty as to which allele is truly ancestral for a given SNP, we used folded SFSs [8]. All sequence data were downloaded from the Drosophila Genome Nexus (DGN, https://www.johnpool.net/genomes.html). For each population a VCF file was constructed by first running the ‘masking package’ (available at DGN) to mask identical-by-descent or admixture tracks when present in one or many genomes, followed by the ‘snp-site’ program [38] to convert the population multi-alignment FASTA to a VCF. Because snp-site assigns the common allele as the reference, a custom script was used to assign the correct reference base. This script also ensures the genotype data conform to the standard VCF format v.4. An additional filter was applied to only keep bi-allelic SNPs genotyped on 50% or more of individuals in each population. With the DGN data in VCF format, we used GATK LiftoverVcf [39] to shift the SNP coordinate positions from D. melanogaster reference genome Dmel 3 to Dmel 6 [40]. We annotated SNPs with SNPEff [41]. For the neutral SFS we used SNPs found in short introns (< 86 bp in length), after removing 8 bp from each side [25,26,42]. To help ensure that the selected and neutral SNPs were as closely matched as possible for local mutational context, we followed Machado et al., [18] by pairing each candidate selective SNP with the nearest short intron (SI) SNP that had the same reference allele and flanking bases. For each candidate selected SNP, a short intron partner was identified by matching a text string that included two nucleotides from the reference genome sequence (1 bp before and after the SNP) and the SNP genotype reference/alternate allele pair (e.g., C/T). As the order of the reference and alternative allele did not matter (all analyses were carried with folded site-frequency spectra) we consider both allele pairs (e.g., C/T and T/C) together when defining the SNP mutational context. For example, if one SNP was C/T, where C was the reference allele, and T was the alternative allele (as in the VCF), and a second SNP was T/C, where T was the reference and C was the alternative, and if they both had an A nucleotide one bp before and after, both SNPs share the same mutational context (e.g., AC/TA or AT/CA). These mutational contexts were defined for every combination of SNP genotypes and flanking bases, giving a list of possible 96 mutational contexts. For each candidate selected SNP, the nearest short intron SNP with matching mutational context, and not previously sampled, was added to the short intron data. Of the three classes of sites nonsynonymous, synonymous, and short intron, the latter class had the fewest SNPs and was the limiting factor for the total number of SNPs included in a data set. We obtained 44,266 and 111,161 short intron and synonymous pairs and 47,208 and 120,777 short intron and nonsynonymous pairs of SNPs with at least 160 genomes genotypes for North Carolina and Zambia, respectively.

Acknowledgments We are grateful to Adam Eyre-Walker for helpful comments, and to Janek Sendrowski for assistance with fastDFE.

[END]
---
[1] Url: https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1011427

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/