Ancestry-agnostic estimation of DNA sample contamination from sequence reads

New model-based methods accurately estimate genetic ancestry

In the absence of contamination, widely used methods such as LASER and TRACE are known to estimate genetic ancestry accurately. Because we propose using a new model-based approach to estimate the genetic ancestry (jointly with contamination rates), we first compared the accuracy of our new method, in the absence of contamination, with LASER and TRACE. We randomly chose 500 ethnically diverse samples from the 1000 Genomes Project low-coverage (4×) genomes and 500 African-American samples from the deeply sequenced (32×) genomes from the InPSYght project. We estimated their genetic ancestries using 100,000 SNPs from the HGDP reference panel (see Methods for details) and compared their genetic ancestry estimates obtained by LASER (using the same sequence data) and TRACE (using the hard-call genotypes). As illustrated in Figure 2, A, C, and E, the estimated PC coordinates of the 1000 Genomes individuals are located close to their corresponding HGDP populations across all three methods. Compared to TRACE and LASER, we observed that the estimated genetic coordinates from verifyBamID2 were the closest to the centroid of the corresponding HGDP population (Table 1) in four of the five populations (all except TSI). These results suggest that our method provides estimates at least as precise compared to those for other state-of-the-art methods.

Genetic ancestry estimates may be confounded by DNA contamination

Next, we constructed in silico–contaminated sequenced data from the 1000 Genomes Project and estimated contamination parameters and genetic ancestries jointly. We observed that, when sequences are contaminated between different continental populations, the genetic ancestry estimates in PC coordinates drift toward the contaminating population when contamination is ignored (Fig. 3A) or when assuming that intended and contaminating samples originated from the same population (Fig. 3C). As the contamination rate increases, drift increases (Fig. 3A,C,E).

However, when we accounted for possible differences in genetic ancestries between the two intended and contaminating samples using our new methods, PC coordinates remained similar to those for uncontaminated samples (Fig. 3E) and contaminated samples constructed from individuals that belong to the same population (Fig. 3B,D,F).

Robust, accurate, ancestry-agnostic estimation of DNA contamination

Next, we evaluated the effect of genetic ancestry misspecification in estimating DNA contamination rates. We constructed contaminated samples between various combinations of populations and compared the accuracy of estimated contamination rates using both the original methods which assume known allele frequencies and the new methods which estimate contamination rate and genetic ancestry jointly.

When contamination happens within the same population, running original methods with correct continental population allele frequencies specified provided accurate contamination estimates (Fig. 4A,E,I). However, using pooled allele frequencies, which would be a default option when it is infeasible to specify population information a priori before sequencing, consistently underestimated contamination rates. Bias was particularly large when intended individuals were of African ancestry.

Specifying incorrect population allele frequencies results in even larger contamination estimation bias. For example, using African allele frequencies on East Asian samples resulted in an average estimate of 2.9% contamination for samples with contamination 10% (Supplemental Table S1), implying that a large fraction of 10% contaminated samples within East Asian ancestry would not have been flagged for contamination-based exclusion at the contamination-exclusion threshold of 1%–3% used by many studies—for example, the Trans-Omics Precision Medicine (TOPMed) study (Natarajan et al. 2018).

Our results consistently demonstrated that the ancestry-agnostic method provides as accurate estimates as the original methods specified with correct population labels (Fig. 4A,E,I; Supplemental Table S1), and the estimates are substantially better than those from pooled allele frequencies or incorrectly specified allele frequencies (Table 2).

When the intended and contaminating populations are different, we observed that contamination is sometimes overestimated due to an increased fraction of heterozygous genotypes than expected by a given contamination rate under the single population model. Our method based on an unequal-ancestry model outperforms all the other methods in terms of overall bias and mean squared error (MSE) (Fig. 4; Supplemental Table S4), correcting for both upward and downward biases in various ancestry combinations. For example, the relative deviation of estimated to intended contamination rate (i.e., $\left|\hat{\alpha }/\alpha -1\right|$) is reduced by 80% (73%–88%) compared to the original verifyBamID with various population allele frequencies, suggesting reduced bias. MSE is also reduced by 92% (86%–97%). This robustness reflects the ability to incorporate differences in population allele frequencies between intended and contaminating individuals (Fig. 4B–D,F–H; Supplemental Table S1).

We also examined the accuracy of our methods for admixed populations by performing a similar experiment using the Mexican population (MXL) and obtained consistent results (Supplemental Table S2).

Results with deep whole-genome sequence data from the InPSYght study

Next, we applied our methods to 500 African-American samples from the InPSYght study (see Methods). Consistent with the results from our in silico contamination studies, we observed that the average contamination rate was 1.1-fold higher with the newer method (0.36% for unequal-ancestry, 0.37% for equal-ancestry) compared to the original method with pooled allele frequency (0.33%) (Fig. 5). The number of samples with an estimated contamination rate >1% increased from 16 (original method with pooled allele frequency) to 21 (unequal-ancestry method) or 23 (equal-ancestry method), suggesting our new method more rigorously screens for contaminated samples.

All 500 deeply sequenced genomes in InPSYght study are reported to be African-Americans, and indeed the estimated PC coordinates for all 500 individuals under all three methods lie between European and African samples. Compared to other methods to estimate genetic ancestry, our estimates resulted in tighter clustering along the European-African segment than LASER and similarly tight clustering to TRACE (Fig. 2B,D,F). For example, the correlation coefficient between the PC1 and PC2 coordinates were 0.927 for LASER, 0.981 for TRACE, and 0.985 for verifyBamID2, corroborating that verifyBamID2 results in a more precise estimate of African ancestry along the European-African segment in PC coordinates.

Impact of number of markers on accuracy, computational cost, and memory requirements

As we have shown previously (Jun et al. 2012), there are trade-offs between computation cost and accuracy of contamination estimates. Using as many as 100,000 variants results in an accurately estimated intended contamination rate. For example, the MSE of relative deviation (i.e., $\left|\hat{\alpha }/\alpha -1\right|$) was 0.02, 0.01, and 0.01 when the intended contamination was 1%, 2%, and 5%, respectively. When we use 10,000 variants, the MSEs modestly increased to 0.11, 0.04, and 0.01, respectively. When we use only 1000 variants, MSEs further increased to 0.69, 0.25, and 0.11, suggesting that the estimates may not be precise for the low contamination rate when using only 1000 variants (Supplemental Table S3).

We also evaluated the computational cost and memory consumption of verifyBamID2 on whole-genome sequence data with various coverages. For the BAM files from the 1000 Genomes whole-genome sequence data (4.3–5.1× coverage), the average wall-clock running time was 5.5 min with a single thread and peak memory consumption was 505 MB when using 10,000 markers in a server with a Xeon 2.27 GHz processor. When using 100,000 markers, the average wall-clock running time was 20.5 min with a single thread and 8.0 min with four threads, and peak memory consumption was 528 MB.

For deep genome data from the InPSYght study (31× coverage) stored in CRAM format, the average wall-clock time was 17.3 min and peak memory consumption was 514 MB when using 10,000 markers. For 100,000 markers, the average wall-clock time was 155.6 min (single thread) or 96 min (four threads) and peak memory consumption was 548 MB.

Likelihood-based mixture model for DNA sequence contamination

In our previous contamination detection methods (Jun et al. 2012), we assumed that the DNA sequence reads from an intended sample are contaminated by sequence reads from, at most, one contaminating sample from the same population and that the population allele frequencies of all analyzed genetic variants are known. For each bi-allelic variant i (1 ≤ i ≤ m), let b_ij ∈ {R,A,O} (1 ≤ j ≤ D_i) be the observed base call representing the reference allele (R), alternate allele (A), or other allele (O) for the j-th read that overlaps the variant; D_i is the observed sequence depth at variant i. Let e_ij ∈ {0,1} be a random variable indicating whether a sequencing error did (1) or did not (0) occur for observed base b_ij; we assume e_ij follows a Bernoulli distribution with success probability ${10}^{-Q_{ij}/10}$ where Q_ij is a Phred-scale base quality score of b_ij. In the absence of contamination, if the true genotype $g_i^s\in \left\{0,1,2\right\}$ represents the count of alternate alleles of the sequenced sample s ∈ {1,2}, then $\mathrm{Pr}\left(b_{ij}|g_i^s,e_{ij}\right)$ can be easily represented as in Table 3, making the simplifying assumption of equally likely errors across four possible nucleotides.

We assume that the observed sequence reads are a (1 − α):α mixture of intended and contaminating reads given a contamination rate 0 ≤ α ≤ 1. Let $g_i^1$and $g_i^2$ represent the true genotypes of the intended and contaminating samples at variant i, respectively. Then, the mixture model likelihood of each observed base becomes

Assuming a homogenous population with known population allele frequency f_i and Hardy-Weinberg Equilibrium, $\mathrm{Pr}(g_i^2;f_i)$ follows a Binomial(2, f_i) distribution. Under the simplifying assumption of independent variants, the likelihood of the contamination rate becomes

The maximum likelihood estimate (MLE) of contamination rate $\hat{\alpha }$ can be obtained using Brent's algorithm (Brent 1973).

As we previously reported (Jun et al. 2012), this model assumes correctly specified population allele frequencies f_i.

Likelihood-based estimation of genetic ancestry (in the absence of contamination)

We extend this model to incorporate genetic ancestry. The key idea of this extension is to use the individual-specific allele frequency (Hao et al. 2015; Conomos et al. 2016) to model the likelihood of the sequence reads. Several methods, including spatial ancestry analysis (SPA) (Yang et al. 2012) and logistic factor analysis (LFA) (Hao et al. 2015), previously proposed modeling allele frequency as a function of genetic ancestry via principal component (PC) coordinates.

Let G be an m × n genotype matrix (where G_ir = 0, 1, or 2 is the number of nonreference alleles at variant i in individual r) of a genetically diverse reference panel of size n, such as 1000 Genomes or HGDP. We define ISAF f_i(0 ≤ f_i ≤ 1) for variant i as a weighted average of genotypes from the reference panel ($f_i=\sum _{r=1}^n{w}_r{G}_{ir}$), where 0 ≤ w_r ≤ 1 and G_ir ∈ {0,1,2} for individual r. For a homogeneous population, w_r = 1/2n results in a pooled allele frequency across all individuals in the reference panel. If each individual can be categorically represented as a one of k mutually exclusive subpopulations, the population-specific allele frequency for the subpopulation s ∈ {1,2, · · · ,k} can be represented as $w_r={\displaystyle \frac{I\left(s_r=s\right)}{2n_s}}$, where s_r ∈ {1,2, · · · ,k} represents the subpopulation that individual r belongs to, and n_s represents the size of subpopulation s. More generally, if individuals’ genetic ancestry is represented as continuous variables (such as PCs, SPAs, or LFAs), the individual-specific allele frequency can be represented as a function of the continuously represented genetic ancestry (Wang et al. 2014; Hao et al. 2015).

The estimated ISAF can be viewed as one-half times the genotype dosages approximated from a fixed number(=K) of factors, such as PCs, SPAs, or LFAs. In our method, we used a linear model to estimate ISAF from PCs, similar to previous studies (Hao et al. 2015; Conomos et al. 2016). Given the reference panel genotype matrix G, let $1/2\hat{G}$ be the ISAF matrix as a function of top K factors. ISAF matrix $1/2\hat{G}$ should well approximate 1/2G. For example, under a linear model, typical principal component analysis takes the singular value decomposition of the mean-centered genotype matrix $\overline{G}=G-2\mathit{\mu }{1}_{\mathit{n}}^{\mathit{T}}=UD{V}^T$, where μ = 1/2nG1_n is the pooled allele frequencies and 1_n is the column-vector of ones. Using the top K eigenvalues and corresponding eigenvectors U^(K),D^(K),V^(K) from the SVD, it is known that $\hat{G}=1/2U^{\left(K\right)}{D}^{\left(K\right)}{[V^{\left(K\right)}]}^T+\mathit{\mu }{1}_{\mathit{n}}^{\mathit{T}}$ minimizes $G-{\hat{G}}_2=\sum _{i,r}{(G_{ir}-{\hat{G}}_{ir})}^2$ among all possible rank K matrices (Pearson 1901), making it a good proxy for the ISAF matrix.

For a new individual s with genetic ancestry represented as x_s ∈ R^K in the PC (eigenvector) space of the reference panel, the ISAF for i-th variant can be modeled as $f_i({\mathit{x}}_{\mathit{s}})={\displaystyle \frac{1}{2}{\mathit{u}}_{\mathit{i}}^{\left(K\right)}{D}^{\left(K\right)}{\mathit{x}}_{\mathit{s}}^{\mathit{T}}+{\mu }_i}$, where ${\mathit{u}}_{\mathit{i}}^{\left(K\right)}$ is i-th row of U^(K) and μ_i is the i-th element of μ. To avoid a boundary condition, we constrain ɛ/2n ≤ f_i(x_s) ≤ 1 − ɛ/2n for a fixed ɛ (we used ɛ = 0.5 in our experiments). Then, the overall likelihood of an individual's genetic ancestry x is

where g_i represents the unobserved genotype of the sequenced sample at variant i. The maximum-likelihood genetic ancestry coordinates can be estimated as ${\hat{\mathit{x}}}_{\mathit{s}}=\mathrm{argma}{\mathrm{x}}_{{\mathit{x}}_{\mathit{s}}\in {\mathbb{R}}^k}L\left({\mathit{x}}_{\mathit{s}}\right)$ using the Nelder-Mead (Nelder and Mead 1965) algorithm, starting with PC coordinates of a randomly selected individual from the reference panel. In all our experiments, we always obtained consistent estimates of ${\hat{\mathit{x}}}_{\mathit{s}}$ regardless of start values with K = 4, which is the default parameter of our implementation. Using K = 4 gave us noticeably more precise estimates of contamination rates and genetic ancestry than smaller K (Supplemental Fig. S1). Using larger values of K (e.g., K = 8) substantially increased the computational time of the Nelder-Mead algorithm and failed to converge occasionally.

Joint estimation of genetic ancestry and DNA contamination

Because our goal is to obtain unbiased estimates of the DNA contamination rate α agonistic to prior knowledge of the genetic ancestry, we propose to jointly estimate α and ancestry by combining the models described in the previous sections. Let x₁, x₂ ∈ R^K be the genetic ancestries of the intended and contaminating samples. Then, the likelihood under the combined model is

When the contamination rate α ≈ 0, the parameters corresponding to x₂ do not contribute (much) to the likelihood, and the estimates of x₂ may be unstable. To address this problem, we initially assume that the intended and contaminating samples are from the same population x₁ = x₂ (‘equal-ancestry’ model) and then repeat the analysis allowing for x₁ ≠ x₂ (‘unequal-ancestry’ model). The dimension of parameter space for the unequal-ancestry model is 2k + 1. We choose final parameter estimates between the two models based on the Akaike Information Criterion (Akaike 1974).

Evaluation on in silico–contaminated data based on 1000 Genomes Project samples

We constructed in silico–contaminated DNA sequence reads using aligned low-coverage whole-genome sequence reads from the 1000 Genomes phase 3 project (The 1000 Genomes Project Consortium 2015). We filtered out unmapped and mark-duplicated reads and then randomly sampled aligned sequence reads proportional to the intended contamination rates α ∈ {0.01, 0.02, 0.05, 0.1, 0.2}. To match the mixing proportion of sequence reads originated from intended and contaminating to be (1 − α):α, each read was sampled with probability (1 − α) and B₁/B₂α from each sample, where B₁ and B₂ are number of aligned bases from unique reads from intended and contaminating samples. We selected four populations, CHS (Han Chinese South), GBR (British in England and Scotland), MXL (Mexican Ancestry from Los Angeles, CA, USA), YRI (Yoruba in Ibadan, Nigeria), and arbitrarily selected 10 pairs of individuals with similar sequencing depths within the same population and across populations. To estimate genetic ancestry and/or contamination rate for these in silico–contaminated sequence reads, we used a reference panel of 938 HGDP (Cavalli-Sforza 2005) individuals across 1000, 10,000, and 100,000 randomly chosen SNPs (pooled MAF > 0.5%), avoiding variants masked by the 1000 Genomes Project (The 1000 Genomes Project Consortium 2015). When we compared estimated genetic ancestry with LASER, we used the same set of selected SNPs and sequence reads as input. For TRACE, we used genotypes from the phase 3 release (for 1000 Genomes) or an interim call set from the GotCloud software tool (Jun et al. 2015) (for InPSYght, see the next section for details) on the same SNP set.

Experiment with real sequence data from the InPSYght study

Next, we applied our method to 500 deeply sequenced (mean depth 32×) genomes from the first two batches of the InPSYght study. For each sample, we evaluated the results from the six models: (1) the original verifyBamID using pooled allele frequencies; the original verifyBamID using (2) African, (3) East Asian, and (4) European allele frequencies; (5) the new verifyBamID2 under the equal-ancestry model; and (6) verifyBamID2 under the unequal-ancestry model. To calculate pooled, population-specific, and individual-specific allele frequencies, we used the 1000 Genomes phase 3 reference panel (n = 2504), randomly selecting 100,000 SNPs among the sites also polymorphic in Illumina Human Omni 2.5 array, with the same filtering criteria (MAF > 5% and 1000 Genomes mask) as above.