Allele age has long been a focus of population genetic research, primarily because it can be an important clue to the fitness effects of an allele. By virtue of their effects on fitness, alleles under directional selection are expected to be younger than neutral alleles of the same frequency. We developed a new coalescent-based estimator of a close proxy for allele age, the time when a copy of an allele first shares common ancestry with other chromosomes in a sample not carrying that allele. The estimator performs well, including for the very rarest of alleles that occur just once in a sample, with a bias that is typically negative. The estimator is mostly insensitive to population demography and to factors that can arise in population genomic pipelines, including the statistical phasing of chromosomes. Applications to 1000 Genomes Data and UK10K genome data confirm predictions that singleton alleles that alter proteins are significantly younger than those that do not, with a greater difference in the larger UK10K dataset, as expected. The 1000 Genomes populations varied markedly in their distributions for singleton allele ages, suggesting that these distributions can be used to inform models of demographic history, including recent events that are only revealed by their impacts on the ages of very rare alleles.
The age of an allele of a given frequency can be reveal the forces acting upon it, with rare alleles being particularly sensitive to recent evolutionary processes [1]. A functional allele that is younger than expected given its frequency is likely to have been under directional selection. This is not surprising for favored alleles, but it is also true for harmful alleles [2], including those with negative impacts on health that are under negative selection [3]. If researchers are able to estimate allele age, they could combine this with other information (e. g. allele frequency, geographical distribution, functional annotation) to improve predictions of an allele’s effect on human health. Alternatively, an allele that is older than expected given its frequency is also a candidate for having an interesting history, as functional alleles older than expected can be the result of balancing or negative frequency-dependent selection [4]. As population genomic samples grow in size, the density of variable sites rises approximately in proportion to the log of the sample size [5], and very large data sets will have large numbers of SNPs and other variants in every gene. If information could be gleaned on the ages of a large number of variants for a functional region of the genome, this could be used to develop a detailed portrait of the history of natural selection specifically on that region. Allele ages are also shaped by processes that act in aggregate across the genome. The overall distribution of ages will be strongly shaped by the demographic history of the population, and for the rarest alleles, that distribution will be acutely sensitive to recent admixture [6]. The age spectrum will also fluctuate spatially along the genome, both stochastically and as a function of the intensity of background selection [7]. In developing a way to study allele age, we considered several constraints. An estimator should not be a function of allele frequency, as we wish to glean information about allele history that is distinct from its frequency. We also prefer an approach that is not a function of the demographic history of the population, as some estimators are [8,9]. An estimator that can get close to the true value of the unknown, regardless of demographic history, enables analyses in cases when the history is not known and it enables comparisons between populations that are not confounded by errors in our knowledge of the demographic history of the populations. We also wish to be able to study the age of the very rarest alleles, including those that appear only once in a sample (singletons). This last criterion leads to an approach that is different from existing methods that are based on the variation observed among copies of an allele, or in flanking regions [8–13]. An estimator should also be applicable for very large sample sizes for which it becomes increasingly possible to find low frequency alleles that arose by multiple mutations [14]. For these cases, as with singletons, we need a method that is applicable for each individual gene copy. Finally, an estimator should not be highly sensitive to the details of the bioinformatics pipeline used to process the data, such as whether the data were statistically phased. We developed a new estimator that focuses, not directly on allele age, but rather on the time when a base position in a particular chromosome first coalesces in the genealogy. The mutation causing a derived allele at that base will have occurred since this first coalescent time, and so first coalescent time can be used as a proxy for allele age. For example, a neutral singleton allele will have a uniform probability of having arisen anywhere along an external branch, and therefore an expected age of half the first coalescent time. We assessed performance using simulated data, and show that it performs well and substantially overcomes the challenges describe above. We also applied it to SNP alleles from the 1000 Genomes Project [15] and the UK10K genome panel [16]. For these data sets, we assessed basic predictions regarding population-specific variation in the ages of rare alleles, the ages of private alleles and alleles shared by populations, and we compared ages of alleles that are expected a priori to have functional impacts with those that are not. Our estimator shares structural similarities with existing methods developed to estimate demographic histories from shared-haplotype tract-length distributions [17–20], but is uniquely able to discern the specific histories of individual rare alleles.
Consider one chromosome (the focal chromosome) from a sample of chromosomes drawn at random from a population, and an individual base position on that chromosome (the focal base). As shown in Fig 1, the focal base can be thought of as the terminal point on a branch A of the genealogy of the sample of chromosomes at that base position. If the focal base is a singleton variant (i. e. a derived allele that occurs only once in the sample), then the mutation causing that allele must have occurred on this branch. Also shown in Fig 1 is that branch A connects the focal chromosome to the most recent common ancestor of that chromosome and a sister branch S, which is ancestral to one or more sister chromosomes of the focal chromosome. Branches A and S share their most recent common ancestor at a time tc generations prior to the sample generation, and branch A is therefore tc generations long. We denote the length of branch S as φ generations. If the sample contains a single unique sister chromosome (that is, the focal chromosome is also the closest relative of the sister at the focal base), then φ is equal to tc. When multiple chromosomes are all equally and most closely related to the focal chromosome at the focal base, they form a clade of sisters, and φ < tc. Consider first the case when there is no recombination and the focal chromosome has just a single sister. In this situation, any differences between the focal and sister chromosomes will have been caused by mutations on branches A and S. We model mutation as a Poisson process, where each base has a constant probability, μ, of mutating in each generation. Treating a chromosome as a continuous line, and considering events just to one side of the focal base (in either the 5’ or 3’ direction), the probability density of distance x from a focal base to the first base that is not identical between the two chromosomes can be approximated with an exponential distribution having a rate of μ2tc: p (x) =μ2tce-xμ2tc. (1) If we knew which chromosome was the sister chromosome, we could compare them and identify the distance from the focal base to the nearest difference between the chromosomes (i. e. x), and use this to estimate tc. This basic idea captures the underlying rationale of our approach. The final formulation takes into account the remaining issues: that we do not know which chromosome is the sister to the focal chromosome; that the focal chromosome may have multiple sisters; and that Eq (1) assumes no recombination (see Materials and methods). Ultimately, an expression that resembles Eq (1), but that differs in replacing x with the longest observed tract of identity between the focal chromosome and each of the other chromosomes in the sample, proves applicable. We call this quantity the maximum shared haplotype (msh), and show that it arises on either branch A or S (Fig 1) with high probability. We use t^c to denote the estimator, and because tc values range over several orders of magnitude and Poisson processes have variances proportional to their means, we focused primarily on the logarithm, log10 (t^c). Performance was assessed in terms of the root mean squared error (RMSE), bias (mean of estimated minus true values), and correlation (Pearson’s r for the true and estimated values of log10 (tc) ) for alleles at all frequencies in a series of large simulated data sets. We varied sample size and recombination rate, and considered three demographic histories that varied in terms of population sizes, exponential growth, historical bottlenecks, and intercontinental migration. We also considered samples of chromosomes with known and with estimated phase. All of these results are summarized in S1 and S2 Tables. For low and intermediate recombination, over a wide range of circumstances, the estimator exhibits an RMSE of about 0. 4 log-transformed generations, corresponding to estimates that are typically within a factor of 2. 5 of the true value. Across the models, correlations of true and estimated values ranged from 0. 4 to 0. 95 with a mode of 0. 9. For recombination rates equal to or less than the mutation rate, bias varies from -0. 4 to 0. 1 with a mode of -0. 2, which corresponds to an average underestimate by a factor of 0. 63. Performance was consistent across the spectrum of allele frequencies, and with respect to particular independent variables, performance was better: when recombination was low; when population size was constant; when phase was known; and when sample sizes were larger. A key factor that determines how informative estimates are is having chromosomes that are much more closely related to their closest relatives than to unrelated chromosomes. Thus, strong recent growth can produce more star-like genealogies that reduces these differences and reduce the quality of the estimates, while larger sample sizes improve them. When recombination rate is appreciably higher than the mutation rate, such that typical msh values are small enough to be within the range of distances between average pairwise SNPs, estimates worsens and the bias shifts from negative to positive (S2 Table). The amount of recombination that is too high relative to the mutation rate will depend on the length of the region of high recombination (e. g. if it is associated with hotspots of short length), and on sample size and the demographic history of the sample. We observed that larger sample sizes exert a greater improvement on the quality of estimates in the case of high recombination than with low to intermediate recombination (S1 and S2 Tables). The change in sign of the bias, with high recombination, and that the observed absolute value of the bias is lower with intermediate levels of recombination, suggests that there are multiple contributions to the bias. Fig 2 shows plots of estimated versus true values for different ranges of allele frequencies and for two different sample sizes (box plots are shown in S1 Fig). Also shown in Fig 2 is an idealized estimator of log10 (tc) that is based on allele frequency and for which all alleles at a given frequency generate the same estimate–shown as a band of red points. We show both the msh-based and the idealized frequency-based estimators to highlight the contrast; the former does not make use of information about a variant’s frequency while the latter does not make use of any information about a variant’s msh value. Unlike the msh-based estimator, where each variant may have a unique msh from which to generate a unique tc estimate, no frequency-based estimator can distinguish among the potentially large number of variants occurring at identical frequencies within a sample. When considered over the full range of allele frequencies, the idealized estimator can explain nearly as much of the variation in log10 (tc) as the msh-based estimator (Fig 2f). However, for rare alleles, the msh-based stimator retains strong performance, whereas the frequency-based estimator explains little to none of the variation in log10 (tc). An important application is assessing upon which of the two chromosomes of an individual a singleton allele correctly resides. As shown in Fig 3, the phasing accuracy for singleton variants in simulated data rises as the ratio of the genetic lengths of the alternative msh tracts diverge. For variants with msh tracts of similar length (and thus high probability of misassigned phase), similar t^c values will be found regardless of how phase is assigned. For non-singleton variants, the estimator is expected to be relatively immune to switch errors introduced by haplotype phasing software. This is because switch errors typically involve low frequency variants among similar haplotypes [21,22], and these do not typically affect the distribution of msh values which are often terminated by relatively common variants or recombination events between unrelated haplotypes. S2 Fig shows that results for statistically phased chromosomes are quite similar to those for the correct chromosomes (r = 0. 941) based on analyses of singleton alleles in male X-chromosome UK10K data. From S1 and S2 Tables we see that mean error is slightly increased when the data (including singletons) are phased statistically. The bias also becomes more negative by a small amount with phasing. Singletons are phased by assigning them to the chromosome that reveals the shorter msh and thus the longer t^c. The proximal effect of this will be to introduce a positive bias that applies in those instances when this phase assignment is not correct. However, the observation that bias becomes slightly more negative with phasing, including for singletons, suggests a greater effect in the other direction. It is possible that the phasing of singletons prior to determining msh values causes them to cluster on fewer chromosomes, thereby lengthening the msh values that are observed. Based on analysis of variation (ANOVA, Table 1, S6 Table), we learned that different populations show different distributions of log10 (t^c) values for singleton variants (d. f. = 25; F = 19248; p < 1 × 10−128). As shown in Fig 4 and S3 Table, populations in Africa have higher geometric means (6353 to 7178 generations) than populations from regions that have not had substantial admixture from African populations (East Asia, South Asia, and Europe), which have lower geometric means (2018 to 3882 generations). Admixed American populations spanned a wide range (5610 to 7194), with values near those of the African populations. The finding of younger ages for singletons from non-African old-world populations is expected under a general Out-of-Africa model in which those populations have passed through a bottleneck and have had an overall lower effective population size and thus shorter coalescent times. The 1KGP data allow us to compare the age of rare variants that are found only in a single population to the age of variants that have the same low frequency within that population, but that also are found in one or more other populations. For all of the 1KGP populations, private singletons had distributions that were shifted to the left (younger) relative to singletons that were shared. Fig 5 shows these distributions for a single population that is representative of those observed from each the five super-populations, each of which showed a characteristic distribution (S9–S13 Figs). In every population, private singleton alleles are younger than singleton alleles that are also shared with other populations (d. f. = 1; F = 896298; p = 0, Table 1), with private singleton variants being 70% younger than comparable shared variants. However, the size of the reduction varies considerably by population (d. f. = 25; F = 3437; p = 0, Table 1), with the largest reduction observed in admixed American populations (82% to 92%). With respect to the functional impact of rare alleles, we expect that the rarest variants in a sample will be enriched (relative to more common variants) for alleles that have an impact on fitness, and therefore that rare alleles will be targeted more by natural selection. This is because such alleles pass through frequency space more quickly than neutral alleles and, if not lost from the population, will reach a given frequency more quickly than a neutral allele [2]. Consistent with this prediction we observed in the 1KGP data that, conditional on population and geographic spread, singleton variants that alter proteins are 3. 6% younger than those that do not (d. f. = 1; F = 24; p = 1 × 10−6, Table 1). There is no significant difference to the effect of protein-changing status in individual populations (d. f. = 25; F = 0. 58; p = 0. 95, Table 1), and there is no significant interaction between protein-changing status and a variant being private to its population (d. f. = 1; F = 2. 86; p = < 9 × 10−2, Table 1). The theory that says that alleles that affect fitness will be younger applies to both harmful and beneficial alleles [2]. However, if it is the case that harmful rare alleles are much more common than beneficial rare alleles, then the age difference between functional alleles and neutral alleles should be greater the rarer are the alleles under comparison. This is because most harmful alleles are rapidly removed from a population and are more likely to be observed at the lowest frequencies in larger samples. We therefore predict that we should observe a greater difference in ages between alleles that change proteins and those that do not in the larger UK10K dataset than in the 1KGP data. In fact, the singletons in the 7,242 samples of UK10k data showed a considerably larger effect of protein-changing status than those identified as singletons in the 100-sample populations of the 1KGP data. The geometric mean t^c value for the 105 protein-changing singleton variants was 12% lower than it was for the 1. 5 × 107 non-protein-changing variants (305 generations vs 347 generations; Wilcoxon-Mann-Whitney p = 0).
Every mutation can be envisioned as occurring on a branch of the genealogy or gene tree for a sample of genomes at the locus where the mutation occurred. Previous estimators of allele age have focused on the time point at which different copies of an allele coalesce with each other, i. e. the time of the top of the gene tree edge upon which the mutation occurred [8–12]. Here we have taken a different approach and focused on the time at the bottom of that same edge, the time of first coalescence for an edge carrying a mutation. These first coalescent times have direct connections to msh values, which are easily measured from aligned genomes. The resulting estimator is not without noise, but estimates covary roughly linearly with true values, with moderate error and bias over a wide range of circumstances. Our estimator also meets a variety of desired criteria established at the outset. It is not a function of allele frequency, and thus can be used to study how an allele came to reach its observed frequency. It is only very weakly a function of demographic history, and thus it can be used to compare the ages of alleles of different populations that have unknown or widely varying histories. It can be applied to alleles that occur only once in a sample, and to genomic data with very large sample sizes. It can be applied individually to each copy of an allele, and thus can be used in cases when an allele has arisen by multiple mutations. And it suffers little degradation in performance when run on statistically phased chromosomes. The estimator is not expected to be highly sensitive to some additional issues that can arise with population genomic data, including sequencing errors. The UK10K data for example are low coverage (7x), however, because new singleton variants are only called when the data are strongly supportive (high false negative rate, low false positive rate) only a small minority of the singletons in the data set are expected to be errors. In fact, the proportion of called variants in the case of the UK10K sample that are true variants was estimated to be quite high (~94% for singleton variants), based on measurements across monozygotic twins. This is shown in the first figure of the extended data for that paper (panel k, row AC = 1, MZ twins section, column" non-ref genotypes %" , divided by two and subtracted from 100) [16]. However the variants that are missing from low coverage data can affect the distribution of msh values in ways not accounted for by the method. For example, many of the mutations that terminate tracts of identify and determine msh are themselves singletons, and so a coverage bias against singletons will tend to lengthen msh values and reduce t^c values. Although the method can be used to study selection, implicit in the method is an assumption of neutrality, such that mutations on branches A and S (Fig 1) do not affect the distribution of tc values. But of course a significant fraction of the mutations in evolutionarily constrained regions are not neutral, and the tc values in these regions are reduced, as we see in our ANOVA results. The question then is, how does this kind of variation in the non-deleterious mutation rate affect msh and t^c values? A closely related question applies to factors such as variation in background selection and flucations in rates of gene flow or admixture, that cause polymorphism levels to wax and wane across chromosomes, and that alter the distribution of msh values that terminate due to recombination events. The basic expectation is that these kinds of variations will have a greater effect on long tc values, when flucations in the detminants of what should be short msh values can have a greater affect. However short tc values are associated with very long msh values that span large portions of chromosomes and can be expected to be more immune to local flucations in the factors that terminate msh tracts. Another potential difficulty is that as sample sizes grow very large, some rare alleles will have been caused by multiple mutations [14]. While this is not a concern for singletons, it is possible that, for example, a doubleton will actually represent two mutations. In these cases, the likelihood estimator will simply generate a composite of two different singleton coalescent times. Unlike methods that estimate the time of most recent common ancestry among individuals carrying a rare variant, t^c is the estimated time since each copy of the rare variant last shared a common ancestor with an individual not carrying the variant. If two copies of an allele are actually the result of independent mutations, the composite t^c estimation will produce an intermediate value and not an aberrantly large one, as would occur with a method based on intra-allelic differences. S3 Fig illustrates the impact of analyzing two independently derived minor variants as copies of a single mutation. Some additional and important benefits of the estimator described here emerged during the course of this study. Estimated tc values, or msh values, can be used to phase singleton alleles (i. e. estimate the correct chromosome for placement in diploid heterozygous individuals). Another benefit is that the method can be applied relatively quickly to very large data sets with the benefit of the PBWT algorithm [23]. Computing time increase linearly with both sample size and numbers of variants, and all the singleton variants for the full UK10K panel for chromosome 22 can be analyzed in under 30 minutes, for example (S4 Table). A third likely benefit, that has not been explored here, is that it should be possible to work with genomic data in which only dispersed portions of the genome have been sequenced. In particular, msh tracts for rare alleles from large samples are typically much longer than the distances between genes, and consequently it should be possible to apply the method to data sets of aligned exome sequences. Our analysis of variants of 1% frequency among 26 populations of the 1KGP data, revealed considerable systematic variation in allele ages depending on the variant’s function, location, and geographic distribution. Different human populations have variants at 1% frequency with very different distributions of ages. Population-specific demographic histories, including population bottlenecks, expansion, and admixture events, have likely all contributed to these broad differences, and therefore changing, in effect, what it means to be a 1% frequency variant on a population-by-population basis. We also observed that the ages of variants found at 1% frequency in one population depend greatly on whether or not they are also found in other populations. To be found in multiple populations, a variant that arose by a single mutation must have either: 1) originated prior to the divergence of those populations and risen to sufficient frequency that it has not been lost in either; or 2) persisted in its population of origin long enough and risen to sufficient frequency that migrants have had an opportunity to bring it to at least 1% frequency in another population. Both of these phenomena are reflected in substantial differences of ages of variants found in a population at 1% frequency depending on their total geographic range. The survey of singleton t^c distributions for 1KGP revealed bimodal distributions in each of the non-African populations (S9–S13 Figs). All of them exhibit a large peak at less than a thousand generations, and another smaller one at more than ten thousand generations, a pattern that is readily interpreted in terms of an older bottleneck associated with the out-of-Africa history of these populations. S14 Fig shows these same patterns arise for simulated data generated with parameter estimates for an Out-of-Africa model [24]. The younger larger peak is characteristic of variants that arose after the population’s ancestors had emigrated from Africa. These are relatively new variants that have not yet had an opportunity to rise to high frequency or spread to other populations. The older peak consists of variants that arose within an ancestral African population and predate the modern human expansion into Eurasia. Although they are identified as singletons in individual populations, these variants are often shared across multiple populations, and they typically occur at frequencies greater than 1% in some populations. In the populations where they are ascertained as singletons, founder events, genetic drift, and natural selection have made these older alleles rare, or even eliminated them entirely before re-introduction by migration from other populations. Independently of geography and demographics, protein-changing variants in the 1KGP populations that are found at 1% frequency are younger on average than comparable variants that do not change proteins, consistent with previous reports [25]. These are expected to be mostly deleterious variants that have not yet been removed from the population, but they may also include some beneficial variants that have not yet been pushed to higher frequencies. At the 1% frequency level, the impact of protein-changing status on the distribution of allele ages is statistically significant, but considerably smaller than the systematic differences introduced by demography and geography. Among the singleton protein-changing variants ascertained in the UK10k data, however, we find a larger difference with respect to variants that do not change proteins. The contrast with the 1KGP data is consistent with the younger, rarer, and more numerous variants of the UK10K data are subject to greater selective forces than the variants at 1% frequency in the 1KGP populations. Our ANOVA did not reveal evidence that local adaptation contributes to the variance in ages of alleles at 1% frequency. If selectively important alleles were disproportionately prevented from migrating between populations or remaining in multiple populations, we would have observed the pool of shared protein-changing variants to be enriched for protein-changing variants without selective function. There would have been less differentiation between non-shared and shared, non-protein-changing variants, and there would have been a significant (negative) interaction term between private status and protein-changing status, which we did not observe. We also did not detect an interaction between populations and protein-changing status. If different populations were broadly experiencing greater or lesser amounts of directional selection on rare variants, or if there were substantial differences between populations in the ability of natural selection to remove harmful variants or raise the frequency of beneficial variants (such as due to variation in effective population size), we would expect to have seen a significant interaction between populations and protein-changing status. However, we did not; and while both of these phenomena may be taking place, they are not major drivers of the distribution of ages of variants found at 1% frequency, in contrast to the contributions of that demography, geographic spread, and global properties of natural selection. These population genetic results are all presented as conditioning on alleles at 1% frequency that have been down-sampled to equal sizes. By considering only alleles of a particular frequency, we are able to draw meaningful comparisons and contrasts among populations. While the formula for the estimator itself is not a function of the frequency of the allele to which it is being applied, the distributions of t^c estimates will vary greatly as a function of the frequencies of the alleles being studied (see e. g. S15 Fig). Indeed, it is because of these very different age distributions that alleles of different frequencies can have very different histories, including the action of natural selection and the amount of movement among populations. Any future work that would make comparisons between populations or estimate demographic models based on t^c estimates will need to explicitly condition on allele frequencies. With an estimate of an allele’s age, an investigator has an important new piece of information to bring to bear on the possible functional impact of a mutation. As shown here, for 1KGP data and even more strongly for the larger UK10K data, functional alleles are younger, and therefore, concomitantly, any allele that is discovered to be especially young, given its frequency in the sample, is a good candidate for having an effect of fitness. In this context, an increase in sample sizes will have multiple important effects. First, the rarest alleles in large samples will be rarer and younger, on average, than those found in smaller samples, and they will thereby be relatively enriched for more alleles of harmful effect and for alleles of more harmful effect (these are the alleles that would not reach those higher frequencies observable with smaller sample sizes). Second, the numbers of alleles in the rarest class rises with sample size, and thus so does the number of very rare alleles observed for a given gene. For example, the average number of autosomal singletons for the 26 1KGP populations (100 genomes) was 3,311,984, while the count for the UK10K data (7242 genomes) was approximately 6 times higher at 19,078,777. Based on those values, and assuming the number of SNPs is a function of the log of the sample size [5], a sample of 2N = 100,000 genomes would have 36 million singletons, and a sample of 1 million genomes would have 45 million singletons. The increasing density of very rare alleles, as sample sizes grow, opens the door to a kind of mapping of functional constraint across a gene that will become increasingly fine-grained. Third, as data sets get very large so grows the potential for studying the impact of selection on variation that has arisen at different times and it will become increasingly possible to assess whether the action of selection has been changing for different genes or regions of the genome. The estimator can also be used to study forces that shape the overall distribution of allele ages across the genome, particularly demographic forces. As shown in Fig 5, the distributions of ages of singletons across 1KGP populations varies greatly, as do the differences between the ages of private and shared alleles. These distributions can be used in principle to develop models of demographic history, including in particular recent events that are only revealed by their impacts on the ages of very rare alleles [26].
We applied the estimator of tc to low frequency alleles in data from 26 human populations, grouped into five super-populations, from the 1000 Genomes Project (1KGP) data set [15]. We sub-sampled each population to achieve a uniform sample size of N = 50 diploid individuals. Chromosomes were taken as phased by the 1000 Genomes Consortium using SHAPEIT [30] with the exception that we re-phased singleton variants in each population based on relative t^c values. To do this, for each heterozygous singleton, we masked all other singleton variants, calculated t^c for both possible phases, and assigned the variant to the chromosome producing the larger t^c. tc values were estimated separately for each population, not collectively as a single pooled sample. We also estimated tc for singleton variants in the mapping sample of 3,621 individuals (7,242 genomes) from the UK10k data set that has been filtered by the UK10k consortium to remove close relatives and individuals of recent non-European ancestry [16]. These include genomes from the ALSPAC cohort, which focused on the Avon region, and the TWINSUK cohort which includes samples from across the UK. Haplotype phase for variants found two or more times was inferred by the UK10k consortium using SHAPEIT. Haplotype phase for singleton variants was determined as for the 1KGP data. For all data we masked TpG/CpG transitions to minimize mutation rate variability. We compared results found using genetic maps based on both linkage disequilibrium [26] and pedigrees [31], and found very little difference in tc estimates (S8 Fig). Results presented here used the pedigree-based map. We assume the (non-CpG) mutation rate is constant and equal to 1 × 10−8 per base per generation [32–34]. We assume a demographic history of constant population size of N = 1 × 104. Data were simulated using the msprime program [35] under three demographic models: a constant population of size N = 1 × 104; a population with recent exponential growth that started from an ancestral population size of Na = 1 × 104 to a population size of N = 5 × 105 over the last 200 generations prior to sampling; and a bottleneck model in which samples were drawn from a ‘European’ population included in an out-of-Africa model with exponential growth and migration using the parameter estimates of Gutenkunst et al. , [24]. We simulated samples of 1000 10-Megabases chromosomes from which we sub-sampled populations of 100 chromosomes for comparison with 1KGP data. We used a fixed mutation rate of μ = 1 × 10−8 per base per generation, as is typical for human populations [32–34]. For recombination rate we used a value of ρ = 1 × 10−8 that is typical of the mean rate per base pair in human populations [26,31], as well as both a lower recombination rate (ρ = 1 × 10−9) and a higher recombination rate (ρ = 1 × 10−7) to assess performance when a different proportion of msh tracts end in a variant that has been introduced to the focal chromosome or a sister chromosome by recombination, and not directly by mutation. We assess error as the root mean squared error of log10 transformed tc values, and bias as the average signed error of the log10 transformed tc values. To assess the impact of phase uncertainty we randomly paired chromosomes and phased the data using SHAPEIT with default parameters. Singletons were phased as for the 1KGP and UK10K data. We used SNPEff [36] to identify variants annotated as missense, stop gained, stop lost, start lost, splice acceptor, or splice donor as ‘protein changing’. For singleton variants in UK10K populations, a variant is labeled as ‘private’ when it is found in only one population, and ‘shared’ when it is also found in other populations (at any frequency). Singleton variant ages in the 1KGP data were analyzed with a type-II ANOVA that examined the effect of three variables, including: (1) the population in which the variant was observed; (2) whether or not a variant is private to a single populations; and (3) whether or not a variant effects a protein sequence. The model has the form log10 (tc) ~ Intercept + PC + Private + Population + (PC × Population) + (Private × Population) + (PC × Private), where PC is a categorical variable indicating protein-changing status, Private is a categorical variable indicating the restriction of a variant to a single population, and Population is a categorical variable indicating the population in which a variant is observed. The only human data used in this study are genome sequences that are publicly available. Two sources of human data were used. The 1000 Genomes data were downloaded from the data server at http: //www. internationalgenome. org/data. The UK10K data are made available to researchers as described (https: //www. uk10k. org/data_access. html) and were made available to Dr. Hey following completion of the UK10K Project Data Access Agreement (signed 14 Oct 2014). These data included dataset IDS EGAS00001000090 (UK10K COHORT ALSPAC) and EGAS00001000108 (UK10K COHORT TWINSUK). These data include no information such that the original subjects could be identified. Because all of the data is publicly available, and because none of the data can be used to identify the original subjects, the study falls under NIH Human Subjects Research Exemption 4. For these same reasons the research was determined to not involve Human Subjects by Temple University IRB.
The program, entitled “runtc”, that implements the estimator for VSF files is available at https: //github. com/jaredgk/msh-python/tree/master/msh_est. The simulation results are available at https: //bio. cst. temple. edu/~hey/nolinks/Platt_etal_SimulationsAndAnalyses. zip.
We developed a way to estimate the time when a copy of a gene most recently shared ancestry with other copies of that gene. This is also an estimate of the upper bound of when a mutation has arisen, and it can be used to study the ages of alleles that are found in a population. The method can be applied to the very rarest alleles found only once in a sample, even in studies of many thousands of genomes. We tested the method extensively, found it performs well, and can be used under a wide variety of conditions. We applied it to 1000 Genomes project data (26 populations) and the UK10K data (over 7000 genomes) and found clear evidence that alleles that change proteins are younger than alleles that do not, as expected. We also observed wide variation in the ages of alleles at low frequency among the 1000 Genome project populations, indicating that our method could be used to study the demographic history of human populations. Going forward, the estimator should be useful for many kinds of questions in population genomics, particularly as sample sizes continue to grow.
alleles
genetic mapping
simulation and modeling
mutation
molecular biology techniques
research and analysis methods
sequence analysis
gene mapping
bioinformatics
sequence alignment
chromosome biology
molecular biology
genetic loci
haplotypes
cell biology
natural selection
heredity
database and informatics methods
genetics
biology and life sciences
evolutionary biology
evolutionary processes
chromosomes
An estimator of first coalescent time reveals selection on young variants and large heterogeneity in rare allele ages among human populations