When sequencing an ancient DNA sample from a hominin fossil, DNA from present-day humans involved in excavation and extraction will be sequenced along with the endogenous material. This type of contamination is problematic for downstream analyses as it will introduce a bias towards the population of the contaminating individual (s). Quantifying the extent of contamination is a crucial step as it allows researchers to account for possible biases that may arise in downstream genetic analyses. Here, we present an MCMC algorithm to co-estimate the contamination rate, sequencing error rate and demographic parameters—including drift times and admixture rates—for an ancient nuclear genome obtained from human remains, when the putative contaminating DNA comes from present-day humans. We assume we have a large panel representing the putative contaminant population (e. g. European, East Asian or African). The method is implemented in a C++ program called ‘Demographic Inference with Contamination and Error’ (DICE). We applied it to simulations and genome data from ancient Neanderthals and modern humans. With reasonable levels of genome sequence coverage (>3X), we find we can recover accurate estimates of all these parameters, even when the contamination rate is as high as 50%.
When sequencing a human genome using ancient DNA (aDNA) recovered from fossils, a common practice is to assess the amount of present-day human contamination in a sequencing library [1,2, 3,4, 5,6]. Several methods exist to obtain a contamination estimate. First, one can look at ‘diagnostic positions’ in the mitochondrial genome at which a particular archaic population may be known to differ from all present-day humans. Then, one counts how many aDNA fragments support the present-day human base at those positions. This is the most popular technique and has been routinely deployed in the sequencing of Neanderthal genomes [7,1]. However, contamination levels of the mitochondrial genome may sometimes differ drastically from those of the nuclear genome [8,9]. A second technique involves assessing whether the sample was male or female using the number of fragments that map to the X and the Y chromosomes. After determining the biological sex, the proportion of reads that are non-concordant with the sex of the archaic individual are used to estimate contamination from individuals of the opposite sex (e. g. Y-chr reads in an archaic female genome are indicative of male contamination) [8,1, 10,4]. Another method uses a maximum-likelihood approach to estimate contamination, but is only applicable to single-copy chromosomes, like the X chromosome in individuals known a priori to be male [11,12]. Finally, one last technique involves using a maximum-likelihood approach to co-estimate the amount of contamination, sequencing error and heterozygosity in the entire autosomal nuclear genome [1,3], using an optimization algorithm such as L-BFGS-B [13]. Afterwards, if the aDNA library shows low levels of present-day human contamination (< ∼2%), demographic analyses are performed on the sequences while ignoring the contamination. If the library is highly contaminated, it is usually treated as unusable and discarded. Neither of these outcomes is optimal: contaminating fragments may affect downstream analyses, while discarding the library as a whole may waste precious genomic data that could provide important demographic insights. One way to address this problem was proposed by skoglund et al. [14], who developed a statistical framework to separate contaminant from endogenous DNA fragments by using the patterns of chemical deamination characteristic of ancient DNA. The method produces a score which reflects the odds that a particular fragment is endogenous or not, based on these chemical patterns. This approach is effective at isolating truly endogenous fragments from contaminant fragments, but at the cost of potentially discarding some fragments that may not have chemical damage and still be endogenous. This becomes more problematic the younger the ancient DNA sample is, because younger samples will tend to have a higher proportion of non-deaminated ancient DNA, and so the method will lead users to discard a larger fraction of endogenous material. Instead of (or in addition to) attempting to separate the two type of fragments before performing a demographic analysis, one could incorporate the uncertainty stemming from the contaminant fragments into a probabilistic inference framework. Such an approach has already been implemented in the analysis of a haploid mtDNA archaic genome [15]. However, mtDNA represents a single gene genealogy, and, so far, no equivalent method has been developed for the analysis of the nuclear genome, which contains the richest amount of population genetic information. Here, we present a method to co-estimate the contamination rate, per-base error rate and a simple demography for an autosomal nuclear genome of an ancient hominin. We assume we have a large panel representing the putative contaminant population, for example, European, Asian or African 1000 Genomes data [16]. The method uses a Bayesian framework to obtain posterior estimates of all parameters of interest, including population-size-scaled divergence times and admixture rates.
We will first describe the probabilistic structure of our inference framework. We begin by defining the following parameters: We are interested in computing the probability of the data given the contamination rate, the error rate, the derived allele frequencies from the putative contaminant population (w) and a set of demographic parameters (Ω). We will use only sites that are segregating in the contaminant panel and we will assume that we observe only ancestral or derived alleles at every site (i. e. we ignore triallelic sites). In some of the analyses below, we will also assume that we have additional data (O) from present-day populations that may be related to the population to which the sample belongs. The nature of the data in O will be explained below, and will vary in each of the different cases we describe. The parameters contained in Ω may simply be the population-scaled times separating the contaminant population and the sample from their common ancestral population. However, Ω may include additional parameters, such as the admixture rate—if any—between the contaminant and the sample population. The number of parameters we can include in Ω will depend on the nature of the data in O. For all models we will describe, the probability of the data can be defined as: P [ a, d | r C, ϵ, w, Ω, O ] = ∏ j = 1 K P [ a j, d j | r C, ϵ, w j, Ω, O ] (1) where P [ a j, d j | r C, ϵ, w j, Ω, O ] = ∑ i = 0 2 P [ a j, d j | i, r C, ϵ, w j ] P [ i | Ω, O ] (2) Here, i is the true (unknown) genotype of the ancient sample, and P[i |Ω, O] is the probability of genotype i given the demographic parameters and the data. We focus now on computation on the likelihood for one site j in the genome. In the following, we abuse notation and drop the subscript j. Given the true genotype of the ancient individual, the number of derived and ancestral fragments at a particular site follows a binomial distribution that depends on the genotype, the error rate and the rate of contamination [1,3]: P [ a, d | i, r C, ϵ, w ] = a + d d q i d (1 - q i) a (3) where q 2 = r C w (1 - ϵ) + (1 - w) ϵ + (1 - r C) (1 - ϵ) (4) q 1 = r C w (1 - ϵ) + (1 - w) ϵ + (1 - r C) (1 - ϵ) / 2 + ϵ / 2 (5) q 0 = r C w (1 - ϵ) + (1 - w) ϵ + (1 - r C) ϵ (6) In the sections below, we will turn to the more complicated part of the model, which is obtaining the probability P[i|Ω, O] for a genotype in the ancient sample, given particular demographic parameters and additional data available. We will do this in different ways, depending on the kind of data we have at hand. First, we will work with the case in which O = y, where y is a vector of frequencies yj from an “anchor” population that may be closely related to the population of the ancient DNA sample. An example of this scenario would be the sequencing of a Neanderthal sample that is suspected to have contamination from present-day humans, from which many genomes are available. For all analyses below, we restrict to sites where 0 < yj < 1. Note that it is entirely possible (but not required) that y = w, meaning that, aside from the ancient DNA sample, the only additional data we have are the frequencies of the derived allele in the putative contaminant population, which we can use as the anchor population too. However, it is also possible to use a contaminant panel that is different from the anchor population (Fig 1A). We will assume we have sequenced a large number of individuals from a panel of the contaminant population (for example, The 1000 Genomes Project panel) and that the panel is large enough such that the sampling variance is approximately 0. In other words, the frequency we observe in the contaminant panel will be assumed to be equal to the population frequency in the entire contaminant population. In this case, Ω = {τC, τA}, where τA and τC are defined as follows: We need to calculate the conditional probabilities P[i|Ω, O] = P[i|y, τC, τA] for all three possibilities for the genotype in the ancient individual: i = 0,1 or 2. To obtain these expressions, we rely on Wright-Fisher diffusion theory (reviewed in Ewens [17]), especially focusing on the two-population site-frequency spectrum (SFS) [18]. The full derivations can be found in the S1 Text, and lead to the following formulas: P [ i = 0 | y, τ C, τ A ] = 1 - y * e - τ C - 1 2 * y * e - τ A - τ C + y y - 1 2 e - τ A - 3 τ C (7) P [ i = 1 | y, τ C, τ A ] = y * e - τ A - τ C + y 1 - 2 y e - τ A - 3 τ C (8) P [ i = 2 | y, τ C, τ A ] = y * e - τ C - 1 2 * y * e - τ A - τ C + y y - 1 2 e - τ A - 3 τ C (9) We generated 10,000 neutral simulations using msms [19] for different choices of τC and τA (with θ = 20 in each simulation) to verify our analytic expressions were correct (Fig 2). The probability does not depend on θ, so the choice of this value is arbitrary. The above probabilities allows us to finally obtain P[i | yj, Ω, O]. Although the above method gives accurate results for a simple demographic scenario, it does not incorporate the possibility of admixture from the ancient sample to the contaminant population. This is important, as the signal of contamination may mimic the pattern of recent admixture. We will assume that, in addition to the ancient DNA sample, we also have the following data, which constitute O: We can then estimate the remaining drift parameters, the error and contamination rates and the admixture time (β) and rate (α) between the archaic population and modern population Y. The diffusion solution for this three-population scenario with admixture is very difficult to obtain analytically. Instead, we use a numerical approximation, implemented in the program ∂a∂i [20]. We incorporated the likelihood functions defined above into a Markov Chain Monte Carlo (MCMC) inference method, to obtain posterior probability distributions for the contamination rate, the sequencing error rate, the drift times and the admixture rate. Our program—which we called ‘DICE’—is coded in C++ and is freely available at: http: //grenaud. github. io/dice/. We assumed uniform prior distributions for all parameters, and the boundaries of these distributions can be modified by the user. For the starting chain at step 0, an initial set of parameters X0 = {rC0, ϵ0, Ω0} is sampled randomly from their prior distributions. At step k, a new set of values for step k + 1 is proposed by drawing values for each of the parameters from normal distributions. The mean of each of those distributions is the value for each parameter at state Xk and the standard deviation is the difference between the upper and lower boundary of the prior, divided by a constant that can be increased or decreased to achieve a desired rate of acceptance of new states [21]. By default, this constant is equal to 1,000 for all parameters. The new state is accepted with probability: P [ a c c e p t ] = m i n 1, P [ a, d | X k + 1 ] P [ a, d | X k ] (10) where P[a, d | Xk] is the likelihood defined in Eq 1. Unless otherwise stated below, we ran the MCMC chain for 100,000 steps in all analyses, with a burn-in period of 40,000 and sampling every 100 steps. The sampled values were then used to construct posterior distributions for each parameter. Fu et al. [5] showed that, when estimating contamination, ancient DNA data can be better fit by a two-error model than a single-error model. In that study, the authors co-estimate the two genome-wide error rates along with the proportion of the data that is affected by each rate. Therefore, we also included this error model as an option that the user can choose to incorporate when running our program. Furthermore, we developed an alternative error estimation method that allows the user to flag transition polymorphisms, which are more likely to have occurred due to cytosine deamination in ancient DNA. These sites are therefore likely to be subject to different error rates than those common in present-day sequencing data [22,23]. Our program can then estimate two error rates separately: one for transitions and one for transversions. Finally, we incorporated an option to include an ancestral state misidentification (ASM) parameter, which should serve to correct for mispolarization of alleles [24]. The standard input for DICE is a file containing counts of particular ancestral/derived base combinations and SNP frequencies (see README file online). As an additional feature, we also developed a module for the user to directly input a BAM file and a file containing population allele frequencies for the anchor and contaminant panels, rather than the standard input. The user can either choose to convert the BAM file to native DICE format using a program provided with the software package and then run the program, or run it directly on the BAM file. In the latter case, instead of calculating genome-wide error parameters, the program will calculate error parameters specific to each sequenced fragment, based on mapping qualities, base qualities and estimated deamination rates at each site (see S2 Text).
We have developed a new method to jointly infer demographic parameters, along with contamination and error rates, when analyzing an ancient DNA sample. The method can be deployed using a C++ program (DICE) that is easy to use and freely downloadable. We therefore expect it to be highly applicable in the field of paleogenomics, allowing researchers to derive useful information from previously unusable (highly contaminated) samples, including archaic humans like Neanderthals, as well as ancient modern humans. Applications to simulations show that the error and contamination parameters are estimated with high accuracy, and that demographic parameters can also be estimated accurately so long as enough information (e. g. a large panel of modern humans) is available. The drift time estimates reflect how much genetic drift has acted to differentiate the archaic and modern populations since the split from their common ancestral population, and can be converted to divergence times in generations if an accurate history of population size changes is also available (for example, via methods like PSMC, [35]). Although we cannot perform proper model testing, we found via extensive simulations that the posterior mode of an MCMC run was a robust heuristic statistic to help detect which panel was most likely to have contaminated the sample. We caution, however, that the fact that a particular panel yields a higher posterior mode than another is no guarantee that it is a better fit to the data for demographic scenarios that may be different from the ones we simulated. We also applied our method to empirical data, specifically to two Neanderthal genomes at high and low coverage, a present-day high-coverage Yoruba genome, and several ancient genome sequences of varying degrees of coverage, some obtained via shotgun-sequencing and some via SNP capture. For the high-coverage Yoruba genome, we infer no contamination, as would be expected from a modern-day sample, and drift times indicating the Yoruba sample indeed belongs to an African population. The contamination and sequencing error estimates we obtained for the Altai Neanderthal are roughly in accordance with previous estimates [4]. The drift times we obtain under the three-population model for the African population (τC + τAfr) are approximately 0. 411 + 0. 009 = 0. 42 drift units. The geometric mean of the history of population sizes from the PSMC results in Prüfer et al. [4] give roughly that Ne ≈ 21,818 since the African population size history started differing from that of Neanderthals, assuming a mutation rate of 1. 25 * 10−8 per bp per generation. If we assume a generation time of 29 years, and use our drift time in the equation relating divergence time in generations to drift time (t/ (2Ne) ≈τ), this gives an approximate human-Neanderthal population divergence time of 531,486 years. This number roughly agrees with the most recent estimates obtained via other methods [4]. Additionally, the Neanderthal-specific drift time is approximately 6. 5 times as large as the modern human drift time, which is expected as Neanderthals had much smaller population sizes than modern humans [36,4]. The admixture rate from archaic to modern humans that we estimate is 1. 72%, which is consistent with the rate estimate obtained via methods that do not jointly model contamination (1. 5–2. 1%) [4]. In the case of the Altai Neanderthal, we observe that the sample was probably contaminated by one or more individuals with European ancestry. When testing modern human and Neanderthal ancient genomes of lower coverage than the Altai Neanderthal, we obtain reasonable parameter estimates for samples of medium to high-coverage. However, we run into problems in estimation when the samples are of low coverage. For these reasons, and from our simulation results, we recommend that our method should be used on nuclear genomes with >3X coverage. The method may converge under certain conditions at coverages as low as 0. 5X (for example, in the case of the Mezmaiskaya genome under the two-population model when using AFR as the anchor and contaminant panel), but, in such cases, we caution the user to check convergence is achieved before drawing any conclusions from the estimates. For SNP capture data, we obtain reliable estimates for samples with a minimum coverage of 500,000 sites that are polymorphic in the anchor panel. The demographic models used in our approach are simple, involving no more than three populations and a single admixture event. This is partly due to limitations of known theory about the diffusion-based likelihood of an arbitrarily complex demography for the 2-D site-frequency spectrum—in the case of the two-population method—and to the inability of ∂a∂i [20] to handle more than 3 populations at a time. In recent years, several studies have made advances in the development of methods to compute the likelihood of an SFS for larger numbers of populations using coalescent theory [37,38,39], with multiple population size changes and admixture events. We hope that some of these techniques could be incorporated in future versions of our inference framework.
When extracting and sequencing ancient DNA from human remains, a recurrent problem is the presence of DNA from the paleontologists, archaeologists or geneticists that may have handled the fossil. If a DNA library is highly contaminated, this will introduce biases in downstream analyses, so it is important to determine the amount of extraneous DNA. Different methods exist for this purpose, but few are applicable to the nuclear genome, and none of them can extract reliable genomic information from highly contaminated samples. Thus, samples with high rates of contamination are usually discarded. Here, we present a method to jointly estimate contamination and error rates, along with demographic parameters, like drift times and admixture rates. Our method can serve to uncover important details about the evolutionary history of archaic and early modern humans from ancient DNA samples, even if those samples are highly contaminated.
sequencing techniques
demography
social sciences
anthropology
ancient dna
contaminants
neanderthals
genome sequencing
materials science
genome analysis
paleontology
dna
paleogenetics
molecular biology techniques
paleoanthropology
research and analysis methods
genomic libraries
genomics
hominids
materials by attribute
molecular biology
people and places
hominins
biochemistry
physical anthropology
nucleic acids
earth sciences
genetics
biology and life sciences
physical sciences
computational biology
dna sequencing
Joint Estimation of Contamination, Error and Demography for Nuclear DNA from Ancient Humans