Demographic model and generation of simulated aDNA

Our general framework is to simulate genomic data using a three-population model:

• - Pop1: the reference haplotype population;

• - Pop2: the target population that will be imputed (i.e., the “ancient population”) and is diverged d generations from Pop1;

• - Pop3: an outgroup population with fixed divergence from Pop1 and Pop2 (Figure 1).

Figure 1 Population demographic model of simulated genomes.

We note that we do not explicitly model an “ancient” Pop2 population that is sampled sometime earlier than Pop1. Instead, we are interested in the effect of the degree of drift between Pop1 and Pop2 on imputation accuracy. Thus, in our model, if present day Pop1 and present day Pop2 diverged 50 generations ago, this would be the equivalent effect of the amount of drift that would occur between present day Pop1 and an ancient Pop2 sampled 50 x 2 = 100 generations ago under a scenario of complete population continuity. This greatly simplifies the parameter space that we are required to explore, as we do not need to consider to what extent ancient Pop2 might have diverged from the lineage from which Pop1 descended before being sampled (of which there are essentially unlimited possibilities). For example, this amount of drift would also be equivalent to both an ancient population sampled 60 generations ago and that diverged from the ancestral population leading to Pop1 80 generations ago, or one sampled 40 generations ago and that diverged 70 years ago.

In order to generate our baseline genomic data (i.e., the truth set), we used fastsimcoal26 [30] to simulate polymorphic haploid genomes roughly equivalent to human chr1 (~220 Mb) using a three-population coalescent model populations under various parameters (Figure 1, Table 1). Effective population size was set to 30,000 for each population. Per site per generation mutation rate was set to 2 x 10^–9, and recombination rates were based on the HapMap GRCh37 genetic maps for chr1. Simulated haploid genomes were randomly paired within populations to make diploid individuals. Due to the computational burden when simulating high recombination rates in fastsimcoal26, regions with rates on the order of 1 x 10^–7 or above were set to 1 x 10^–7. Simulations were performed under five different divergence (d) times for Pop1 and Pop2: d = 50, 100, 500, 1500 and 2000 generations (the last being roughly equivalent to human non-African/African population divergence ~50 thousand years ago assuming 25 years per generation). Pop3 has a fixed divergence of 3,000 generations from the ancestor of Pop1 and Pop2. Our simulations generated ~100,000 biallelic variants above 5% frequency in the reference population, which is within the range expected for current human paleogenomic applications (~400,000 SNPs >5% frequency are found in European populations from the 1KGP based on whole genome shotgun sequencing [31] and 60,000 when filtering down to the most commonly applied 1240K in-solution capture SNP set [6,7]).

Table 1 Simulation parameters.

Furthermore, in order test the impact of a reference panel made up of haplotypes of varying divergence times, we simulated a scenario where the reference panel consisted of N = 1,000 individuals from Pop 1 that diverged, d = 2000 generations from Pop2, while also including five additional individuals that diverged d = 25 generations from Pop2.

In addition, to explore a more complex demographic scenario, we simulated a population expansion in the target population characteristic of the demographic patterns of Mid-Holocene European farmer populations. Using parameters estimated from Veeramah et al. [32] we modeled a population that diverged from the reference population 250 generations ago (~6 Ky) and expanded at a rate of 2% per generation. The target population (Pop2) was sampled 200 generations in the past, resulting in a total of 50 generations of divergence between the target sample and the reference panel (d = 50).

Following simulation by fastsimcoal26, the data from our model were then partitioned and converted into three types of observable datasets: (a) true phased genomes from Pop1 and Pop3 for use as the reference panel for imputation and downstream population genetic analysis; (b) true phased genomes from Pop2; and (c) unphased, low coverage aDNA-like genomes from Pop2. Resulting Arlequin files were converted to VCF files with custom python scripts for all three populations. Low coverage datasets for Pop2 with 0.1x, 0.5x, 1x, and 5x mean coverage were generated from the original true genotypes by first drawing the total number of reads for each site from a random Poisson distribution with γ = 0.1, 0.5, 1, 5, respectively. The number of reads from each allele at true heterozygous sites was then drawn from a binomial distribution with p = 0.5. Finally, to account for sequencing errors, the cycle position of each read was drawn from a random uniform distribution on the interval 1–81 (81 being the total read length and noting that we do not explicitly model variable ancient fragment size) and Phred quality scores were simulated based on an 8th century Viking genome (Figure S1) [22]. We chose not to explicitly introduce aDNA damage, such as deaminations, as the several potential ways to introduce post-mortem DNA damage would have significantly increase the complexity of our data generation, while most common uses cases in paleogenomics largely mitigate such effects through UDG treatment, terminal base clipping, and/or conditioning analyses to transversion sites [2]. SNP calling and genotype likelihood estimation were performed using custom scripts using the expressions described by DePristo et al. [33]. Non-biallelic sites (predominantly triallelic) were removed from any downstream analysis. We note that if a user were to utilize data with aDNA damage, a useful strategy would be to estimate genotype likelihood using methods that incorporate patterns of DNA damage, such as ATLAS [34].

Imputation strategies

In order to assess how different imputation strategies affect genotype calling accuracy, we examined a number of different of analytical pipelines to impute the ancient-like Pop2. Our strategies differed mainly in the size of reference panels, and choice of imputation software.

Traditionally, most imputation studies for ancient DNA have used large reference panels, such as the 1000 Genomes Project [31]. It has been suggested that rare variants are better imputed by increasing reference size with more diverse populations [35]. By using different reference panel sizes, we aim to shed light on how performance is affected by increasing reference size without necessarily introducing other populations. The reference population sizes used were N = 25, N = 100 and N = 1,000. Additionally, to circumvent the issues that arise with low coverages at fixed genotype calls, we use genotype likelihoods as input to BEAGLE and GLIMPSE, two of the most popular methods for aDNA imputation.

Imputation pipeline

For every tested combination of imputation strategies, divergence times, and coverages, we used genotype likelihoods as input to BEAGLE and GLIMPSE. Phased Pop1 was always used as a reference panel. We imputed the equivalent of chr1 using default imputation window size as defined in BEAGLE v4 [17] and GLIMPSE 1.1.0 [18]. For both software, burn-in iterations were set to 5 (default), and phasing and imputation iterations were set to 10 to increase imputation and phasing accuracy. All other parameters were set to default. Sites with a minor allele frequency (MAF) below 5% in the reference population were removed before imputation, as lower frequency variants are consistently shown to be poorly imputed due to their scarcity in reference panels, which tend to be highly enriched for higher frequency variants. Custom-made genetic maps based on HapMap GRCh37 recombination maps with capping at 1 x 10^–7 as described above were used for GLIMPSE imputation, while BEAGLE v4 does not utilize genetic maps. We also briefly compared the results of GLIMPSE 1.1.0 to GLIMPSE v2 [36] with N = 1000, as the latter is optimized for large reference panels with more than 2000 reference haplotypes.

Performance evaluation

For different combinations of reference and target datasets, we evaluated imputation accuracy by measuring the proportion of matching genotypes between the imputed and the true datasets for Pop2. We also measured accuracy within the context of MAF in the reference population and genotype probabilities (GP). Unless otherwise stated, results are shown for the most recent (d = 50) and the most divergent (d = 2000) reference panels made up of N = 1000 samples to evaluate the accuracy of imputation at our lower and upper bounds. We also examined the potential effect of imputation strategies on population genetic inference using principal component analysis (PCA) via smartpca [10,37], with Pop2 true and imputed genotypes projected onto PC space determined from Pop1 and Pop3 with the parameters lsqproject: YES and no outlier removal (outliermode: 2). Euclidean distance between true and downsampled target individuals was calculated using PC1 and PC2, and scaling the latter by the relative size of the eigenvectors (i.e., percentage of variation explained by each PC). We also explored the accuracy of phasing haplotype in the target Pop2 individuals with the vcftools [38] function diff-switch-error, which assesses the switch error between the true genotypes and those phased by BEAGLE and GLIMPSE.

Imputation accuracy by reference panel

Figure 2a and Figure S2a summarize the overall proportion of correctly imputed genotypes (our measure of accuracy) in the target population, Pop2, across a range of mean depth of coverages, divergence times (d) from Pop1, and reference population sample sizes (N). In addition to overall concordance, we calculated accuracy separately for true heterozygous sites (Figure 2b, Figure S2b). For samples with a mean depth of coverage of 0.1x overall accuracy for all reference panel sizes and divergence times is within the range of 60%–67%. From 0.5x onwards, the divergence times have a greater effect on accuracy, with more recently diverged target samples being better imputed, as expected due to longer shared IBD stretches between the target and reference samples. This effect is amplified with increase in the reference panel size and coverage, such that the highest accuracies are observed for the smallest divergence (d = 50), largest reference size (N = 1,000) and highest coverage (5x).

Figure 2 (a) Overall proportion of correctly imputed genotypes; (b) proportion of correctly imputed heterozygous genotypes.

As accuracy at homozygous sites is expected to be higher than at heterozygous sites (especially at lower coverages where often only one sequencing read might be sampled), we also conditioned accuracy at only true heterozygous sites. Heterozygous accuracy is substantially lower than overall accuracy for all combinations of variables, with only the single combination of d = 50, N = 1000 and 5x via GLIMPSE providing somewhat comparable results. At 0.1x, all target sample genotypes have an accuracy of ~40%, regardless of reference population size and divergence time. For 0.5x and above accuracy again predictably increases with coverage and reference panel and decreases with divergence.

Both overall and heterozygous accuracies are similar for d = 1500 compared to d = 2000, perhaps suggesting a saturation point. In addition, when considering all sites, at depths of coverage below 5x, BEAGLE slightly outperforms GLIMPSE, although both software perform similarly. However, at heterozygous sites GLIMPSE is consistently more accurate, especially with increasing coverage and smaller d.

We also note that for all coverages ≥ 0.5x, overall accuracy and accuracy at true heterozygote sizes is consistently higher for small reference sizes and divergence times (N = 25, d = 50) than for a large reference size and divergence time (N = 1,000, d = 2,000). Similarly, while increases in coverage consistently improve accuracy, the more diverged the reference and target, the slower the rate of increase in accuracy that accompanies increases in coverage (across the range of N tested accuracy generally increases by ~15% from 0.5x to 5x for d = 50, but less than 10% for d = 2000). Indeed, greater overall and heterozygous accuracy is consistently achieved when the coverage is 1x and d = 50 compared to when coverage is 5x and d = 1500 or 2000.

Additionally, we compared the performance of GLIMPSE v1 for N = 1,000 to a more recent version of GLIMPSE (v2) optimized for larger reference panels. Both versions of GLIMPSE performed similarly and reached similar levels of accuracy overall and at heterozygous sites (Figure S3).

Finally, we examined a scenario where the primary reference panel (N = 1000) for d = 2,000 was supplemented by five individuals d = 25 generations divergent from the target individual. Despite the addition of these more closely related haplotypes, we observed similar overall and heterozygous imputation accuracy when using a reference N = 1,000 without the individuals of varying divergence (Figure S4).

Imputation accuracy over the allele frequency spectrum

To examine imputation performance as a function of the allele frequency spectrum, we measured accuracy of imputation genotypes in three bins (0.05-0.1, 0.1-0.3, 0.3-0.5) based on the minor allele frequency (MAF) in the reference population (Figure 3). We report imputation accuracy for the least and most divergent target samples (d = 50 and d = 2000, respectively) to evaluate performance at both upper and lower bounds. We observed that heterozygote imputation accuracy is heavily dependent on MAF, with the best imputed variants being unsurprisingly those in the highest bin (Figure 3b). This effect is particularly pronounced with lower coverage, regardless of imputation method or d. At 0.1x the heterozygote accuracy is ~50% larger for the highest MAF bin (0.3–0.5 ~65%) compared to the least frequent bin (0.05–0.1 < 14%). As coverage increases, rarer variants are generally better imputed. Variants are generally better imputed at d = 50 compared to d = 2000, although the difference in accuracy between both divergence times is again much more evident at rarer variants. While GLIMPSE v1 and GLIMPSE v2 tend to perform similarly overall for N = 1000 for most parameter combinations, when imputing with more divergent reference panels (d = 2000) GLIMPSE v2 did outperform GLIMPSE v1 for the higher frequency variants bin (MAF 10%–50%) at 5x (74% vs. 65%) (Figure S5).

Figure 3 (a) Overall proportion of correctly imputed sites partitioned by MAF in a reference population of N = 1000 and d = 50 generations divergent; (b) proportion of correctly imputed heterozygous sites partitioned by MAF in a reference population of N = 1000 and d = 50 generations divergent.

When evaluating overall accuracy (Figure 3a), we somewhat surprisingly see the opposite pattern: rare variants are imputed correctly more frequently than common variants. A likely explanation for this seemingly counter intuitive observation is that higher MAFs in the reference panel allow for heterozygous sites in the target sample to be better discriminated from homozygous sites. However, when reference MAFs are low, sites could be essentially randomly imputed as homozygous for the major allele and be correct often simply by chance. Such concordance will affect overall imputation accuracy, given that rare variants represent roughly 30% of markers in the reference panel in our simulations (Figure S6).

Genotype probabilities

Once missing sites have been imputed, it is important to decide which sites will remain for downstream analyses. One possibility is to filter sites based on the average probability that a genotype call is correct via estimated genotype probabilities (GP) emitted by the imputation software and restrict analyses to variants with high certainty (e.g., >99%). We examined to what extent these emitted genotype probabilities reflected true error rates (e.g., are sites with a 99% genotype probability correctly imputed 99 times out of a 100?). To avoid the scenario of homozygous variants correctly imputed by chance observed above as a function of MAF, we restrict results to true heterozygous sites. Figure 4 summarizes the proportion of correctly imputed true heterozygous sites relative to the total number of sites imputed as heterozygous within genotype probability bins for d = 50 and d = 2000. While both methods seem to overestimate confidence in genotype calling at all population sizes and coverages, BEAGLE qualitatively tends to outperform GLIMPSE in assigning more representative probabilities for GPs greater than 60%. Although GLIMPSE correctly imputes slightly more heterozygous sites than BEAGLE (Figure 2), the probabilities associated with genotype calls are much less informative. We note that for d = 2000, the genotype probabilities emitted by GLIMPSE are particularly highly discordant (in the direction of being large overestimates) with genotype accuracy, regardless of coverage. While both versions of GLIMPSE were generally similarly discordant for N = 1000, a notable improvement was observed in assigning representative high (>95%) genotype probabilities when using GLIMPSE v2 with d = 2000 compared to GLIMPSE v1 (Figure S7).

Figure 4 Proportion of correctly imputed heterozygous genotypes partitioned by genotype probability scores (a) d = 50, N = 1000; (b) d = 2000, N = 1000.

Genotype phasing

Prior to imputation, haplotypes are resolved by assigning non-missing sites to a parental chromosome. As a result, imputation of missing variants will be affected by phasing accuracy. We evaluated phasing accuracy by measuring the proportion of ‘switches’ between the known maternal and paternal haplotypes from the true Pop2 dataset. In Figure 5, we show the proportion of correctly phased sites. Similar to imputation, accuracy is improved when using larger reference panels to phase higher coverage target samples, with a greater accuracy as divergence to the reference panel decreases, and phasing of our most divergent samples does not improve with reference population size, again suggesting a divergence threshold in which haplotype estimation does not improve.

Figure 5 Proportion of correctly phased genotypes based on switch-error. (a) Reference N = 25. (b) Reference N = 100. (c) Reference N = 1000.

Principal component analysis

In order explore how the kind of discordances between the true genotypes and the imputed genotypes might actually affect downstream analysis, we examined Pop2 imputed samples in a diploid PCA, projecting them on top of Pop1 and Pop3 populations that are 3000 generations diverged (Figure 1) and calculating the Euclidean distance between the true and imputed samples weighted by the relative percentage of variance explained by each principal component (Figure 6). Although we note that PCA can often be performed with pseudo-haploid calls, our analysis focuses on how reliable genotype imputation is for population genomic inference analyses. Figure 6 (outgroup Pop3 not shown for scaling) and Figure S8 show that higher coverage imputed genomes unsurprisingly map more closely to the actual true genotypes and demonstrate smaller variances. Similarly, imputed genomes generally map closer to the reference population as a function of decreasing coverage, with this difference becoming increasingly apparent with increasing d. Noticeably BEAGLE outperforms GLIMPSE for d ≥ 500 in terms of mapping closer to the true population and further from the reference population given the same coverage, but for d ≥ 100 shows the curious pattern of lower coverage samples mapping away from both the reference and true target population and moving towards the outgroup populations, and is thus outperformed by GLIMPSE. Interestingly, imputed samples do not seem to cluster more closely to the true samples when applying post-imputation genotype probabilities filters (GP > 0.99) (Figure S8A) and are systematically at greater Euclidean distances from the true sample than the non-GP filtered samples (Figure S8B). In addition, lower-coverage samples show a more widespread distribution across the PCA when imposing this filter. These observations suggest that genotype probability filtering discards a large number of sites, introducing more noise and without yielding substantial improvements for diploid-based inference.

Figure 6 Principal component analysis of imputed and true target samples (outgroup Pop3 not shown). (a) PCA projection; (b) Euclidean distance between true and imputed samples.

Imputation performance under a scenario of recent population expansion

While we cannot explore all possible demographic scenarios within the scope of this study, it is generally thought that Eurasian populations grew significantly within the Holocene era. Therefore, to gain some insight into how our results were affected by imposing constant population sizes, we also simulated a scenario in which the target haplotypes originated from a population that underwent population expansion, characteristic of the demographic patterns of Mid-Holocene European farmer populations. We modeled a population that diverged from the reference population 250 generations ago (~6 Ky) and expanded at a rate of 2% per generation [32]. The target population (Pop2) was then sampled 200 generations in the past, resulting in a total of 50 generations of divergence between the target sample and the reference panel (d = 50). Compared to the scenario of d = 50 with no population expansion, we observe a decrease of between 5% to 10% with regard to both overall accuracy of imputation and at true heterozygous sites (Figure S9). When examining this affect across the allele frequency spectrum it is notable that there is a much larger drop off in accuracy for variants in the lower MAF bin under a population expansion (0.05–1, with a difference in error rate of 22% and 24% for GLIMPSE and BEAGLE respectively for 5x) compared to the highest bin (MAF 0.3-0.5, 11% and 2% decrease for GLIMPSE and BEAGLE respectively) (Figure S10). While filtering for SNPs with a MAF > 5% in the reference population might be expected to mitigate most of the potential error introduced by an increase in newer low frequency variants in the expanding target population, there clearly are still residual effects on the haplotype patterns of more established lower frequency variants, likely due to distortion of patterns of linkage disequilibrium.

Effects of reference and target divergence

As would be expected, target samples that were more recently diverged from the reference population were better imputed overall. The greatest impact of increasing divergence on accuracy appears to be at heterozygous sites at low MAF in the reference population. This is in line with the work by Sousa da Mota et al. [29], where imputation of low-frequency variants was negatively impacted by the age and divergence between the African target samples and the reference panel. Interestingly, imputation accuracies for the most diverged target samples (i.e., oldest) we simulated (d = 1500 and d = 2000) were usually within the same range and behaved similarly across every imputation strategy and coverage, suggesting a possible divergence threshold in which imputation performance will not significantly improve with sequencing read depth. We also found that when divergence between the reference haplotypes and target genome was high (d ≥ 1500), this could significantly affect accuracy even when other conditions such as reference population size and coverage were optimum. For example, our results suggest that when assembling reference haplotypes sets, smaller numbers of samples that are more closely related to the target population to be imputed may be a more favorable strategy compared to compiling large numbers of very highly diverged samples. This was also the case even when integrating a few less divergent reference haplotypes into larger and more divergent reference panels. In addition, imputation may perform poorly even in samples with good (in the context of most existing ancient genomes) coverage if the reference sample is distantly diverged from the target (for example imputing an ancient Africa genome at ~5x using non-African reference haplotypes).

Generally, our results suggest that for humans, imputation with modern populations will likely produce usable results for historical-era and young prehistoric populations (for example those that straddle the Neolithic and Paleolithic era), but caution should be applied to the imputation of ancient genomes from the deep prehistoric era (~50,000 years ago), or that are deeply diverged from the reference populations, such as using non-African populations to impute African paleogenomes.

Even the addition of a small number of high coverage ancient genomes to the reference set may significantly improve imputation performance, particularly for inference of heterozygous genotypes with low MAF, as suggested by Ausmees et al. [27].

The impact of low coverage

While accuracy was benefited by increasing the size of reference panels, improvements in accuracy were less pronounced as coverage simultaneously decreased. Most strikingly, imputation performance for the lowest coverage samples did not considerably change with any of the imputation strategies we used. Samples with 0.1x coverage were consistently imputed with ~40% accuracy at heterozygous sites, or ~60% total accuracy, regardless of the imputation strategy or divergence between the reference and target datasets. This can be an important consideration, given that several ancient samples present low coverage (<1x). While our results suggest that imputation of missing data in very low coverage samples may not reach enough accuracy to perform population genetic analyses that require high SNP density, imposing post-imputation GP filtering has been suggested to be a viable option to discard a large number of incorrect calls. While our observations suggest this might be a viable option when imputing with BEAGLE, GLIMPSE filtering may require an alternative approach, since we found it greatly overestimates high GP calls, especially at lower coverages and high divergence. A possible way to circumvent this issue is that proposed by Hui et al. [28], where a pre-imputation GP filter is imposed, at the expense of discarding a large number of SNPs.

Beagle vs. Glimpse

Although both software perform similarly, BEAGLE’s overall accuracy is generally slightly higher than GLIMPSE. Accuracy at heterozygous sites may be a more interesting metric to evaluate imputation performance, given that they are less prone to random agreement and will affect overall imputation accuracy. GLIMPSE v1 and GLIMPSE v2 generally outperform BEAGLE at heterozygous sites, especially for more recent samples with higher coverages. When imputing lower frequency variants, GLIMPSE tends to outperform BEAGLE across every simulated coverage. The choice of which method to use may depend on the interest of the user. For example, one might prefer GLIMPSE if the reference and target datasets are not very distantly diverged. If filtering by GPs is of interest, one should consider that GLIMPSE v1 will likely keep a large number of sites with non-representative GPs. Even computational resources may be of consideration, as GLIMPSE constantly outperformed BEAGLE in terms of running time by several orders of magnitude, especially with larger reference panels [18].

Downstream analysis of imputed ancient genomes

When analyzing our imputed diploid genotypes using PCA, higher coverage samples showed closer mapping to true samples, while lower coverage samples clustered with the reference panel. The disparity between true and imputed samples increased with greater divergence times, especially with GLIMPSE imputation at lower coverages. Even at 5x, imputation may not be a particularly suitable method for generating accurate genome-wide diploid genotypes for downstream population genetic analysis if the target and reference population are too far diverged, and may introduce cryptic biases. Interestingly, applying genotype probability filters did not improve clustering to true samples significantly, while the clustering of lower-coverage samples was more noisy. These results suggest that the common practice of genotype probability filtering may not always yield substantial improvements for downstream diploid inference, and instead introduce greater noise.

In addition, PCA presents somewhat of a baseline for population genetic inferential methods, in that it does not rely on signals that correlate amongst sites (i.e., haplotypes). More sophisticated haplotype-based analysis such as IBD inference and chromosome painting will likely result in even larger biases, as our results suggest divergence between the reference and target population has a significant impact on haplotype switch-error, with both increases in coverage and reference size having very little effect once d is ≥ 1500.

Ethics Statement

Not applicable.

Consent for Publication

Not applicable.

Competing Interests

Krishna R. Veeramah is a member of the Editorial Board of the journal Human Population Genetics and Genomics. The author was not involved in the journal’s review of or decisions related to this manuscript. The author has declared that no other competing interests exist.

Availability of Data and Material

All simulation scripts are available at github.com/marianaer/imputation.

Author Contributions

MER and KRV conceived the study and wrote the paper, MER performed the analysis.