Model for triallelic loci

The diffusion approximation we used is based on a triallelic extension to the standard Wright–Fisher (WF) model for allele frequency dynamics, which assumes nonoverlapping generations and random mating. The two derived alleles have selection coefficients, and so their fitnesses relative to the ancestral allele are and If the two derived alleles have frequencies in generation t in a diploid population of size N, then their frequencies in generation are sampled from a trinomial distribution, such that the probability of sampling is(1)whereand is the trinomial coefficient From here on, we focus on relative allele frequencies and

Most applications of the biallelic WF model assume infinite sites, so each new mutation is unique, and new mutations enter the population at a rate proportional to Here is the population-scaled mutation rate, is the ancestral effective population size, and μ is the per-generation mutation rate. Mutations begin at frequency and are assumed to evolve independently. Given these assumptions, the density function for derived allele frequencies in a population can be approximated by diffusion theory (Kimura 1964), such that the expected total number of alleles with frequency between and is a key result from Poisson random field theory (Sawyer and Hartl 1992). The expected sample allele frequency spectrum F with n samples is then(2)where is the binomial coefficient. The likelihood of an observed allele frequency spectrum under this model is then a product of Poisson likelihoods for each entry in the spectrum (Sawyer and Hartl 1992).

Whereas new biallelic mutations begin at frequency triallelic loci are created when a novel mutation occurs at a locus that is already biallelic. The new derived allele initially has frequency and the existing derived allele has a frequency drawn from the population distribution of biallelic frequencies in that generation. The net rate at which triallelic loci arise is thus(3)where is the rate for mutations that hit existing biallelic sites and produce a third allele. Triallelic sites then evolve under the three-locus WF model, and we denote the density function for frequencies of triallelic loci as The triallelic frequency spectrum summarizes sequence data from a sample of individuals by storing the counts of triallelic loci with each set of observed derived allele frequencies (Jenkins et al. 2014) (Figure 1, E and F). The expected triallelic frequency spectrum T with n samples is proportional to the integral of the density function φ against the trinomial sampling distribution:(4)Because the net triallelic mutation rate is sensitive to mutation rate heterogeneity, in our triallelic analyses we focused on the normalized triallelic frequency spectrum, which does not depend on the overall rate of creation. Similarly, because the order in which the two derived alleles arose is often unknown, we considered only counts of major and minor derived alleles, which have respectively higher or lower sample frequencies (Figure 1). That is, for given major and minor derived allele frequencies i and j, with we collapsed the and counts together into the bin. If in a sample we observe counts of independent triallelic frequencies Poisson Random Field theory shows that the data are Poisson distributed with mean enabling likelihood calculations.

Diffusion approximation to the triallelic frequency spectrum with selection

To obtain the expected sample frequency spectrum for a given model of selection and demography, we numerically solved the corresponding diffusion equation. First described by Kimura (1955, 1956), the triallelic diffusion equation models the evolution of the density function for the expected number of loci in the population with derived allele frequencies such that and (Figure 1B):(5)Time τ is measured is units of generations, where is the ancestral effective population size. The spatial second-derivative terms account for genetic drift, which is scaled by the relative population size and the mixed derivative term accounts for the covariance in allele frequency changes. The population-scaled selection coefficient is where s is the relative fitness of the derived vs. ancestral allele. Here that selection coefficient must be adjusted to to account for competition between the two segregating derived alleles, dependent on their allele frequencies. For example, if their selection coefficients are roughly equal, they will be effectively neutral when at high frequency. In general,(6)with a similar expression for

Like the biallelic diffusion method Equation 5 does not account for recurrent mutation, which would tend to increase derived allele frequencies. Recurrent mutation could be accounted for in the first-derivative terms, but at the cost of additional model complexity. If it is common, neglecting recurrent mutation can bias inferences of mutation rate, population size, and selection (Desai and Plotkin 2008; Mathew et al. 2013). Applying our present theory thus requires that the mutation rate be high enough to create a substantial number of triallelic sites for inference, but not so high that a large fraction of biallelic or triallelic sites are affected by recurrent mutation. For most eukaryotes, including humans and Drosophila, mutation rates are low enough that recurrent mutation is negligible in most applications (Desai and Plotkin 2008).

Some analytic results are known for triallelic diffusion (Tier and Keller 1978; Tier 1979; Spencer and Barakat 1992), but we solved Equation 5 numerically. We used a finite-difference method similar to that in (Gutenkunst et al. 2009). To integrate the diffusion equation forward in time, we used operator splitting to separately apply the nonmixed and mixed derivative terms each time step (Supplemental Material, File S1). We integrated the nonmixed terms, using a conservative alternating direction implicit (ADI) finite difference scheme (Chang and Cooper 1970). We integrated the mixed term, using a standard explicit scheme for mixed derivatives. We used uniform grids in x and y with equal grid spacing so that grid points lie directly on the diagonal boundary of the domain, which readily allowed the diagonal boundary to be absorbing. Although these integration schemes worked well in the interior of the domain, application at the diagonal boundary led to an excess of density being lost (File S1 and Figure S1). To avoid this excess loss, we did not apply the ADI and mixed derivative schemes at the closest grid points to the diagonal boundary. Instead, at each time step we calculated the amount of density at each grid point that would fix along the diagonal boundary, and we directly removed that amount from the numerical density function and added it to the boundary.

To inject density into φ for new triallelic loci, at each time step we added density to the first interior rows of grid points based on the expected background biallelic frequency For example, we added to the row of grid points with weight for point proportional to the biallelic population allele density at frequency x. We directly coupled with to track To obtain the expected sample frequency spectrum T from the population frequency spectrum φ, we numerically integrated against the trinomial distribution with sample size n, using Equation 4. Our code implementing these methods is integrated into available at https://bitbucket.org/gutenkunstlab/dadi.

Calculating frequency spectra under a DFE

Given a DFE, the expected sample frequency spectrum can be obtained by integrating over the expected frequency spectrum for each selection coefficient, weighted by the DFE. For biallelic sites, the DFE is a univariate distribution. For triallelic sites, the DFE is a two-dimensional joint distribution, because there are two derived alleles. Moreover, the two marginal distributions are identical, because we assume no knowledge of which allele arose first.

For our primary analysis, we used a lognormal model for the deleterious triallelic DFE (Figure 1, C and D), plus a point mass of positive selection. The lognormal distribution readily generalizes to an arbitrary number of dimensions, and the bivariate lognormal distribution includes a correlation coefficient ρ that characterizes the correlation between selection coefficients. If the selection coefficients of the two derived alleles at a single triallelic locus are independent, whereas if they are equal. For a fixed marginal DFE, as the correlation coefficient ρ increases, more segregating triallelic loci are expected, particularly at moderate and high derived allele frequencies (Figure 1, C–F). We quantified the relative importance of identity and location for protein mutation fitness effects through ρ; low correlation suggests that identity is more important, whereas high correlation suggests that location within the protein is more important.

To numerically integrate over the univariate DFE, we used a logarithmically spaced grid with 2000 grid points ranging from to along with and a point mass of positive selection Biallelic spectra were cached for each resulting in 2001 cached spectra. We assumed that alleles with were effectively lethal and did not contribute to the sample frequency spectrum. We also assumed that alleles with were effectively neutral, and we used the cached spectrum for for contributions from this range of the DFE (Figure S2A).

To integrate over the bivariate DFE we used a logarithmically spaced grid with 50 grid points ranging from to along with and determined by the univariate DFE fit. We cached spectra for each possible pair yielding cached spectra. A pair of selection coefficients could fall into four quadrants, depending on the sign of and The overall frequency spectrum was calculated by summing over the weighted frequency spectra for each quadrant based on the DFE parameters and ρ. The weights were for both for one selection coefficient positive and the other negative, and for both These weights were found by taking the distribution of two point masses (one for positive selection, and one for negative selection, ) and extending it to a bivariate distribution of point masses with correlation coefficient ρ (File S1). To integrate over the continuous distributions with one or both of the selection coefficients negative, we used the trapezoid rule. We approximated as effectively neutral and as effectively lethal (Figure S2B).

Genomic data

We extracted SNPs from phase 3 of the Drosophila Population Genomics Project (DPGP3) population of fruit flies from the Drosophila Genome Nexus Data (Lack et al. 2015). The data we used consist of 197 sequenced genomes from a Zambian population obtained through high-coverage haploid embryo sequencing. This population has high genetic diversity, and it did not experience the out-of-Africa bottleneck or New World admixture that other D. melanogaster populations have experienced (Lack et al. 2015). We used Annovar (Wang et al. 2010) to determine the transcript and codon position of each coding SNP. The ancestral state of each codon was determined using the aligned sequences of D. melanogaster (April 2006, dm3) and D. simulans (droSim1) downloaded from the University of California, Santa Cruz genome database, by assuming that the D. simulans allele was ancestral. We excluded loci with no aligned D. simulans sequence. We downloaded the reference transcript sequences from Ensembl Biomart (Flicek et al. 2014) and used the ancestral states determined by the droSim1 alignment to determine the ancestral codon state.

Inferring the selection correlation coefficient

In our application to D. melanogaster, we used biallelic synonymous data to infer the single-population demographic history and then used nonsynonymous data to infer the parameters of the DFE. Using the unfolded synonymous allele frequency spectrum, we fitted a neutral three-epoch demographic model. This model has two instantaneous size changes, at times and in the past, with constant population sizes, and relative to the ancestral population size. We also included a parameter to account for ancestral state misidentification, which creates an excess of high-frequency derived alleles (Baudry and Depaulis 2003). Specifically, we compared the data not with the expected true unfolded frequency spectrum under the demographic model, but rather with the expected observed unfolded frequency such that where n is the sample size. We chose to include misidentification in our model rather than adjusting the data spectra (Hernandez et al. 2007), because adjusting the data leads to violations of the Poisson random field assumption, most obviously when the adjustment leads to negative entries in the data spectrum. The population-scaled mutation rate was an implicit free parameter. We used the built-in optimization routines in (Gutenkunst et al. 2009) to fit the model to the data. We fixed this demographic model for all future inferences.

The unfolded biallelic nonsynonymous allele frequency spectrum was used to infer the marginal DFE. As described above, we used a lognormal distribution for negative selection combined with a point mass of positive selection. This yielded a total of four parameters, μ and σ for the lognormal portion and and proportion for the point mass. As in the fits for demography using synonymous data, we also included a parameter to model ancestral state misidentification. In this fit, the population-scaled mutation rate was fixed to and we again used ’s optimization routines to fit the DFE to the data.

Finally, we used triallelic data with two mutually nonsynonymous derived codons to infer the correlation coefficient ρ. We fixed the demography to that inferred from the biallelic synonymous data, and we fixed the DFE parameters μ, σ, and to the values inferred from the biallelic nonsynonymous data. This left the correlation coefficient ρ as the only free parameter of the bivariate DFE, and we also included a free parameter to account for ancestral misidentification. Assuming that the two observed derived alleles were equally likely to be the true ancestral allele, we calculated the expected observed triallelic spectrum from the expected true spectrum by We also left the overall population-scaled mutation rate for triallelic loci as an implicit free parameter, so our fit considered only the distribution of triallelic codons among frequency classes, not the overall number of such codons. We did this because the overall number of triallelic codons can be strongly affected by mutation rate heterogeneity, and imperfect modeling of that heterogeneity could bias our results.

We estimated model parameters by maximum composite likelihood. Following the Poisson random field framework, likelihoods of the data D given the model parameters were calculated by assuming that each entry in the observed triallelic frequency spectrum was an independent Poisson random variable with mean (Sawyer and Hartl 1992), where T is the expected triallelic frequency spectrum generated under (7)Because our SNP data are not actually independent, is not the true likelihood, but rather a composite likelihood. To account for this, we calculated parameter uncertainties for each model fit, using the Godambe information matrix (Coffman et al. 2016), which adjusts the composite-likelihood statistic to account for the effects of linkage. To do so, we generated 1000 bootstrap data sets by dividing the D. melanogaster autosomal genome into 1000 regions of equal length and resampling among these regions.

Tests on simulated data

To generate simulated data for tests of statistical power, we first calculated the expected frequency spectrum under each model considered, using our diffusion method. To generate an observed frequency spectrum with exactly n entries, we generated n multinomial samples of frequencies, weighted by the expected frequency spectrum. To generate an observed frequency spectrum with a given mutation rate θ, we scaled the expected frequency spectrum by θ, treated the bin weights as Poisson random variables, and sampled independently for each bin.

Mutational scanning data

For comparison with our population genetic inference, we considered data from three mutational scanning studies (Roscoe et al. 2013; Firnberg et al. 2014; Starita et al. 2015). Each study assayed a different protein from a different organism, using a different proxy for fitness. In all three experiments, the distribution of fitnesses was bimodal, with peaks of moderately and strongly deleterious mutations, although the relative sizes of these peaks differed markedly between experiments (Figure S3, A–C). To calculate the fitness correlation coefficient, we sampled a pair of mutually nonsynonymous mutations from each site in the protein (excluding mutations without reported fitness) and calculated the Pearson correlation of those fitnesses. The confidence intervals in Table 1 are 2.5% and 97.5% quantiles from 10,000 repetitions of this sampling. To visualize the correlations, we calculated the proportion of mutually nonsynonymous mutation pairs within each possible bin of joint fitness effects (Figure 4B and Figure S3, D–I). Because our population-genetic analysis is not sensitive to strongly deleterious mutations, we focused our analysis on moderately deleterious mutations (shaded regions in Figure S3, A–C, joint distributions in Figure S3, D–F). For details on each data set, see File S1.

Data availability

The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article.