Dating genomic variants and shared ancestry in population-scale sequencing data

1Each generation, a human genome acquires an average of about 70 single-nucleotide changes through mutation in the germline of its parents. Yet while, at a global scale, many millions of new variants are generated each year, the vast majority are lost rapidly through genetic drift and purifying selection. Consequently, even though the majority of variants themselves are extremely rare, the majority of genetic differences between genomes result from variants found at global frequencies of 1% or more, which may have appeared thousands of generations ago. Genome sequencing studies have catalogued the vast majority of common variation (estimated to be about 10 million variants), and, at least within coding regions and particular ancestries, to date, more than 660 million variants genome-wide have been reported, many of them at extremely low frequency.

2Despite the importance of genetic variation in influencing quantitative traits and risk for disease, as well as providing the raw material on which natural selection can act, relatively little attention has been paid to inferring the evolutionary history of the variants themselves, with notable exceptions of evolutionary importance, particularly those affecting geographically varying traits such as skin pigmentation, diet, and immunity. Rather, attention has focused on the indirect use of genetic variation to detect population structure, identify related samples, infer parameters of models of human demographic history, and estimate the time of ancestral population divergence. The evolutionary history of classes of variants contributing to polygenic adaptation (for example, those affecting height) or causing potential loss of gene function has received attention, though rarely at the level of specific variants. Previous work on rare variants has identified ancestral connections between individuals and populations and demonstrated evidence for explosive population growth. Nevertheless, to date, no comprehensive effort has been made to infer the age, place of origin, or pattern of spread for the vast majority of variants.

3We have developed a method for estimating the age of genetic variants; that is, the time of origin of an allele through mutation at a single locus. Our approach, which we refer to as the Genealogical Estimation of Variant Age (GEVA), is similar to existing methods that involve coalescent modeling to infer the time to the most recent common ancestor (TMRCA) between individual genomes. However, these methods typically operate on a discretized timescale, utilize only a fraction of the information available in larger sample data, or employ approximations to overcome computational complexity. Alternate approaches, in particular those that have been used to indicate allele age, are often strongly parameterized and make assumptions about demographic history, selection, or the shape of genealogical trees. Recently proposed approaches for genealogical inference are also able to estimate allele age by placing mutation events on the branches of reconstructed trees at variable sites. Our work is motivated by the desire for a fast and scalable nonparametric estimator of allele age that makes no assumptions about the demographic or selective processes that shaped the underlying genealogy and that is robust to the frequency and types of error found in modern whole-genome population sequencing studies. We learn about the age of a mutation by combining probabilistic distributions of the TMRCA between hundreds or thousands of genomes on a continuous timescale, irrespective of contemporary allelic distributions or historic population boundaries and without the need to fully reconstruct genealogies.

4The approach used within GEVA is outlined in Fig 1. A copy of the piece of the ancestral chromosome on which the mutation occurred is still present today in the individuals carrying the derived allele (Fig 1A). Over time, additional mutations have accumulated along the inherited sequence (haplotype), and its length has been broken down by recombination during meiosis in each generation. We locate this ancestral segment (often referred to as an identity-by-descent [IBD] region) relative to the position of a given target variant using a hidden Markov model (HMM) constructed empirically from sequencing data to provide robustness to realistic rates of sequencing and genotyping error (Fig 1B). By measuring the impact of mutation and recombination on the segments shared between pairs of haplotypes, we infer the TMRCA using probabilistic models to accommodate the stochastic nature of mutation and recombination processes (Fig 1C). Moreover, we make full use of the information available in whole-genome sequencing data to perform comparisons between pairs of chromosomes that both carry the mutation (concordant pairs) and pairs in which one carries the mutation and the other carries the ancestral allele (discordant pairs), thereby considering genealogical relationships that are both younger and older, respectively, than the time of mutation. Information from hundreds or thousands of haplotype pairs is then combined within a composite-likelihood framework to obtain an approximate posterior distribution on the age of the derived allele (Fig 1D). One benefit of our method is that we can increase the number of pairwise TMRCA inferences incrementally to update age estimates or to combine information across many data sources to improve the genealogical resolution from a wider distribution of independently sampled genomes. We additionally use a heuristic method for rejecting outlier pairs to improve robustness to low rates of data error and recurrent mutation. Full details are given in S1 Text.

5To validate GEVA, we performed coalescent simulations under different demographic models; see Fig 2. Using a standard coalescent model with constant mutation and recombination rates, we found low bias (relative error, ε = 0.268; see S2 Text) for allele age estimates and high correlation between true and inferred age (Spearman's ρ = 0.953; Fig 2A). We compared our approach for estimating the TMRCA to the computationally more demanding pairwise sequentially Markovian coalescent (PSMC) methodology, which forms the basis of many applications in ancestral inference. PSMC estimates a model of the demographic history between pairs of chromosomes (over a discretized grid of time intervals) and can, for every position in the genome, return the inferred posterior distribution on the TMRCA, thus enabling a composite-likelihood estimator of allele age as in GEVA.

6We found that PSMC-based age estimations performed similarly well to GEVA (ρ = 0.952), though the time discretization increased bias (ε = 0.530) and, in particular, led to overestimation of the age for the youngest variants. We note that PSMC was not designed strictly for this purpose and hence is not optimized for estimating allele age. Conversely, pairwise estimates of the TMRCA between concordant haplotypes were highly correlated with true TMRCA in both GEVA (ρ = 0.922) and PSMC (ρ = 0.919), but the correlation for discordant pairs was lower in GEVA (ρ = 0.586) compared to PSMC (ρ = 0.766); see S1 Fig. Such differences in relation to estimating allele age with high accuracy are tolerated because the time of mutation is estimated from the composite distribution of TMRCA posteriors from many pairwise comparisons performed at a single locus. Events that occurred distant in time (relative to the time of mutation) will have little influence on the estimate.

7Under a complex demographic model that recapitulates the human expansion out of Africa and with empirical and variable recombination rates, age estimation in GEVA maintained a similarly high level of accuracy (ε = 0.198, ρ = 0.937; Fig 2B). In this situation, although PSMC modeled the more dynamic demographic histories between haplotypes with higher accuracy (correlation of true and inferred TMRCA for discordant pairs: ρ = 0.915) compared to GEVA (ρ = 0.775), the time discretization resulted in artifacts at more recent times (TMRCA for concordant pairs: ρ = 0.892) that were not present in GEVA (ρ = 0.932), leading to worse performance when estimating the age (ε = 0.409, ρ = 0.900; S2 Fig), with the addition of substantial computational cost. We next introduced realistic data complications by reproducing empirically estimated genotype errors in simulated data (see S2 Text), as well as errors arising through in silico haplotype phasing (Fig 2B). We found age estimation in GEVA to remain largely unbiased and strongly correlated with true age after the inclusion of data error (ε = 0.346, ρ = 0.925; S3 Fig) and after phasing (ε = 0.430, ρ = 0.921; S4 Fig). The PSMC-based approach continued to show higher bias and reduced correlation at the same set of variants, both after error (ε = 1.042, ρ = 0.882) and after phasing (ε = 1.009, ρ = 0.880). Reduced data quality resulting from sequencing errors may introduce false signals of pairwise differences seen between haplotypes, and phasing errors may lead to an underestimation of haplotype lengths at variants that are relatively young, for which we overestimate TMRCA and hence allele age (particularly for alleles younger than approximately 100 generations).

8To evaluate the performance of GEVA on empirical data, we first considered variants affecting the well-studied lactase persistence (LP) trait, for which numerous approaches, including the use of archaeological data, genetic data, and a biological understanding of the functional and evolutionary impact of previously associated variants, have resulted in consensus expectations for the age. The LCT gene encodes the lactase enzyme but is regulated by variants in an intron of the neighboring MCM6 gene (Minichromosome Maintenance Complex Component 6). We estimated the age of the derived T allele of the rs182549 variant (G/A-22018), which is at a frequency of approximately 50% in European populations and forms part of a haplotype associated with LP. Under a model that jointly considers mutational and recombinational information, we estimated the allele to be 688 generations old (Fig 3A), originating approximately 14,000 to 21,000 years ago, depending on assumptions about generation time in humans. Our estimate is based on data from two different sources, the 1000 Genomes Project (TGP) and the Simons Genome Diversity Project (SGDP), which, when estimated separately, give very similar ages (692 and 687 generations, respectively). The full result data set for this variant is available online: https://human.genome.dating/snp/rs182549. We obtained a similar age estimate of 693 generations for the derived A allele of the rs4988235 (C/T-13910) variant (see https://human.genome.dating/snp/rs4988235), which is also strongly associated with LP and in near perfect association with rs182549, though we note that there is evidence for multiple origins of the variant. Previous estimates of the age of these variants range between 2,200 and 21,000 years, putting our estimates on the higher end of this range. Multiple sources of information suggest that these variants only achieved high frequency in European populations within the last 10,000 years (approximately 400 generations) and that LP alleles were rare until the advent of dairy farming in Europe. Our results therefore suggest that the mutation conferring the strongly selected phenotype (estimated to have a selection coefficient of up to 15% in European and up to 19% in Scandinavian populations) was present for hundreds of generations before its rapid sweep through the population.

9We next considered the protein-coding missense variant rs3827760 in the Ectodysplasin A Receptor (EDAR) gene, where the derived G allele (Val370Ala substitution) is found at high frequency in East Asian populations (87% in the TGP, 82% in the SGDP) and American populations (39% and 80%, respectively) and is associated with sweat, facial and body morphology, and hair phenotypes. We estimated the variant to be 1,456 generations old, approximately 29,000 to 44,000 years (Fig 3B; https://human.genome.dating/snp/rs3827760), again with strong concordance between the TGP (1,513 generations) and the SGDP (1,346 generations). Our estimate is consistent with previous estimates and limited evidence from ancient DNA studies. Our results further show that most individuals carrying the allele share a common ancestor close to the time the allele arose through mutation. This may suggest that the variant rapidly rose in frequency following its origin, which is consistent with previous findings of strong positive selection of this variant in East Asia. We compared this result to variants of similar age (1,450 ± 100 generations) in the Atlas of Variant Age (see further below): we found 1,043,376 variants dated in the TGP, of which 2,073 have reached a frequency higher than 30% globally and only 130 above 80% in East Asian populations, demonstrating how unusual such a rapid rise in frequency is.

10Finally, we considered the variant rs80194531, in which the derived allele causes an Asn78Thr substitution in the Zinc finger E-box–binding homeobox 1 (ZEB1) gene. The variant is reported as pathogenic for corneal dystrophy but is present at 6% in African ancestry samples within the TGP or the SGDP. We estimated the age of the variant to be 5,866 generations old (120,000 to 180,000 years), again with consistency between the TGP and the SGDP (5,879 and 5,854 generations, respectively; Fig 3C; https://human.genome.dating/snp/rs80194531). Such an ancient age seems inconsistent with the reported dominant pathogenic effect. Moreover, of the 1.2 million variants found at comparable frequencies (5%–7%) in African ancestry individuals within the TGP, we found that 46% were estimated to be younger than the rs80194531 allele, suggesting that this variant is in no way unusual.

11We next sought to characterize the age distribution of genetic variation across the human genome, for which we applied GEVA to more than 45 million variants identified in the TGP or the SGDP. More than 32 billion haplotype pairs were analyzed to estimate shared haplotype segments and TMRCAs. For variants present in both data sources (13.7 million), we additionally estimated the age by combining pairwise TMRCA distributions that were inferred independently in each sample after confirming that separately obtained age estimates agreed (Spearman’s ρ = 0.862; see S5 Fig). We make this information, referred to as the Atlas of Variant Age for the human genome, publicly available as an online database (https://human.genome.dating). A breakdown by chromosome of the number of variants dated and haplotype pairs analyzed is given in S1 Table. Further details are given in S3 Text.

12We find substantial variation in the relationship between estimated age and allele frequency, depending on the population in which frequency is measured and the geographical distribution of the variant (Fig 4A). Variants in African ancestry groups are typically older than in other groups and also have the greatest variance in age for a given frequency. For example, variants below 0.5% (within a given population) have a median age of 670 generations in African ancestry groups, 377 generations in East Asian ancestry groups, and 488 generations in Europeans (see S2 Table). The age distribution of variants restricted to a particular ancestry group (or shared between them) indicates the degree of connection between populations. For example, there are many variants up to 5,000 generations old (100,000 to 150,000 years) that are restricted to African ancestry groups yet are observed at frequencies up to 10%, but variants in this frequency range that are restricted to East Asian or South Asian ancestry groups are typically under 1,000 generations (20,000 to 30,000 years) or 1,300 generations old (26,000 to 39,000 years), respectively. Conversely, cosmopolitan variants that are shared among every ancestry group are typically older than 2,000 generations (40,000 to 60,000 years) despite being observed at global frequencies below 0.5% (S6 Fig). Variants restricted to American ancestry groups are typically younger than 750 generations (15,000 to 22,500 years), consistent with existing knowledge about the settlement of the Americas via the Bering land bridge that connected Asia and North America during the last glacial maximum around 15,000 to 23,000 years ago. We note, however, that recent admixture and the sampling strategies of the different data sets can have a strong impact on age distributions. For example, variants at high frequency within American populations but that are nevertheless restricted to just American and African populations are, on average, younger than lower-frequency variants (within American populations) with the same geographical restriction (S6 Fig). These variants likely arose recently within Africa and entered American populations through admixture, rising to high frequency through population bottlenecks. Similarly, variants found only within European and African populations but that have a frequency below 0.5% in Europeans are, on average, older than variants observed at higher frequencies in Europeans and older than variants restricted to only Europeans in the same frequency range, suggesting recent gene flow from Africa into Europe.

13Such heterogeneity in the relationship between allele age and frequency, coupled with heterogeneous and unknown sampling strategies, complicates the use of frequency as a means of assessing variants for potential pathogenicity during the interpretation of individual genomes. The Atlas of Variant Age potentially offers a more direct approach for screening variants, given the high probability of elimination of nonrecessive deleterious variants within a few generations. To assess the value of allele age in the interpretation of potentially pathogenic variants, we estimated the ages of variants that had effects predicted as damaging by Polymorphism Phenotyping v2 software (PolyPhen-2) or deleterious by Sorting Intolerant From Tolerant software (SIFT) in the TGP (Fig 4B). Of the approximately 70,000 variants analyzed, 50% of damaging and 49% of deleterious variants were estimated to have arisen within the last 500 generations (10,000 to 15,000 years), compared to 41% of benign (PolyPhen-2) and 42% of tolerated (SIFT) variants (S7 Fig). Compared with control sets of variants (those annotated as benign or tolerated and matched for allele frequency within the focal ancestry group), variants annotated as damaging or deleterious had a notable dearth of older variants (>1,000 generations) for a given frequency, consistent with theoretical expectations and previous findings. Our results suggest that old alleles can largely be excluded from consideration of pathology (though recent origin is not evidence in favor of pathogenicity).

14Finally, we investigated the extent to which patterns of sharing of variants of different ages could power approaches for learning about genealogical history. Previous work has highlighted the descriptive value of genetic variants in identifying individuals who share recent common ancestry and patterns of demographic isolation and migration, though it has also highlighted the challenges of interpreting the output of approaches such as principal component analysis (PCA). Conversely, numerous model-based approaches have been developed that use patterns of variant and haplotype sharing to infer underlying demographic parameters, though these typically make strong simplifying assumptions about the space of possible histories.

15Here, we present a nonparametric approach for combining descriptive and inferential approaches to learn about ancestral connections between individual genomes (and groups of individuals) based on variant age information. Unlike existing methods that assign time-invariant ancestry proportions to individual genomes by reference to contemporary populations, we can estimate the fraction a given genome shares because of common ancestry with any other genome at different points in time, referred to as the cumulative coalescent function (CCF). We use a fast dynamic-programming approach to estimate a maximum likelihood CCF between any pair or group of individuals (Fig 5A; see S4 Text). Using simulated data, we found that inference of ancestral relationships between individuals (CCFs inferred from estimated allele ages) correctly reflected differences of relatedness among individuals within and between different ancestry groups and revealed patterns qualitatively consistent with past demographic events (see S5 Text), though we note that uncertainty in variant age estimates may cause oversmoothing of coalescent profiles, in particular in the distant past (>20,000 generations ago).

16To illustrate the value of this nonparametric approach in describing the history of individuals and groups, we first considered the coalescent history between a single individual of American (Puerto Rican) ancestry from the TGP (ID: HG00733) and all others in the TGP sample, using GEVA age estimates for variants on chromosome 5 (Fig 5B). As a positive control, we included the parents of HG00733 (HG00732 and HG00731), who reach a CCF of near 1 in the most recent epoch (though note that the parents were used for haplotype phasing, which estimates transmitted haplotypes; hence the CCF reaching 1 rather than the expected one-half). We show the ancestry of HG00733 shared with every individual in the TGP sample in S1 Movie; also see https://human.genome.dating/ancestry/HG00733. Within the first 100 generations, we see additional coalescence with the untransmitted parental chromosomes and other individuals from the Puerto Rican sample. The earliest common ancestry outside Puerto Rico is seen with a European individual in Spain (paternal side; approximately 60 generations ago) and a Mexican ancestry individual (maternal side; approximately 80 generations ago). Coalescence with individuals sampled from outside the Americas occurs further back in time (>100 generations ago), initially with European individuals (predominantly around 300 to 600 generations ago), then uniformly with non-African individuals around 1,500 to 4,000 generations ago, more strongly with African individuals around 5,000 to 15,000 generations ago, and uniformly with all individuals around 20,000 generations ago. Because of the impact of data errors on age estimation of recent variants (highlighted above), the absolute timings of the early events are likely substantially overestimated. However, we expect the relative ordering of events to be robust and consistent among sample comparisons (see S5 Text).

17The CCFs to all other members of a reference panel (averaged across all chromosomes in both haploid genomes) provide an overview of the genealogical relationships for a target individual. As an example, we inferred the CCF profiles of a Siberian Eskimo to all other individuals in the SGDP (Fig 5C, top), showing common ancestry to other Central Asian and Siberian individuals within a few hundred generations, substantial common ancestry with American individuals before 1,000 generations, and typically more recent common ancestry with East Asians than West Eurasians or Africans; see https://human.genome.dating/ancestry/S_Eskimo_Sireniki-1. Notably, relatively little additional coalescence is seen during a period around 2,000 to 10,000 generations ago, which is a pattern shared among non-African individuals (also see S1 Movie) and agrees with previous findings of a period of reduced coalescence, peaking 100,000 to 200,000 years ago.

18The CCF can also be represented as a coalescent intensity function (CIF; see S4 Text), which measures the rate of change of common ancestry over time (Fig 5C, middle). The CIF reveals additional structure; for example, around 3,000 to 20,000 generations ago, those parts of the Siberian Eskimo's genome that have not yet coalesced with other genomes sampled from the same ancestry group have a very low CIF, while the CIF to the African ancestry samples (which have had very little coalescence until this point) is relatively high (though note the absolute rate remains very low over this period).

19We can further summarize the coalescent profiles for the target individual by computing the maximum CIF over the sample. This statistic captures properties analogous to the effective population size parameter, Ne, in population genetics modeling, in which the expected coalescent intensity is inversely related to population size. We refer to the maximum CIF in the following as the Ne equivalent (Fig 5C, bottom), though we note that (in the case in which there are genuine populations) it is likely to be downward biased compared to existing methods; for example, we would expect PSMC to infer absolute values of ancestral population size more accurately. We therefore use Ne equivalents to provide a relative summary of genealogical histories across the entire cohort.

20We estimated all pairwise CIFs in TGP (S2 Movie), which we aggregated across autosomes to generate a coalescent profile between each pair of individuals (see S6 Text). Likewise, we estimated the ancestry shared between each individual pair in SGDP (S3 Movie); CIFs were further aggregated among the 130 population groups (Fig 6). These reveal how the rates and structure of coalescence have changed over time, with the most recent epoch (up to approximately 200 generations ago) dominated by coalescence within populations, but also identify recent connections between groups such as between southern Siberian and northern East Asian populations (Fig 6A). Several populations such as the Kusunda (Nepal), Saami (Finland), and Negev Bedouins (Israel) show strong within-group coalescence up to this point, though by 500 generations ago they are coalescing primarily with other populations. The epoch around 800 generations ago is dominated by structure broadly corresponding to the continental level (Fig 6B), though some African populations—for example, the Mbuti (Congo), but more dramatically the Khomani San (South Africa) and Ju'hoansi (Namibia)—remain isolated up to approximately 1,500 generations ago (30,000 to 45,000 years), which overlaps with previous findings, and we see these two populations to be strongly connected for an extended period further back in time (S3 Movie). Around 800 generations ago, there is very little remaining structure among West Eurasian populations, but many additional intercontinental connections are now identified. For example, we see a north-to-south gradient of decreasing coalescence between American populations and Siberian or East Asian populations. In particular, we identify strong coalescence of American ancestry individuals with Siberian Eskimos, Aleutian Islanders, and Tlingit people in a period between 500 and 1,000 generations ago and very little structure among American, Siberian, and East Asian populations as a whole further back than 1,000 generations ago (S3 Movie), which agrees with previous results regarding the human migration into the Americas, extended isolation, and subsequent dispersal across the continent. By 4,000 generations ago, we see high levels of coalescence between non-African and African populations (Fig 6C) and essentially no structure in the epoch around 20,000 generations ago (Fig 6D). We note that the exact timings of demographic events (signified by periods of intense or reduced coalescence), while showing consistent patterns within and among population groups, may carry additional noise due to uncertainty in allele age estimates.

21The maximum CIF profiles (Ne equivalents; Fig 6E) highlight several features, including differences among modern ancestry groups in the intensity of coalescence within the last 1,000 generations (particularly intense for American and Oceanic populations); a major period of intense coalescence among all non-African ancestry individuals around 1,000 to 2,000 generations ago, following the migration of modern humans out of Africa; a weaker, but still marked, increase in coalescent intensity for African ancestry samples around 2,000 generations ago; and an older reduction in coalescent intensity, peaking around 5,000 to 8,000 generations ago, potentially driven by ancient population structure within Africa and (for non-African populations) possible admixture with archaic lineages. We find minor quantitative, but not qualitative, differences among chromosomes (S8 Fig). We note that the bottleneck and steep change in population size around 60,000 to 80,000 generations ago, although consistent across ancestry profiles inferred for all individuals in the sample, is likely to be an artifact possibly caused by decreased accuracy of age estimates in the distant past, along with a breakdown of assumptions regarding the statistical independence of very old mutations.

22We have demonstrated how allele age estimates can provide insight to a range of problems in statistical and population genetics. However, there are several important assumptions and limitations of the approach. First, a key assumption is that of a single origin for each allele. Given the size of the human population and the mutation rate, it is likely that every allele has arisen multiple times over evolutionary history. Yet, unless the mutation rate is extremely high, it is still probable that most individuals carrying the allele do so through common ancestry. Moreover, multiple origins can potentially be identified through the presence of the allele on multiple haplotype backgrounds, as has, for example, been seen for the rs4988235 variant at LCT (though we note that a previous study concluded that the allele of variant rs4988235 was brought into African populations through historic gene flow, possibly through the Roman Empire), the O blood group, or alleles in the Human Leukocyte Antigene (HLA) region. A variant lying in a region with high rates of noncrossover (gene conversion) may similarly be found on multiple haplotype backgrounds. However, for genomes with very high mutation rates, such as HIV-1, recurrence is sufficiently high to make estimates of allele age meaningless. In addition, while we have shown GEVA to be robust to realistic levels of sequencing and haplotype phasing error, the actual structures of error found in reference data sources, such as the TGP, have additional complexity whose effect is unknown.

23Our approach also assumes a known and time-invariant rate of recombination. For most species, only indirect estimates of the per-generation recombination rate are available, and in humans and mice, there is evidence for evolution in the fine-scale location of recombination hotspots through changes in the binding preferences of PRDM9 (PR domain zinc finger protein 9). However, because broad-scale recombination rates evolve at a much lower rate than hotspot location and because our approach for detecting recombination events is driven largely by the presence of recombinant haplotypes, we expect GEVA to be relatively robust for recent variants. Older variants may be more affected, but for such variants, most information comes from the mutation clock, which is likely to have been more stable over time.

24We have shown that the ages of variants are highly correlated when estimated separately in independent samples, which agrees with the assumption that, for the majority of alleles, mutations occurred only once in the history of the population, such that the true age of an allele refers to a fixed point in time. In principle, the age can be estimated irrespective of the frequency distribution observed in different study cohorts. Our results show substantial heterogeneity in the relationship between frequency and age, revealing the impact of often unknown demographic variables on the distribution of alleles within and across different ancestry groups, which may further be confounded by the mode and strength of selection on particular alleles. The estimation profiles of variants shown in Fig 3 provide only a few examples of cases in which the frequency may not be reliable as a proxy for the age. Our method and the Atlas of Variant Age may therefore provide more robust measures for analyses that otherwise rely on population or sample-specific allele frequencies.

25The Atlas of Variant Age also has multiple applications beyond statistical and population genetics. For example, recent variants provide a natural index when searching for related samples in population-scale data sets. We have demonstrated how variant age information can be leveraged more comprehensively, using a nonparametric approach, to learn about genealogical history, relatedness between individuals, and ancestral connections between populations. Our approach to infer coalescent profiles (pairwise CCF and sample-wide CIF) captures ancestry proportions and coalescence rate variations as a function of time, from which we can distinguish recent from past demographic effects or estimate changes of relative relatedness over time. As illustrated for whole-genome sequencing data from the TGP and SGDP, our methodology accesses a time dimension to study relatedness between individual genomes and among entire sample cohorts. Applications to other data sets may use information from the Atlas of Variant Age without the need to re-estimate ages, thus demonstrating the value of this resource for exploratory data analysis.

26Finally, as presented here, it is possible to combine information from multiple, potentially even distributed data sets by estimating coalescent time distributions for pairs of concordant and discordant haplotypes in each data resource separately or to update age estimates by the inclusion of additional samples. Estimation of the age of particular variants may gain additional accuracy by increasing the number of pairwise comparisons, depending on the history of the allele and the genealogical resolution attainable from the diversity of ancestral backgrounds inherent to available sample data. Future studies may reveal new variants, particularly within under-represented ancestry groups, that can also be added to the Atlas of Variant Age. The methods are also applicable to nonhuman species as long as estimates of mutation rate, recombination rate, and generation time are available. Future extensions to infer location of origin or the ancestral haplotype, integrating the growing wealth of genome data from ancient samples, will be an important step towards reconstructing the ancestral history of the entire species.

27Estimation of allele age and shared ancestry was conducted on publicly available data sets: the TGP and the SGDP. We used phased haplotype data of chromosomes 1–22 from the final release TGP panel (Phase 3; GRCh37), available for 2,504 individuals from 26 populations worldwide (5 continental population groups). Additional data were available from the TGP for 31 related individuals, which we included in our shared ancestry analysis. We used phased haplotype data of chromosomes 1–22 from the publicly available SGDP panel (PS2; GRCh37), consisting of 278 individuals from 130 populations worldwide (7 continental population groups). Recombination rates were determined for each chromosome using the genetic maps available from the International HapMap Project (Phase 2; GRCh37). Genotype data from the Illumina Platinum Genomes Project (IPG) (GRCh37; chromosomes 1–22) were used as a reference to measure genotype error in a matched subsample from the TGP. We used information from the Ensembl database (GRCh37; release 92, version 20180221) to determine the ancestral and derived allelic states for variants in both the TGP and SGDP panels, as predicted through multispecies alignments in the Ensembl EPO pipeline.

28The Atlas of Variant Age is publicly available as an online database: https://human.genome.dating/. It contains the age estimation results for more than 45 million variants in the human genome, including the results of all pairwise analyses conducted at each focal site to estimate shared haplotype segments and infer TMRCA between haplotype pairs (more than 32 billion). We used an HMM to infer shared haplotype segments, for which we constructed an empirical error model (through comparison between genotypes from the IPG and TGP) to generate frequency-dependent emission and initial state probabilities; these are available online: https://github.com/pkalbers/geva/tree/master/hmm. Furthermore, we performed shared ancestry analyses for every pair of individual chromosomes in the TGP (including related individuals) and the SGDP, which are also available at https://human.genome.dating/.

29The source code of the GEVA method is available online at https://github.com/pkalbers/geva. The source code of the CCF dynamic programming algorithm is available at https://github.com/pkalbers/ccf. We modified the original source code of MSMC2 to optimize the performance of the PSMC algorithm in our simulation analysis, available at https://github.com/pkalbers/msmc2.