Estimating the Strength of Selective Sweeps from Deep Population Diversity Data
key words: selective sweep, deep diversity data, selection coefficient
running title: Estimating the Strength of Selective Sweeps
MPI for Developmental Biology
Tel: +49 (0)7071 601 1345
Selective sweeps are typically associated with a local reduction of genetic diversity around the adaptive site. However, selective sweeps can also quickly carry neutral mutations to observable population frequencies if they arise early in a sweep and hitchhike with the adaptive allele. We show that the interplay between mutation and exponential amplification through hitchhiking results in a characteristic frequency spectrum of the resulting novel haplotype variation that depends only on the ratio of the mutation rate and the selection coefficient of the sweep. Based on this result, we develop an estimator for the selection coefficient driving a sweep. Since this estimator utilizes the novel variation arising from mutations during a sweep, it does not rely on preexisting variation and can also be applied to loci that lack recombination. Compared with standard approaches that infer selection coefficients from the size of dips in genetic diversity around the adaptive site, our estimator requires much shorter sequences but sampled at high population depth in order to capture low-frequency variants; given such data, it consistently outperforms standard approaches. We investigate analytically and numerically how the accuracy of our estimator is affected by the decay of the sweep pattern over time as a consequence of random genetic drift and discuss potential effects of recombination, soft sweeps, and demography. As an example for its use, we apply our estimator to deep sequencing data from HIV populations.
The frequency and strength of positive selection are key parameters of the evolutionary process yet reliable estimates are often very difficult to obtain (Nielsen et al., 2005; Eyre-Walker, 2006). As it becomes increasingly practicable to analyze the genetic diversity of natural and experimental populations at high depth, we can hope to obtain better estimates from a detailed analysis of the signatures positive selection is expected to leave in such data.
The standard model describing the population genetic signatures of positive selection is the so-called selective sweep (Maynard Smith and Haigh, 1974). Selective sweeps produce distinct effects on the patterns of genetic diversity in the vicinity of the adaptive site: (i) Hitchhiking of linked neutral polymorphism with the adaptive allele generates dips in diversity that level off with increasing genetic distance from the selected site (Maynard Smith and Haigh, 1974). (ii) Because different lineages from a population sample taken after the sweep can coalesce almost instantaneously during the early phase of the sweep, there should be an excess of singletons in the neutral diversity around the selected site (Slatkin and Hudson, 1991; Barton, 1998; Durrett and Schweinsberg, 2004). (iii) Derived alleles can be brought to high frequencies (Fay and Wu, 2000), which is unlikely under random genetic drift alone. (iv) The adaptive haplotype should be longer than expected under neutrality since recombination has had only a short time to degrade it (Hudson et al., 1994). These hallmark signatures underlie the majority of approaches to scan for recent adaptation in population genetic data (Hudson et al., 1987; Tajima, 1989; Wiehe and Stephan, 1993; Fay and Wu, 2000; Sabeti et al., 2002; Kim and Stephan, 2002; Przeworski, 2003; Nielsen et al., 2005; Voight et al., 2006; Sabeti et al., 2007; Andolfatto, 2007).
In addition to detecting selective sweeps, one would often like to know the strength of selection that drove the sweep (Macpherson et al., 2007; Hernandez et al., 2011; Sattath et al., 2011). One common approach is to infer selection coefficients of the adaptive allele from the size of the dip of approximate length around the adaptive site, where is the rate of recombination per length in the corresponding region (Maynard Smith and Haigh, 1974; Kaplan et al., 1989; Kim and Stephan, 2002). Approaches based on this signature rely on the interplay between recombination and ancestral variation and are thus not applicable for scenarios where recombination or ancestral variation is poorly characterized. Furthermore, certain scenarios of adaptation do not always generate clear dips in diversity. Examples are incomplete sweeps where the adaptive mutation is not fixed in the population, and soft sweeps, where more than one haplotype has swept through the population (Pritchard et al., 2010).
Here, we develop an estimator for selection coefficients at candidate loci where a selective sweep has occurred recently. Our estimator is based on the analysis of the novel haplotype variation that arises from neutral mutations on the sweeping adaptive haplotype. These early variants are very different in kind from the variation due to neutral mutations occurring after a sweep, since they have experienced an initial phase of exponential amplification. Mutations after the sweep will quickly outnumber the few that happened early during the sweep, but they will drift to higher population frequencies comparatively slowly.
The frequency spectrum of the early haplotype variants is determined by the distribution of their seeding times during the sweep, as illustrated in Figure 1. We show analytically that their rank-frequency spectrum (the frequencies ordered by decreasing abundance) decays according to a simple power-law, which differs markedly from the approximately exponential decay expected under neutral evolution. This power-law is characterized by only a single parameter: the ratio of the rate at which new haplotypes arise and the selection coefficient of the sweep.
We use this result to develop an estimator of that compares the strength of selection to the rate at which novel haplotypes are produced. Novel haplotypes can be produced by mutation or recombination. In many organisms, recombination is rare and the mutation rate is larger or similar to the recombination rate (Roach et al., 2010; Ossowski et al., 2010; Haag-Liautard et al., 2007). Hence our estimate is much less sensitive to poorly characterized recombination rates, thereby overcoming several problems of estimators based on dips in diversity. In particular, our estimator should also be applicable to organisms that exchange genomic segments via horizontal transfer (many bacteria), and populations where levels of ancestral variation are very low, for instance in experimental evolution with clonal populations.
Our estimator relies on deep diversity data for accurate measurements of the population frequencies of rare haplotypes. With large-scale sequencing projects such as the 1000 genomes project in humans (1000 Genomes Project Consortium et al., 2010) and similar projects in flies (Drosophila Population Genomics Project, 2011) and plants (Cao et al., 2011) currently under way, such data will become available soon.
Materials and Methods
Forward simulations of selective sweeps
Simulations were performed with custom written C++ programs that represent each haplotype in the population by a bit string of fixed length or . Populations are initialized with a sample of size individuals obtained from the neutral coalescent using the program ms (Hudson, 2002). Note that this initial variation is only needed to allow diversification by recombination. We constrained ms to return a sample with (or , respectively) segregating sites, leaving the site in the center of the locus for the beneficial mutation. Each genotype sampled from the coalescent is copied times into the initial population. The beneficial mutation is introduced into one randomly chosen individual in generation 0. If the beneficial mutation is lost due to random drift, the population is reset to the neutral sample and the beneficial mutation is introduced again until a successful sweep is observed. Note that our initial condition with each site polymorphic is not supposed to imply that the product of the population size and the per site mutation is much larger than one. Instead, those loci correspond to the subset of polymorphic loci scattered along a genomic region. The results do not depend on the number of segregating sites in the initial sample, as long as the pairwise difference between two randomly chosen genomes is almost always much larger than one. We repeated the simulations using only a tenth of the segregating sites and obtained similar results (see Supplement).
The Malthusian fitness of a haplotype is given by or , depending on whether or not the haplotype carries the beneficial allele at position . In each generation, a pool of gametes is produced to which each individual contributes a Poisson distributed number of copies with mean , where serves as a carrying capacity to keep the population size approximately constant at and is the mean fitness in the population. A fraction of the newly produced gametes are paired at random and crossed over at a randomly chosen position (since , there is no crossover in most gametes). At the neutral loci, mutations are introduced at random positions into the gametes with the specified mutation rate . When simulating the sweeps of different strength, the mutation rate is changed accordingly to keep the ratio of the strength of selection and the total mutation rate in the range specified in the figure. This is akin to changing the size of the locus simulated, and allows more efficient simulations. The code is available on request.
Analysis of HIV data
Sequences from the relevant samples were read by a python script and translated. Sequences corresponding to identical amino acid sequences were grouped together. Within each of these groups, the number of transitions and transversion that do not change the amino-acid sequence were determined. In case the sequence overlapped with the end of the gene and extended into the long terminal repeat (LTR), only the aminoacid sequence up to the stop codon was used to group sequences and mutations in the LTR treated as neutral mutations. Sequences from the study on early immune system escape mutations by Fischer et al. (2010) are available from the NCBI short read archive under accession number SRA020793. Sequences from the study on drug resistance evolution by Hedskog et al. (2010) were obtained from the authors.
Haplotype frequency spectrum of a selective sweep
Consider a new adaptive mutation with selection coefficient in a panmictic haploid population of constant size . We assume that selection is strong (). Once an adaptive mutation becomes established in the population by reaching a population frequency that assures its escape from future stochastic loss (Maynard Smith, 1971), its frequency trajectory can be modeled by logistic growth
Let us first assume that recombination is infrequent compared to mutation (we will discuss recombination below). If neutral mutations occur during the early phase of the sweep on the adaptive haplotype, they can generate new variants of the adaptive haplotype that can also rise in frequency (Figure 1A). We adopt an infinite haplotypes model where every such mutation creates a distinct haplotype. When the frequency of the beneficial allele is still small, novel adaptive haplotypes become established in the population at an approximate rate , where is the rate at which neutral mutations occur on the sweeping haplotype per generation. The factor accounts for the establishment probability of those new haplotypes. After the sweep is completed, a novel haplotype variant will be the more frequent the earlier it arose in the sweep. Thus, to understand the spectrum of haplotype frequencies, we have to study the distribution of times at which these haplotypes are seeded.
Since novel haplotypes are seeded in independent events, the total number of established new haplotypes up to time is Poisson distributed with parameter
The approximation assumes that the relevant haplotypes are seeded while the sweeping allele is still rare and increases exponentially, which is appropriate as long as , and that , which can be ensured by reducing the size of the locus under investigation.
The probability density that haplotype is seeded at time is given by the rate of establishing new haplotypes multiplied by the probability of having haplotypes at time
where we have used the same approximations as in Equation (2) and dropped factors independent of that ensure normalization. This distribution is shown in Figure 1B for different . The intervals between the peaks of and become smaller with increasing , while the width of each peak decreases. More precisely, we derive from Equation (3) that the most likely seeding time of the th haplotype is given by
and that the peak of has a width of (if ). The joint distribution of seeding times as well as the unordered frequency spectrum are derived in Supplementary Information, section 2.
Assuming that the th haplotype establishes at and has same fitness as the founding haplotype, it will also rise in frequency, albeit with a time-lag . Together with Equation (4) we obtain an approximation for the expected frequency of the th haplotype
The limit corresponds to the asymptotic behavior as the beneficial allele approaches fixation. Note that this expression holds for haplotypes . The zero-th haplotype, i.e., the haplotype the mutation inititally arose on, grows as and asymptotes to . Each of the haplotype variants carrying the adaptive allele thus grows at a slightly smaller rate, , than the overall rate at which the adaptive allele rises, accounting for the founding of new haplotype variants by mutation.
The interplay between the exponential amplification of the adaptive allele and the generation of new adaptive haplotype-variants thus gives rise to a simple power-law decay of the rank-frequency spectrum of adaptive haplotypes: the most abundant adaptive haplotype, on average, is typically times more frequent than the 2nd most abundant adaptive haplotype, times more frequent than the 3rd most abundant adaptive haplotype, and so forth. This power-law spectrum with exponent differs markedly from that of neutral evolution, where haplotype frequencies are expected to decay as with haplotype rank (; see Supplementary Information, section 1).
The distribution of haplotype frequencies after a sweep is related to the distribution of “family sizes” studied by Barton (1998). Barton presented numeric results for the distribution of the size of groups of individuals that share a common ancestor at the beginning of the sweep. Haplotype frequencies, however, refer to the size of groups identical by state, rather than descent, and characterize diversification that happened after the sweep, rather than the ancestral variation that survived the sweep.
In order to compare the haplotype counts in a sample of size to Equation (5), we need an estimate of . The simplest estimate of is the sample frequency of the most abundant haplotype, , whose deterministic value would be precisely . However, can be quite variable because the variances of the first few seeding times are still large. To circumvent this problem, one can construct a more stable estimate of based on the first haplotypes and compare the sum of the sample frequencies to the sum of the deterministic expectation of the frequencies:
where . This effective is insensitive to as long as the first few haplotypes are included. Figure 2A shows a comparison of the simulated haplotype spectra with
The logarithmic axes of the plot are chosen such that the algebraic decay shows up as a slope . The observed slope is correct initially, but the curves become slightly flatter after haplotype rank . This is due to degradation of the sweep pattern by random genetic drift, which we will discuss next.
In deriving Equation (7), we have neglected fluctuations in the seeding and establishment times of novel haplotype variants as well as genetic drift after the sweep is complete. Relaxing these assumptions requires a stochastic calculation, which is given in the Supplementary Information, sections 3A and B. Specifically, we calculate the probability of finding haplotypes to be present in more than copies in the population, assuming the logistic trajectory of the beneficial allele given in Equation (1). These calculations are lengthy, but the resulting effects can be understood by the simple arguments given below.
Variations in the establishment times will cause the haplotype spectra to fluctuate and result in spectra below or above Equation (7) if new haplotypes were seeded fortuitously early or late. As expected from the distribution of seeding times in Figure 1, the frequencies of the first few haplotypes fluctuate substantially, while the later ones do not. Since the number of haplotypes seeded by time is Poisson distributed, so is the number of haplotypes above a certain frequency after the sweep is complete (neglecting genetic drift, there is a one to one mapping between establishment time and eventual frequency). Hence, the fluctuations due to random seeding times are small if the expected number of haplotypes above the chosen frequency is large. In accordance with these arguments, variation of haplotype spectra in Figure 2A decreases with increasing .
Once the beneficial allele approaches fixation, the strength of selection goes to zero, and exponential amplification ceases. The frequencies of the rare haplotypes thereafter change primarily due to random drift, while the frequencies of the common haplotypes decrease due to additional mutations that produce new rare haplotypes. Both of these processes lead to a homogenization of haplotype frequencies, i.e., common haplotypes becoming rarer and rare haplotypes becoming more common, which ultimately wipes out the sweep signature. This relaxation to the neutral haplotype spectrum is shown in Figure 2B.
Since the time required to drift to frequency is approximately , the haplotype rank-frequency spectrum will soon develop a bulge of drift dominated haplotypes at low frequency and high rank. This accumulation of rare haplotypes due to genetic drift can be calculated perturbatively (see Supplementary Information, section 3). For the expected number of haplotypes with frequencies above , one obtains
where is the age of haplotype which is of the same order as the age of the sweep. Note that Equation (8) is essentially Equation (5) solved for with the “drift contribution” subtracted from the frequency . Since the age of the haplotype is similar to that of the beneficial allele, , we conclude that for frequencies less than , the sweep spectrum is eroded, while for frequencies much larger that , it prevails. After a time of order , the entire spectrum has relaxed to the neutral haplotype frequency spectrum as shown in Figure 2B. Note that the time it takes for the beneficial allele to go from frequency 0.9 to complete fixation can be long, and substantial erosion of the sweep spectrum can happen during this time interval.
In general, the time that has elapsed since the beginning of the sweep is unknown. To obtain an estimate of the age of the sweep, it is useful to reconsider Equation (5), which states that the frequency of the founding haplotype is expected to asymptote to . This behavior, and in particular the more accurate estimate of given in Equation (6), allows one to obtain a rough estimate of the sweep’s age.
In addition to genetic drift, limited sampling of the population results in a deviation of the observed from the expected sweep spectra. The detected haplotype counts are a convolution of the their true frequencies with the distribution expected from sampling the population. After ranking of haplotype counts, a large number of rare variants leads to a flattening of the spectrum, and in particular an excess of singletons.
Estimating selection strength
According to Equation (7), the expected rank-frequency spectrum of haplotypes in a selective sweep is determined by the single parameter – the ratio of the mutation rate over the locus and the strength of selection. Given an independent estimate of , one can therefore use the haplotype frequency spectrum to estimate the selection coefficient of a sweep, for example by simply counting the number of different haplotypes present above a specified frequency cutoff in the sample. Let be the overall number of different haplotypes which are present in at least copies each. The estimator for the strength of selection is then
where is either set to the observed or determined via Equation (6).
The estimator allows us to either fix and then determine the count from the data, or fix and then measure . In either case, the cut-off should be chosen to achieve maximum accuracy. Since the number of haplotypes above a frequency threshold is Poisson distributed, fluctuations of decrease with smaller and larger . To low , however, will include parts of the frequency spectrum that has already been degraded by drift, which predominantly affects the rare haplotypes. In addition, limited sampling depth limits . Hence should be chosen as large as possible, but such that and the frequency spectrum is still of power-law form down to haplotype . In this case, the relative error of is of order (assuming the model assumptions are met).
It is generally advisable to fix and determine because spectra come flatter than due to the confounding effect of genetic drift and limited sampling of the population. Figure 3 shows the performance of our estimator when applied to simulated sweeps (see Methods) for different selection coefficients and mutation rates as well as its dependence on the choice of . The simulations confirm that our estimator performs accurately over the range of moderate selection () to strong selection () given the parameters used. However, there is a systematic downward bias for small . This bias is due to genetic drift and limited sampling depth.
In the Supplementary Information, section 3, we derive an estimator that accounts for genetic drift perturbatively, see Equation (24) of the supplementary information. In essence, this estimator contains the same correction present in Equation (8), i.e., it subtracts the drift contribution from the frequency of the rare haplotypes. Importantly, this drift correction sets the limits of applicability of : (i) After the completion of the sweep, we can estimate only for a limited time. Specifically, for a given sample frequency , we require that . (ii) The range of that can be reliably estimated by the method is limited: even if we catch the sweep right at the time of fixation of the adaptive allele, it still needs to have occurred faster than the time scale of its own degradation by drift. Given that the time it takes to complete a sweep is approximately and that , we require . Together with , we find that the range of that can be estimated is bounded from below by . This estimate is consistent with the deviations in Figure 3 for .
Up to now, we have neglected any potential effects of recombination on the frequency distribution of the new adaptive haplotype variants arising in a selective sweep. In fact, it is one of the advantages of our estimator that it does not exclusively rely on recombination for its inference of selection coefficients and can thus also be applied to systems which lack recombination or where recombination rates are not well known. However, if recombination occurs on the sweeping haplotype, this can generate new variants of the adaptive haplotype in a manner similar to mutation. This becomes important whenever the recombination rate over the locus is comparable to its mutation rate . In this case, the rate at which new haplotype variants are generated during the sweep will actually be higher than assumed based on the mutation rate , and we thus expect to underestimate the strength of the sweep.
We propose two ways to incorporate recombination: First, one can treat recombination analogously to mutation and assume that every recombination event generates a new variant of the adaptive haplotype that will be different from all other existing variants. Under this assumption, the mutation rate simply needs to be replaced by the sum . Figure 4A shows how our estimator performs for different ratios of mutation to recombination rate under this approach. The infinite haplotypes assumption is appropriate if ancestral diversity is high and independent recombination events are thus unlikely to yield equal outcomes. If, however, ancestral diversity is low, the infinite haplotypes model will overestimate the rate at which new haplotypes arise.
Alternatively, we can try to effectively exclude most recombined adaptive haplotype variants. Recombined variants are likely to differ from the founding haplotype at more than one site, whereas variants originating from mutation should differ at only one site. By filtering out the former, we can treat the remainder as originating primarily from mutation events. As shown in the Supplementary Information, section 4, such a set of pruned “mutation-only” haplotypes has a slightly different rank-frequency spectrum: The exponent is exactly 1, rather than . Besides that, Equation (9) can be applied with being the mutation rate and being the abundance of the haplotype in the pruned set of haplotypes. The estimates obtained by this procedure are shown in Figure 4B. While this modified estimator still seems to underestimate selection coefficients slightly, it performs consistently across different recombination rates.
The difference in the exponent of the rank frequency spectrum arises from slightly different rates of seeding and amplification of new haplotypes. Without pruning, the rate of establishment of novel haplotypes is proportional to the frequency of the beneficial allele, which follows a logistic growth with rate . Haplotype frequencies, however, only grow with rate . This slight difference in rates is responsible for the deviation of the exponent from by or , where the latter also accounts for recombination. When restricting the analysis to haplotypes that descend directly from the founding haplotype, the rate of establishment of novel haplotypes is proportional to the frequency of the founding haplotype. Hence, establishment and amplification happen with the same rate and the exponent is exactly 1.
We note another minor difference in the dynamics of haplotype frequencies in the presence of recombination. When the beneficial allele reaches fixation, the most abundant haplotype will most likely recombine with itself. Thus it diversifies more slowly after the sweep is complete. Minor haplotypes, however, continue to diversify through recombination, which has the effect of slowly decreasing their frequency and producing novel low frequency haplotypes. This slow decrease opposes the effects of genetic drift, which has the tendency to increase the frequency of minor alleles. In simulations, estimates of often become more accurate due to this effect.
Application of to deep sequencing data from HIV
Deep population diversity data of the kind required for the application of our estimator have recently been obtained for HIV populations by sequencing viral RNA from plasma samples (Tsibris et al., 2009; Fischer et al., 2010; Hedskog et al., 2010). We will present three examples from such studies to validate our method and discuss its applicability.
As a first example, we investigated a sample taken shortly after infection when HIV evolution is primarily driven by mutations that allow the virus to escape the immune system. For each of the patients studied in Fischer et al. (2010), samples from several time points were investigated and amplicons sequenced using 454 technology. In one of the epitopes studied (gene nef in patient CH40), a mutation spread sufficiently slowly such that it was possible to observe the mutation rise from low frequency at day 16 to high frequency at day 45 (Figure 5A). From this time series Fischer et al. (2010) estimate a selection coefficient of 0.3, assuming a generation time of HIV between 1.5 and 2 days (Perelson et al., 1996),
To validate our estimator, we can compare this direct estimate of the selection coefficient obtained from time-course data with that obtained from the frequency spectrum of haplotype diversity at the locus, which is shown in Figure 5B. The figure contains the frequency of haplotypes with the dominant amino acid sequence which differ only by putatively neutral synonymous mutations. The haplotype rank-frequency spectrum according to Equation (9) suggests a ratio of . There is considerable uncertainty in the mutation rate estimates for HIV, ranging from (Mansky and Temin, 1995) to per site and generation (O’Neil et al. (2002) for mutations in the LTR). The majority of these mutations are transitions. In our example the sequenced locus tolerates a total of 93 transitions that do not change the amino-acid sequence of Nef or fall into the LTR which overlaps the sequenced region. Using per generation as an averaged estimate for the transition rate, we arrive at
The uncertainty stemming from the random times at which the haplotypes are seeded is expected to result in relative errors , which for used here amounts to . Additional uncertainty (probably larger) in the mutation rate needs to be added to these error bars. Given these uncertainties we consider our estimate in excellent agreement with the independent estimate from time-course data. Note that recombination is not expected to make a large contribution to haplotype diversification since the effective HIV recombination rate is less than half the mutation rate (Neher and Leitner, 2010; Batorsky et al., 2011). Furthermore, the HIV population has undergone several sweeps prior to the date this sample was taken such that the ancestral diversity was very low.
We also checked whether the haplotype variation in our sample is consistent with the characteristic star-phylogeny expected for mutations that arose on the background of the sweeping founder haplotype. In particular, under the assumption of an infinite sites model where each such mutation gives rise to a new haplotype, we would expect that descendants should always be one mutation away from the founding haplotype, while any two descendants should be two mutations away from each other. Our example shows precisely this pattern, as can be seen in the inset of Figure 5A. This is further evidence that recombination did not contribute substantially to haplotype diversification.
As a second example, we consider deep sequencing data of the reverse transcriptase (RT) locus of HIV (120 bp, amino acids 180-220) (Hedskog et al., 2010). The patients analyzed in this study were under anti-retroviral treatment, which implies strong selection pressure for resistance mutations at the RT locus. During the course of treatment resistance has evolved repeatedly in several patients. We therefore expect to see signatures of strong selective sweeps in these viral populations.
As before, we restricted our analysis to synonymous differences because our model assumes selectively neutral mutations. The 120 bp analysis windows typically contained enough variation such that in most samples more than 8 different haplotypes were sampled. Hence we choose for estimating and then measured the cutoff frequency of the least abundant haplotype to evaluate Equation (9).
Figure 6A shows the measured haplotype rank-frequency spectrum for one of the patient samples (time point 4, patient 4). The spectrum decays with the characteristic power-law of a selective sweep. As discussed above, we determined the number of synonymous transitions in the founding sequence and used as the rate of individual transitions. With this mutation rate estimate, we arrive at per generation. Note that the causative mutation, or the combination of mutations, that drove this sweep does not necessarily have to lie inside the analyzed window but could also be located elsewhere in the genome.
Figure 6B shows haplotype data from another population sample (time point 4, patient 5). In this third example, the pairwise mutational distances between different haplotypes reveal a more complex pattern (left inset) that does not seem to be compatible with the simple star phylogeny observed in the example from Figure 6A. However, when we focus on the two most abundant haplotypes 0 and 1, we see that these two haplotypes differ at more than one site from each other. The majority of the remaining haplotypes differs by only one mutation from either haplotype 0 or 1. This pattern can be interpreted as the result of two independent partial sweeps where the two founding haplotypes (0 and 1) differ at two sites.
We can disentangle these two hard sweeps by assigning the rare haplotypes to either of the two founding haplotypes if they from the latter by only one mutation. When reordering the distance matrix according to this assignment, we recover the two individual sweep pattern as blocks on the diagonal. This is shown in the right inset of Figure 6B. Both reordered components now exhibit the characteristic power-law decay in their individual rank-frequency spectra with similar selection coefficients, again in the range of a few percent.
A scenario as the one shown in Figure 6B, with several haplotypes carrying the same adaptive mutation rising in frequency simultaneously, is commonly referred to as a “soft sweep” (Hermisson and Pennings, 2005). Soft sweeps can occur, for example, if an adaptive mutation arises repeatedly on different haplotypes, which is expected to be common in viruses with large population sizes and high mutation rates such as HIV, where every single point mutation can arise many times each generation. Alternatively, soft sweeps can result when a sudden change in environment renders a previously neutral or deleterious allele selectively beneficial and adaptation involves several different adaptive haplotypes from the standing genetic variation.
An unambiguous decomposition of the haplotype distance-matrix into individual sweep components of a soft sweep requires that the founding haplotypes differ by several mutations. This is likely to be the case for a diverse ancestral population. Note that the above approach also suggests how can be applied to incomplete sweeps or sweeps restricted to a subpopulation. In those cases, the reordered distance matrix should display the characteristic sweep pattern only for a sub-block on the diagonal that is embedded into the diverse remainder of the population.
Intriguingly, we see very little evidence of degradation of the haplotype spectrum by genetic drift. This might have several causes: (i) a very large population size, (ii) a very recent sweep, or (iii) a scenario where sweeps are so frequent that they start overlapping such that exponential amplification of the most fit variants never ceases. The latter is expected in large continuously adapting populations (Neher and Shraiman, 2011). To investigate this matter further, one would need denser time-course data and information on genetic diversity across the genome.
We have investigated the pattern of haplotype variation that arises from mutations occurring during the early phase of a selective sweep on the sweeping haplotype. We found that a selective sweep leaves a characteristic signature in the frequency spectrum of such novel haplotypes: when ordering all adaptive haplotype variants by their abundance, the frequency of the th haplotype variant is approximately times the frequency of the first, where is the mutation rate of the locus and is the selection coefficient of the sweep. The power-law decay of the rank-frequency spectrum is a consequence of exponential amplification (via selection) competing with diversification (through mutation and recombination), similar to the mechanism of preferential attachment originally proposed by Yule (1925). We applied this finding to construct an estimator for the selection coefficient at a candidate locus where a selective sweep is suspected to have occurred recently. Such loci could be either identified by standard approaches to localize sweeps (Fay and Wu, 2000; Sabeti et al., 2002; Kim and Stephan, 2002; Nielsen et al., 2005; Voight et al., 2006; Sabeti et al., 2007), or a priori information could suggest that a sweep has occurred, such as in our HIV examples, where the failure of anti-retroviral therapy implied an adaptive change at the RT locus.
The classic signature used to infer the strength of a recent sweep has been the size of the dip in diversity around the adaptive site (Maynard Smith and Haigh, 1974). Underlying this approach is the rationale that recombination breaks up linkage with increasing genetic distance from the adaptive site: The further away from the sweeping locus, the higher the probability that pre-sweep variation has become unlinked from the adaptive allele and thus remains polymorphic after the sweep. Hence, the width of the dip in genetic diversity scales with the ratio of the selection coefficient and the recombination rate.
While successfully applied in many studies (Kim and Stephan, 2002; Andolfatto, 2007; Macpherson et al., 2007; Sabeti et al., 2002, 2007; Sattath et al., 2011), this approach suffers from a number of shortcomings that can limit its applicability. First, it relies on recombination and thus cannot be applied in organisms which frequently self, such as many plants or yeast, as well as organisms that recombine via horiontal gene transfer, such as bacteria. Recombination rates can also fluctuate strongly along the genome and in time (Winckler et al., 2005; Coop and Przeworski, 2007), rendering their precise estimation difficult, especially in the regions of reduced genetic diversity around a selective sweep. Furthermore, the approach relies on the presence of pre-existing variation at the sweep locus and assumes that we know how much variation was present originally. Ancestral diversity is literally absent in evolution experiments where populations start from a single clone, or in populations where adaptation occurs so frequently that genetic diversity is not fully restored between recurrent selective sweeps (Gillespie, 1994). Finally, for strong sweeps, the dips in diversity can potentially span very large regions extending up to entire chromosomes (Andersen et al., 2012). In such cases, dip-sizes can no longer be assessed accurately.
Our estimator is less dependent on ancestral diversity or accurate recombination rates estimates, since we investigate the new variation that emerges during the sweep and compare selection strength to rate of haplotype diversification (mutation and recombination) rather than to the recombination rate alone. It is essential for our estimator that one can accurately determine haplotype population frequencies on the order of , requiring deep population samples. For example, if we assume a sweep locus with and base our analysis on the five most frequent adaptive haplotype variants (), it would be necessary to measure haplotype frequencies of 2% accurately. This calls for a sample size of roughly . Population genetic data of such depth is already available for HIV, where a large number of sequences can be obtained from plasma samples of infected patients (Tsibris et al., 2009; Fischer et al., 2010; Hedskog et al., 2010). Several efforts are currently being made to achieve a comparable in-depth characterization of the population diversity in several eukaryotes including humans (1000 Genomes Project Consortium et al., 2010), flies (Drosophila Population Genomics Project, 2011), and plants (Cao et al., 2011).
In practice, one typically has a data set (as in the HIV example a sample of size sequences of length 120-200bp) and wants to estimate . The size of the locus corresponds to a certain mutation rate, which will limit the range of selection coefficients that give rise to an observable sweep spectrum: If is too small (say smaller than 10 times the inverse sample size), very little variation will be observed. On the other hand, if is close to one, the founding haplotype diversifies as rapidly as it is amplified and will no longer be the dominant haplotype in the sample. Without a dominant haplotype, the method cannot be applied. Furthermore, our calculations always assume that .
Too large values of can be circumvented by restricting the analysis to only a fraction of the locus, which effectively reduces . Hence, if a sample contains a large number of rare alleles and haplotypes but no dominant haplotype, one can reduce the window size to see whether a dominant haplotype with a trailing sweep spectrum emerges. Similarly, if one analyzes long genomes using a windowing approach and observes a region almost devoid of diversity (corresponding to ), one can increase by increasing the window size until several low frequency haplotypes are observed. One can thus tune the sensitivity of the estimator to different ranges of by adjusting through the length of the locus. A selection coefficient of and a mutation rate of per site, for example, would require a locus of length 10 kb to achieve .
We used the popular program sweepfinder (Nielsen et al., 2005) to compare the performance of our estimator to a traditional approach where the selection coefficient of a sweep is inferred from its signature in the surrounding ancestral diversity. As expected, our estimator consistently outperforms the traditional approach if ancestral diversity is low (Supplementary Fig. 2). Interestingly, even when ancestral diversity is rather high (e.g. , comparable to the level of neutral diversity observed in D. melanogaster), the estimates from our method still have substantially lower variance than those obtained from sweepfinder. Note, however, that the comparison between the two methods is based on rather different data. While the sweepfinder requires long sequences () at moderate coverage, our methods works with much shorter regions () at very deep coverage.
In order for our estimator to be applicable, haplotype variants that were produced after a sweep and rose in frequency by random genetic drift need to be still at low-enough frequency such that they have not yet degraded the sweep signature. The extent of this degradation can be estimated from a simple argument: the lowest population frequency entering our estimator is of the order , while drift-dominated haplotypes will typically be at frequencies , where is the time that has elapsed since the sweep. We thus require , which translates into the condition that population sizes have to be sufficiently large and that sweeps are not too old. As shown in Figure 2, drift results in a growing bulge at low frequencies where the rank spectrum is exponential rather than a power-law. Since the strength of genetic drift is very poorly known, we propose to inspect the ranked haplotype spectrum for deviations from the power-law and choose small enough such that drift dominated haplotypes are excluded from the analysis. In addition, one should check whether the sample is compatible with a near star phylogeny. Failure to exhibit these features should indicate either the absence of a selective sweep, or a sweep that has already been degraded.
While we have focussed on mutation as a source of novel haplotypes, recombination can produce new haplotypes as well. Provided each recombination event results in a unique haplotype, all of the above formulae hold with instead of alone. If, however, recombination rates are not known, recombinant haplotypes can be filtered out by restricting the analysis to only those haplotypes that differ by a single mutation from the founding haplotype. This removes the majority of recombinant haplotypes since recombination with the diverse ancestral population typically incorporates several polymorphisms at once.
Recent studies suggest that in many selective sweeps the adaptive allele has not actually become fixed in the population (incomplete sweep), or that several different haplotypes, all carrying the adaptive allele, have swept simultaneously through the population (soft sweep) (Hermisson and Pennings, 2005; Pritchard et al., 2010; Karasov et al., 2010; Burke et al., 2010). We have demonstrated how incomplete sweeps and soft sweeps can be analyzed by our method (Figure 6B), making use of the fact that the novel haplotype variants our estimator is based on are related to the founding haplotype by a simple star phylogeny, which allows to easily differentiate them from other haplotypes that are not descendants of the founding haplotype.
Our results on the haplotype pattern of a selective sweep were derived under the assumption of a panmictic population of constant size, raising the question how they might be affected by past demographic events such as population expansions or bottlenecks. Most current approaches to investigate selective sweeps rely on the specific patterns a sweep leaves in ancestral neutral diversity. These approaches can thus be very sensitive to past demographic events that shape the patterns of ancestral diversity over a time scale of neutral coalescence, which can include quite ancient demographic events. The approach presented here, however, is fundamentally different. In contrast to ancestral neutral variation, which is typically old, we focus on the very recent variation that arises during the early phase of a selective sweep. Demographic events that occurred prior to the onset of the sweep are thus irrelevant. Only very rapid changes in population size that happen while the adaptive allele is sweeping can cause significant deviations in the haplotype frequency spectra.
The key scenario to be discussed in this context is that of a recent population expansion, since our estimator measures specifically the rate at which new haplotype variants were amplified. Indeed, the haplotype frequency spectrum in an expanding population will resemble that of a selective sweep if the expansion lasted long enough for all coalescence to happen during the expansion and if the expansion rate is large enough such that the spectrum has not yet been eroded by drift, i.e., . In this case, our estimator turns into an estimator of the expansion rate that might be applicable to scenarios such as the expansion of an HIV population in a newly infected individual or the spread of novel strains of influenza. If, on top of an expansion, a beneficial mutation is spreading through the population, the haplotype carrying the beneficial mutation (and its descendants) is expanding faster than the ancestral haplotypes. The estimate of the selection coefficient will then in fact be an estimate of .
The spread of a beneficial mutation in the population generally reduces genetic diversity in the vicinity of the adaptive site. That a selective sweep can also amplify new diversity at very low population frequencies is thereby often overlooked. We have shown that the spectrum of this new variation records the exponential amplification of the novel beneficial allele in a clock-like fashion and can thus be used to estimate its selection coefficient. With the recent advances in sequencing technologies, the required information about low-frequency genetic variation is no longer elusive, making our estimator applicable for a wide range of analyses.
This research was supported by the European Research Council under Grant No. 260686 (RAN). PWM was an HFSP postdoctoral fellow in the Petrov lab at Stanford University. Part of this work was conducted during the program “Microbial and Viral Evolution” at KITP supported by the National Science Foundation under Grant No. PHY05-51164. We thank Dmitri Petrov, Nick Barton, Boris Shraiman, Ben Callahan, Taylor Kessinger, and members of the Petrov lab for helpful feedback.
- 1000 Genomes Project Consortium et al. (2010) 1000 Genomes Project Consortium, R. M. Durbin, G. R. Abecasis, D. L. Altshuler, A. Auton, et al., 2010 A map of human genome variation from population-scale sequencing. Nature 467: 1061–73.
- Andersen et al. (2012) Andersen, E. C., J. P. Gerke, J. A. Shapiro, J. R. Crissman, R. Ghosh, et al., 2012 Chromosome-scale selective sweeps shape caenorhabditis elegans genomic diversity. Nature genetics 44: 285.
- Andolfatto (2007) Andolfatto, P., 2007 Hitchhiking effects of recurrent beneficial amino acid substitutions in the Drosophila melanogaster genome. Genome Res 17: 1755–62.
- Barton (1998) Barton, N., 1998 The effect of hitch-hiking on neutral genealogies. Genet Res 72: 123–133.
- Batorsky et al. (2011) Batorsky, R., M. F. Kearney, S. E. Palmer, F. Maldarelli, I. M. Rouzine, et al., 2011 Estimate of effective recombination rate and average selection coefficient for hiv in chronic infection. Proceedings of the National Academy of Sciences of the United States of America 108: 5661–6.
- Burke et al. (2010) Burke, M. K., J. P. Dunham, P. Shahrestani, K. R. Thornton, M. R. Rose, et al., 2010 Genome-wide analysis of a long-term evolution experiment with Drosophila. Nature 467: 587–90.
- Cao et al. (2011) Cao, J., K. Schneeberger, S. Ossowski, T. Günther, S. Bender, et al., 2011 Whole-genome sequencing of multiple arabidopsis thaliana populations. Nature genetics .
- Coop and Przeworski (2007) Coop, G., and M. Przeworski, 2007 An evolutionary view of human recombination. Nat Rev Genet 8: 23–34.
- Drosophila Population Genomics Project (2011) Drosophila Population Genomics Project, 2011 Drosophila Population Genomics Project. www.dpgp.org.
- Durrett and Schweinsberg (2004) Durrett, R., and J. Schweinsberg, 2004 Approximating selective sweeps. Theoretical Population Biology 66: 129–38.
- Ewing and Hermisson (2010) Ewing, G., and J. Hermisson, 2010 Msms: a coalescent simulation program including recombination, demographic structure and selection at a single locus. Bioinformatics 26: 2064–5.
- Eyre-Walker (2006) Eyre-Walker, A., 2006 The genomic rate of adaptive evolution. Trends Ecol Evol (Amst) 21: 569–75.
- Fay and Wu (2000) Fay, J. C., and C. I. Wu, 2000 Hitchhiking under positive Darwinian selection. Genetics 155: 1405–13.
- Fischer et al. (2010) Fischer, W., V. V. Ganusov, E. E. Giorgi, P. T. Hraber, B. F. Keele, et al., 2010 Transmission of single HIV-1 genomes and dynamics of early immune escape revealed by ultra-deep sequencing. Plos ONE 5: e12303.
- Gillespie (1994) Gillespie, J. H., 1994 The Causes of Molecular Evolution. Oxford University Press, Oxford.
- Haag-Liautard et al. (2007) Haag-Liautard, C., M. Dorris, X. Maside, S. Macaskill, D. L. Halligan, et al., 2007 Direct estimation of per nucleotide and genomic deleterious mutation rates in drosophila. Nature 445: 82–5.
- Hedskog et al. (2010) Hedskog, C., M. Mild, J. Jernberg, E. Sherwood, G. ran Bratt, et al., 2010 Dynamics of HIV-1 quasispecies during antiviral treatment dissected using ultra-deep pyrosequencing. PLoS ONE 5: e11345.
- Hermisson and Pennings (2005) Hermisson, J., and P. Pennings, 2005 Soft sweeps: Molecular population genetics of adaptation from standing genetic variation. Genetics 169: 2335–2352. 10.1534/genetics.104.036947.
- Hernandez et al. (2011) Hernandez, R. D., J. L. Kelley, E. Elyashiv, S. C. Melton, A. Auton, et al., 2011 Classic selective sweeps were rare in recent human evolution. Science 331: 920–924.
- Hudson (2002) Hudson, R. R., 2002 Generating samples under a Wright-Fisher neutral model of genetic variation. Bioinformatics 18: 337–8.
- Hudson et al. (1994) Hudson, R. R., K. Bailey, D. Skarecky, J. Kwiatowski, and F. J. Ayala, 1994 Evidence for positive selection in the superoxide dismutase (sod) region of drosophila melanogaster. Genetics 136: 1329–40.
- Hudson et al. (1987) Hudson, R. R., M. Kreitman, and M. Aguadé, 1987 A test of neutral molecular evolution based on nucleotide data. Genetics 116: 153–9.
- Kaplan et al. (1989) Kaplan, N. L., R. R. Hudson, and C. H. Langley, 1989 The “hitchhiking effect” revisited. Genetics 123: 887–99.
- Karasov et al. (2010) Karasov, T., P. W. Messer, and D. A. Petrov, 2010 Evidence that adaptation in Drosophila is not limited by mutation at single sites. PLoS Genet 6: 1000924.
- Kim and Stephan (2002) Kim, Y., and W. Stephan, 2002 Detecting a local signature of genetic hitchhiking along a recombining chromosome. Genetics 160: 765–77.
- Macpherson et al. (2007) Macpherson, J. M., G. Sella, J. C. Davis, and D. A. Petrov, 2007 Genomewide spatial correspondence between nonsynonymous divergence and neutral polymorphism reveals extensive adaptation in Drosophila. Genetics 177: 2083–99.
- Mansky and Temin (1995) Mansky, L. M., and H. M. Temin, 1995 Lower in vivo mutation rate of human immunodeficiency virus type 1 than that predicted from the fidelity of purified reverse transcriptase. J Virol 69: 5087–94.
- Maynard Smith (1971) Maynard Smith, J., 1971 What use is sex? J Theor Biol 30: 319–335.
- Maynard Smith and Haigh (1974) Maynard Smith, J., and J. Haigh, 1974 The hitch-hiking effect of a favourable gene. Genet Res 23: 23–35.
- Neher and Leitner (2010) Neher, R. A., and T. Leitner, 2010 Recombination rate and selection strength in HIV intra-patient evolution. PLoS Comput Biol 6: e1000660.
- Neher and Shraiman (2011) Neher, R. A., and B. I. Shraiman, 2011 Genetic draft and quasi-neutrality in large facultatively sexual populations. Genetics 188: 975–996.
- Nielsen et al. (2005) Nielsen, R., S. Williamson, Y. Kim, M. J. Hubisz, A. G. Clark, et al., 2005 Genomic scans for selective sweeps using snp data. Genome research 15: 1566–75.
- O’Neil et al. (2002) O’Neil, P. K., G. Sun, H. Yu, Y. Ron, J. P. Dougherty, et al., 2002 Mutational analysis of hiv-1 long terminal repeats to explore the relative contribution of reverse transcriptase and rna polymerase ii to viral mutagenesis. J Biol Chem 277: 38053–61.
- Ossowski et al. (2010) Ossowski, S., K. Schneeberger, J. I. Lucas-Lledó, N. Warthmann, R. M. Clark, et al., 2010 The rate and molecular spectrum of spontaneous mutations in Arabidopsis thaliana. Science 327: 92–4.
- Perelson et al. (1996) Perelson, A. S., A. U. Neumann, M. Markowitz, J. M. Leonard, and D. D. Ho, 1996 HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time. Science 271: 1582–6.
- Pritchard et al. (2010) Pritchard, J. K., J. K. Pickrell, and G. Coop, 2010 The genetics of human adaptation: Hard sweeps, soft sweeps, and polygenic adaptation. Current Biology 20: R208–R215.
- Przeworski (2003) Przeworski, M., 2003 Estimating the time since the fixation of a beneficial allele. Genetics 164: 1667–76.
- Roach et al. (2010) Roach, J., G. Glusman, A. Smit, C. Huff, R. Hubley, et al., 2010 Analysis of genetic inheritance in a family quartet by whole-genome sequencing. Science 328: 636.
- Sabeti et al. (2007) Sabeti, P. C., B. Fry, J. Lohmueller, E. Hostetter, C. Cotsapas, et al., 2007 Genome-wide detection and characterization of positive selection in human populations. Nature 449: 913–8.
- Sabeti et al. (2002) Sabeti, P. C., D. E. Reich, J. M. Higgins, H. Z. P. Levine, D. J. Richter, et al., 2002 Detecting recent positive selection in the human genome from haplotype structure. Nature 419: 832–7.
- Sattath et al. (2011) Sattath, S., E. Elyashiv, O. Kolodny, Y. Rinott, and G. Sella, 2011 Pervasive adaptive protein evolution apparent in diversity patterns around amino acid substitutions in drosophila simulans. PLoS Genet 7: e1001302.
- Slatkin and Hudson (1991) Slatkin, M., and R. R. Hudson, 1991 Pairwise comparisons of mitochondrial dna sequences in stable and exponentially growing populations. Genetics 129: 555–62.
- Tajima (1989) Tajima, F., 1989 Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585–95.
- Tsibris et al. (2009) Tsibris, A. M. N., B. Korber, R. Arnaout, C. Russ, C.-C. Lo, et al., 2009 Quantitative deep sequencing reveals dynamic HIV-1 escape and large population shifts during ccr5 antagonist therapy in vivo. PLoS ONE 4: e5683.
- Voight et al. (2006) Voight, B. F., S. Kudaravalli, X. Wen, and J. K. Pritchard, 2006 A map of recent positive selection in the human genome. PLoS Biol 4: e72.
- Wiehe and Stephan (1993) Wiehe, T. H., and W. Stephan, 1993 Analysis of a genetic hitchhiking model, and its application to DNA polymorphism data from Drosophila melanogaster. Mol Biol Evol 10: 842–54.
- Winckler et al. (2005) Winckler, W., S. R. Myers, D. J. Richter, R. C. Onofrio, G. J. McDonald, et al., 2005 Comparison of fine-scale recombination rates in humans and chimpanzees. Science 308: 107–11.
- Yule (1925) Yule, G. U., 1925 A mathematical theory of evolution, based on the conclusions of Dr. J. C. Willis, F.R.S. Philos Trans R Soc Lond, B, Biol Sci 213: 21–87.
Appendix A Neutral haplotype spectrum
The haplotype spectrum expected in a haploid neutral Fisher-Wright model without recombination can be calculated from the Ewens sampling formula (Ewens:1972p35481). Ewens showed that the probability of a sample of size is
where is the number of allele classes that are sampled times and with . The expectation of is therefore given by
where the last two approximate inequalities are accurate if and , respectively. Hence the expected number of allele classes with more than members is roughly , where cutting off the sum at approximately accounts for the exponential. With this approximation, the th abundant allele class is expected to contain
copies of the allele. A more accurate expression of the spectrum is obtained by determining the such that , using the exact expression given above. This numerical solution for the haplotype spectrum is plotted in Figure 2B of the main text.
Appendix B The distribution of haplotype frequencies
In the main text, we calculated the distribution of the establishment time of the th haplotype and the frequency of the corresponding haplotype. Here, we show how the joint distribution of all seeding times and the resulting frequency spectrum can be calculated assuming that the novel haplotypes are rare and evolve independently, which is justified if they constitute a small share of the total population, i.e., if . In this case, the probability that haplotypes are present in frequencies is given by
where is the probability that a haplotype has frequency at time given it became established at time . The distribution of establishment times is given by
where is the rate of establishing novel adaptive haplotypes (main text Equation (1) and below). Note that the defined in Equation (15) are not ordered. They are distributed according to a Poisson point process with density . Assuming that established novel haplotypes increase in frequency logistically according to Equation (5) of the main text, we have
where is the Dirac -function (the stochastic analog is calculated below, see also (Desai:2007p954)). Substituting into Equation (14) and integrating over , we obtain
with . Haplotypes that are common after the sweep are most likely seeded early during the sweep. Furthermore, we showed in Equation (5) of the main text that their relative frequencies stay approximately constant during the amplification phase. Hence we can determine the joint distribution of frequencies at early times while is still exponential. After substituting the and simplifying, we find
where we dropped factors independent of which ensure normalization. A very similar result was found in (Desai:2007p954). At large , the form of the prefactor changes due to the saturation of the allele frequency at , but the distribution of the frequencies of the haplotypes that were seeded early during the sweep remains of this form until the spectrum is eroded by genetic drift.
The haplotype spectrum therefore decays with a power , which is consistent with the power obtained for the cumulative or rank spectrum (integrating yields ). More importantly, this result tells us that the distribution of haplotype frequencies conditional on the number of haplotypes observed is approximately independent of if . Hence, given that a sweep occurred, all information about the strength of the sweep is contained in the number of haplotypes and the precise values of their frequencies do not contain any additional information if . However, whenever there are deviations from the assumptions made here, the haplotype frequencies will contain additional information.
Appendix C Stochastic derivation of the haplotype spectrum
The dynamics of rare haplotypes are strongly influenced by random genetic drift and we have to ascertain the deterministic arguments made in the main text by a more careful stochastic calculation. While hard in general, an approximate analytic calculation of the frequency spectrum of rare haplotypes is feasible in our case for the following reasons: (i) The dynamics of a beneficial allele are essentially deterministic since it is much more frequent than haplotypes that arise through secondary mutations. (ii) The dynamics of rare haplotypes can be described by a linear branching process since they are always a small fraction of the total population.
As already done in Equation (14), we decompose the distribution of haplotype frequencies into the distribution of times when the novel haplotypes arise and probability that a haplotype is present in copies at time , given it arose at time . We will derive first and consider the spectrum due to the superposition of several independent seeding events below.
c.1 Distribution of rare variants arising in a logistic sweep
To model the stochastic dynamics of rare haplotypes, we use a continuous time branching process in which individuals produce identical copies of themselves with rate and die with rate , i.e., the unit of time is chosen to be the generation time. The average number of offspring of a given individual in this model is . Hence, is the growth rate of the haplotype carrying the beneficial allele. In the case of a sweep, we have
where the first term accounts for selection ( is the frequency of the beneficial allele) and the second term accounts for mutations that change the state of the haplotype. The dynamics of are described by the forward Master equation
which accounts for replication (first term) and death (second term). To solve for , it is useful to consider the generating function , which obeys the equation
with initial condition . This equation can be solved via the method of characteristics, with the result
where we have used along the way. The latter is a good approximation if selection is weak in one generation and amounts to neglecting terms of order . We will now substitute the explicit expression for , where it will be convenient to parametrize the frequency of the beneficial allele as with . Using this form of , we find for the generating function
where . Any haplotype that is abundant enough to be sampled with high probability most likely originated in the early phase of the sweep (), which allows for the approximation where is the sample size. Furthermore, we will typically observe the spectrum at times when the sweep is almost complete. Hence we can approximate where is the frequency of the deleterious wild type allele at the time of sampling. Using these simplifications, we obtain
This expression is straightforwardly expanded into a geometric series in whose coefficients are . For large , one finds
with relative corrections being on the order of . The quantity is the mean copy number of the haplotype conditional on non-extinction, and the two terms contributing to have a straightforward interpretation: The first term is the contribution of selection, which amplifies the haplotype before the fixation of the beneficial allele. The second term is the contribution of random genetic drift, which evaluates simply to in the limit of . The latter is the analog of the well known fact that a non-extinct neutral allele in a neutral Moran process is on average present in copies after generations. The expression for exhibits a crossover from an early regime where selection dominates to a random drift dominated regime at large . In the limit , we have
This crossover will inform us below about how long the contribution of random drift can be neglected when applying our estimator of the strength of the selective sweep.
c.2 The haplotype frequency spectrum
Having calculated the copy number distribution of a haplotype that originated at time , we now have to determine the distribution of seeding times and calculate the resulting spectrum of haplotype frequencies. New haplotypes that contain the beneficial allele are produced at rate
where, as before, is the frequency of the sweeping allele. Note that this differs from the rate of establishment of novel variants by a factor , which will reemerge from the stochastic calculation. The deterministic approximation for is valid if it is unlikely that new variants are seeded before establishment of the founding variant, which requires (see Desai:2007p954). Since novel haplotypes are seeded and evolve independently to a good approximation, the number of haplotypes present in copies at time is Poisson distributed with mean
Due to the exponential nature of , it is convenient to sum over and calculate the expected number of haplotypes with copy numbers greater than . The sum is well approximated by the integral , and we have
where the last approximation assumes that novel haplotypes are seeded while the beneficial allele is still expanding exponentially, i.e., .
It does not seem possible to evaluate the integral over in Equation (29) analytically. However, the integral is dominated by contributions from a well defined time interval and can be evaluated perturbatively. For a very large population where drift is negligible, and for , this integral simplifies to
which is very sharply cutoff for . Hence we can send the upper integration boundary to infinity without loss of accuracy and evaluate this integral exactly. One finds . The contribution to this integral come from a narrow peak of width and height . Genetic drift and mutation predominantly change only this height and width, leaving the shape of the integral approximately invariant. Hence we can evaluate this integral by calculating where the integral peaks and how wide this peak is (Laplace’s method).
Including the correction due to drift and finite term corrections, the integrand peaks when , which translates into as opposed to without the corrections. The second derivative of the logarithm of the integrand is approximately given approximately by
where is the age of the haplotypes. Hence the peak dominating the integral becomes wider by a factor . With these corrections to the “height” and the “width” of the integrand, we obtain the approximate expression
for the integral. This expression is accurate as long as and . The age of the haplotype evaluates approximately to , which is of order . The additional factor accounts for the additional time the older haplotypes have been degraded by mutations, while the accounts for the contribution of drift to the frequency of the haplotypes.
After the sweep, the frequency of the most common haplotype is and the expected number of haplotypes above frequency is given by
This result tells us that the mean number of haplotypes with frequencies greater than is approximately linear in and decreases with approximately as . Furthermore, the expected number of haplotypes above is increasing with time, since rare haplotypes increase in frequency due to random drift.
Given that we observe haplotypes at frequencies higher then , we can use Equation (33) to estimate :
This equation differs from Equation (9) of the main text by a reduction of due to random drift, which has been ignored in the simple deterministic derivation given in the main text. This reduction allows us to correct for the effects of genetic drift as long as drift is not to strong. Obviously, the correction fails as approaches .
Equation (34) informs us about the regimes where the proposed methods to estimate the selection coefficient is likely to work. Random genetic drift will degrade the signature of the sweep for haplotype frequencies smaller than . Since the time needed for completion of the sweep is on the order of , we require . Since is itself small, we need for the method to work. The breakdown of the method is clearly seen in Figure 3 of the main text once falls below 100.
Appendix D Pruning recombinant haplotypes
Haplotypes that arise by mutation from the founding haplotype differ at exactly one position from the founding haplotype, while haplotypes that result from a recombination event with a member of the diverse ancestral population typically differ at several positions. Furthermore, haplotypes that are mutants of mutants will differ at two positions from the founding haplotype. In the main text, we argued that one can restrict the haplotype spectrum to those haplotypes that differ only at a single site from the founding haplotype, thereby removing most recombinant haplotypes and mutants of mutants. Here, we show that frequency spectrum of such a restricted set of haplotypes differs slightly from that of all haplotypes.
As before, the frequency of the beneficial allele will typically follow Equation (1) of the main text. The frequency of the founding haplotype, however, will remain below this frequency due to loss through recombination and mutation.
where we have abbreviated the initial growth rate of a haplotype by . Mutations on this haplotype establish with rate . Haplotypes that establish at time then typically follow a frequency trajectory
The most likely seeding time of the ith haplotype is given by
Hence we obtain for the ratios of haplotype frequencies at times when the beneficial allele is near fixation
This differs from Equation (7) of the main text in that the ratio is proportional to , rather than . This difference is due to the fact that here haplotypes grow with the same rate as the rate at which they are seeded. In the previous case where all haplotypes are considered, haplotypes grow with rate , while the seeding rate is proportional to the frequency of the beneficial allele which grows with rate .