GBE HuConTest: Testing Human Contamination in Great Ape Samples Martin Kuhlwilm 1,*, Claudia Fontsere 1 , Sojung Han 1 , Marina Alvarez-Estape 1 , and Tomas Marques-Bonet 1,2,3,4 1 Institut de Biologia Evolutiva, CSIC-Universitat Pompeu Fabra, PRBB, Barcelona, Catalonia, Spain 2 CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Barcelona, Spain 3 Institucio Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Catalonia, Spain 4 Institut Catala de Paleontologia Miquel Crusafont, Universitat Auto noma de Barcelona, Edifici ICTA-ICP, Cerdanyola del Vallès, Barcelona, Spain *Corresponding author: E-mail: martin.kuhlwilm@upf.edu. Accepted: 19 May 2021 Abstract Modern human contamination is a common problem in ancient DNA studies. We provide evidence that this issue is also present in studies in great apes, which are our closest living relatives, for example in noninvasive samples. Here, we present a simple method to detect human contamination in short-read sequencing data from different species: HuConTest. We demonstrate its feasibility using blood and tissue samples from these species. This test is particularly useful for more complex samples (such as museum and noninvasive samples) which have smaller amounts of endogenous DNA, as we show here. Key words: contamination, nonhuman primates, next-generation sequencing, fecal DNA, ancient DNA. Significance Human contamination can be a confounding factor in genomic studies, especially in the case of fecal, museum, or ancient DNA from great apes. It is important for quality assessment, screening purposes, and prioritization to identify and quantify such contamination. The tool presented here is a simple and versatile method for this purpose and can be applied to a wide range of sample types. Main Text devise a strategy based on diagnostic sites dispersed across Contaminationfrom exogenous sources is a problem common the autosomes which can help detecting human contamina- in ancient DNA, where multiple tools exist (reviewed by tion in an unbiased manner and with sparse data available. Peyregne and Prüfer (2020)), as well as in studies of nonhu- man primates (Prado-Martinez et al. 2013). Specifically, hu- Determination of Diagnostic Sites man contamination may occur in great ape samples of various origin and quality. Previously, differences in the mitochondrial We used previously published diversity data on high-coverage genome between species were used to assess contamination genomes from all great apes and modern humans (supple- (Prado-Martinez et al. 2013), which is a sensible strategy for mentary table S1, Supplementary Material online, fig. 1A), high-coverage data. However, this approach is of limited use specifically, genomes from 58 chimpanzees and 10 bonobos for shallow shotgun sequencing, especially of samples with (Pan clade) (Prado-Martinez et al. 2013; De Manuel et al. low endogenous DNA content, such as fecal, historical, or 2016), 43 gorillas (Gorilla clade) (Prado-Martinez et al. ancient samples, as well as sequencing data obtained after 2013; Xue et al. 2015), 27 orangutans (Pongo clade) enrichment through capture (Fontsere et al. 2020). Here, we (Prado-Martinez et al. 2013; Nater et al. 2017), and 19 ß The Author(s) 2021. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Genome Biol. Evol. 13(6). doi:10.1093/gbe/evab117 Advance Access publication 26 May 2021 1 Kuhlwilm et al. GBE modern humans from Africa and western Eurasia (Mallick et allele in the test species. We then calculate the number of al. 2016). The latter was chosen to represent the most com- reads matching the target species complex allele and subtract mon variation across modern humans. All genomes were this value from the total read depth, thus retrieving the num- processed as described previously (De Manuel et al. 2016): ber of reads matching the human reference allele (more Sequencing data were mapped to the human genome (hg19) strictly, not matching the target species complex allele). We using BWA-MEM 0.7.7 (Li and Durbin 2009), PCR duplicates perform this calculation for each chromosome separately in were removed using samtools (Li et al. 2009), and reads were order to obtain the standard variation and report the genome- locally realigned around indels using the GATK IndelRealigner wide point estimate, one standard deviation, and the number 3.4-46 (McKenna et al. 2010). Genotypes were obtained in- of positions observed by the test. The latter value is useful to dividually using GATK UnifiedGenotyper with the assess the reliability of the test at extremely shallow sequenc- EMIT_ALL_SITES parameter, and GVCFs from individuals ing. The test can be applied to files with a bam or cram ex- were merged with GATK CombineVariants. The three species tension, containing short sequencing reads mapped to the complexes Pan, Gorilla, and Pongo were then filtered sepa- human genome (hg19 or GRCh38). The basic filtering at rately: biallelic SNPs within each species complex together this step can be simple, but it is advisable to remove adapter with humans were retrieved and filtered to exclude repetitive sequences (Schubert et al. 2016) and PCR duplicates to assess regions of the genome and regions with low mappability the unique contaminant fraction, as well as unmapped reads, (35mer mappability). Finally, for each individual, genotypes nonprimary alignments, and sequences with a low mapping were set to missing at sequencing coverage lower than 6 quality (<30). We specifically recommend filtering the and higher than 100, and with a mapping quality lower sequences on fragment/insert length (35 base pairs) to avoid than 20. spurious alignments, which may happen at a high rate in the We retrieved SNPs where at least 98% of the chromo- case of samples with large amounts of bacterial DNA (Meyer somes in the species complex showed the alternative allele et al. 2016). (different from the human reference), with <5% of missing We tested the contamination test by artificially introducing genotypes, and where all modern human chromosomes in- modern human sequencing reads into bam files from the cluded in this study carried the human reference allele, with- other species (in silico contamination), using eight human out allowing for missing genotypes. We allowed for residual individuals that were not part of the reference panel (supple- amounts of human-reference-like genotypes in the great ape mentary table S2, Supplementary Material online) (Auton species, to account for residual sources of error in the refer- et al. 2015), and great ape samples from other studies ence set that might result in erroneous rare variation, and we (Locke et al. 2011; Prüfer et al. 2012; Besenbacher et al. deemed tolerating these a conservative strategy for determin- 2019). First, each human bam file was downsampled to ing diagnostic sites. Across the whole genome, we find 1.14 M reads and merged with a chimpanzee bam file 4,460,987 diagnostic sites for Pan species, 6,981,108 diag- (ERR032960), to simulate 5% of human contamination. nostic sites for Gorilla species, and 7,518,570 diagnostic sites Since the read length differs between sequencing libraries for Pongo species. The differences between species are par- from different studies, we account for the expected amount tially explained by their evolutionary divergence to humans of human contamination by using the percentage of human (Kaessmann and P€a€abo 2002; Prado-Martinez et al. 2013; base pairs added to the final bam file. After running the hu- Kuhlwilm et al. 2016) but also the number of individuals man contamination test in each file, we detect an average of used, as well as sequencing quality and coverage. We used 5.5% human contamination (supplementary table S2, the R package rtracklayer (Lawrence et al. 2009) to perform a Supplementary Material online), with minimal differences be- liftover of these positions to the human genome version 38 tween humans from different world regions. When testing a (GRCh38). gradient of increasing amounts of introduced human sequen- ces from 0.1% to 39% to a chimpanzee bam file (sup- plementary table S3, Supplementary Material online, fig. 1C), Contamination Assessment and the contamination is estimated correctly. The test is perform- Performance ing well for in silico contamination from modern humans in Contamination is defined here as the proportion of observed each of the great ape species (supplementary table S4, human reference alleles at diagnostic positions in short se- Supplementary Material online). quencing reads (fig. 1B). The test itself is wrapped in an R We also determined the inferred amount in the case of script (confirmed to work with R versions 3.2.0, 3.4.4, 3.5.0, crosstesting, that is, performing the test of diagnostic sites 3.6.0, and 4.0.1 (R Core Team 2015)), to directly process the for one clade in samples from the other clades (supplementary number of reads carrying the reference or the alternative al- table S5, Supplementary Material online). Practically, this lele. We use samtools mpileup (tested for version 1.0 and 1.9) could be the relevant when, for example, a fecal sample to obtain read depth and alternative read depth at diagnostic was attributed to a chimpanzee but originated from a gorilla. sites, and join these data with information on the alternative We find contamination estimates of 44–82%, depending on 2 Genome Biol. Evol. 13(6): doi:10.1093/gbe/evab117 Advance Access publication 26 May 2021 HuConTest GBE A B C D E FIG. 1.—Summary of this study. (A) Schematic tree of the great ape species, with approximate divergence times (Besenbacher et al. 2019). (B) Schematic representation of the method. (C) Performance of detection of in silico-contamination in a gradient from 0.1 to 39%, point estimate 6 one standard deviation. (D) Performance when downsampling sequencing data from fecal samples with 1–3.5% of human contamination. Point estimates 6 one standard deviation. (E) Contamination estimates of blood samples for sequencing libraries from all species (violet; bonobo N ¼ 52, chimpanzee N ¼ 15; gorilla/orangutan N ¼ 47), fecal samples before and after capture (brown; N ¼ 109, without sample N42003 which has high levels of non-great ape contamination), museum samples (grey; N ¼ 8) and RNA-sequencing data (green; N ¼ 4). Genome Biol. Evol. 13(6): doi:10.1093/gbe/evab117 Advance Access publication 26 May 2021 3 Kuhlwilm et al. GBE the species combination. This is a consequence of the shared competitive mapping of sequences when large amounts of ancestry between humans and the other clades, where, for contamination are detected. This will help to determine the example, some of the chimpanzee-diagnostic alleles are the species of origin, for example using BBSplit (https://source- same in humans and gorillas, and other the same in chimpan- forge.net/projects/bbmap/, last accessed May 27, 2021) with zees and gorillas. This demonstrates that the test is clade- a reference panel of great apes, and possibly other primate specific, and large amounts of reads that do not carry species living in the same habitat. We applied this method to clade-specific alleles will be detected when a different primate these four samples (N42003_Shotgun1, N31908_Shotgun1, species are present. N33104_Shotgun1, and N41207_Shotgun1), and find that the main contaminant in one sample is most likely another primate rather than human (supplementary table S7, Application to Other Sample Types Supplementary Material online). It is known that chimpanzees We first applied the test to blood samples from all great ape hunt other primates (Boesch and Boesch 1989), and DNA species, which are generally expected to contain at most small from primate prey can persist in the feces. We conclude amounts of human contamination. For 67 randomly chosen that the design of the contamination test presented here is sequencing libraries from seven chimpanzees and four bo- able to identify reads carrying mutations that differ from the nobo individuals (Prüfer et al. 2012), we found an average target species, even if these are not human-specific. When of 0.28% (0.13–0.61%) of reads that are putatively due to applying BBsplit method to in silico-contaminated samples, human contamination (fig. 1E). Four tissue samples from we confirm humans as the source of the contamination— chimpanzees (White et al. 2019) show low estimates of con- although with less precision regarding the amount of contam- tamination (0.03–0.067%), as expected for samples likely not ination when compared to our method—while the majority containing true contamination. Similar results are obtained for of unambiguously mapped sequences align to the target spe- four libraries from gorilla (0.033–0.159%, on average cies (supplementary table S7, Supplementary Material online). 0.075%) and 43 libraries from orangutan (0.08–0.35%, on Our analysis shows that DNA extracts/libraries from fecal average 0.22%) blood samples (Locke et al. 2011; samples are occasionally contaminated, and may need to be Besenbacher et al. 2019). We conclude that traces of putative removed from certain downstream analyses. Hence, it is ad- human contamination are observed, if at all, only at very small visable to perform a contamination test for sequencing data amounts in sequencing data from great ape blood samples. from this type of sample, comparable to ancient and historical These estimates are conservative, since sequencing errors, samples. We assessed the power to detect human contami- mapping reference bias, and variation in these individuals nation with very shallow sequencing, by downsampling the may contribute to these numbers, especially considering sequencing reads of the three fecal samples from White et al. that error rates of these sequencing technologies were de- (N31908_Shotgun1, N33104_Shotgun1, and creasing after the publication of some of these studies (Locke N41207_Shotgun1) with 1–3.5% human contamination. et al. 2011; Prüfer et al. 2012). We also note that results for We downsampled these in several steps down to 1,000 data mapped to hg19 and hg38 are almost identical (supple- production reads (supplementary table S8, Supplementary mentary table S6, Supplementary Material online). Material online), and calculated the estimated amount of hu- We then applied the contamination test to noninvasive man contamination. These results (fig. 1D) confirm that our samples which usually contain small amounts of host DNA method is robust in confidently detecting human contamina- and may require target hybridization methods to obtain suf- tion even in the case of very shallow sequencing, as low as ficient data (Hernandez-Rodriguez et al. 2017; Fontsere et al. 1,000 reads aligned to the human reference genome, al- 2020). We applied our method to shotgun and exome cap- though with a high standard deviation. In the case of fecal ture sequencing data that were obtained from the same 109 samples with around 5% of estimated host DNA, this could sequencing libraries from chimpanzee fecal samples (White et be as little as 20,000 production reads, making the test al. 2019). We found an average of 0.35% (0–24.6%) human applicable to shallow data from an initial screening procedure contamination for the precapture (shotgun) and 0.32% (Fontsere et al. 2020). (0.05–21%) human contamination in the postcapture We also applied the test to published sequencing data (enriched) sequencing data (supplementary table S6, from eight museum bone samples from gorillas (van der Supplementary Material online, fig. 1E), with strong correla- Valk et al. 2019). Here, we find an estimated contamination tion for the same samples (r ¼ 0.99, P value < 2.2  1016). of on average 0.68% (0.55–0.72%), which is slightly lower We find one sample with an estimate of 24.6% and three than the reported estimates which were based on mitochon- more samples with more than 1% of human contamination drial diagnostic loci (0.28–1.67%, on average 1%), and (supplementary table S6, Supplementary Material online). In slightly higher than estimates for blood samples, as expected the case of fecal samples collected from the field that may for museum specimens that have been handled by humans. contain other mammalian DNA than the target species Contamination estimates from mitochondrial and nuclear loci through diet or misidentification, it is advisable to perform a from the same sample have been found to not be identical in 4 Genome Biol. Evol. 13(6): doi:10.1093/gbe/evab117 Advance Access publication 26 May 2021 HuConTest GBE hominin samples (Prüfer et al. 2014), and at shallow sequenc- Besenbacher S, Hvilsom C, Marques-Bonet T, Mailund T, Schierup MH. ing coverage, a small number of reads would overlap with 2019. Direct estimation of mutations in great apes reconciles phylo- genetic dating. Nat Ecol Evol. 32(2): 86–292. mitochondrial diagnostic loci. Still, the differences between Boesch C, Boesch H. 1989. Hunting behavior of wild chimpanzees in the these methods are minor, and results on data mapped to Taı̈ National Park. Am J Phys Anthropol. 78(4): 547–573. hg19 and hg38 are almost identical (supplementary table Brawand D, et al. 2011. The evolution of gene expression levels in mam- S6, Supplementary Material online), as is the case for blood malian organs. Nature 478(7369):343–348. samples. Finally, we performed the contamination test on Fontsere C, et al. 2020. Maximizing the acquisition of unique reads in non- invasive capture sequencing experiments. Mol Ecol Resour (21):745– RNA-sequencing data from great ape tissue samples 761. doi: 10.1111/1755-0998.13300. (Brawand et al. 2011), mapped using tophat2 (Kim et al. Hernandez-Rodriguez J, et al. 2017. The impact of endogenous content, 2013). We find slightly higher amounts of contamination replicates and pooling on genome capture from faecal samples. Mol (supplementary table S6, Supplementary Material online), ei- Ecol Resour (18):319–333. doi: 10.1111/1755-0998.12728. ther due to real contamination in the samples or to higher Kaessmann H, P€ a€ abo S. 2002. The genetical history of humans and the great apes. J Intern Med. 251(1): 1–18. error rates and mapping bias in transcriptome data compared Kim D, et al. 2013. TopHat2: accurate alignment of transcriptomes in the to genome sequencing data. presence of insertions, deletions and gene fusions. Genome Biol. 14(4): R36. Kuhlwilm M, et al. 2016. Evolution and demography of the great apes. Supplementary Material Curr Opin Genet Dev (41):124–129. doi: 10.1016/j.gde.2016.09.005. Supplementary data are available at Genome Biology and Lawrence M, Gentleman R, Carey V. 2009. rtracklayer: an R package for interfacing with genome browsers. Bioinformatics 25(14): 1841–1842. Evolution online. Li H, 1000 Genome Project Data Processing Subgroup, et al. 2009. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25(16): 2078–2079. Method Availability Li H, Durbin R. 2009. Fast and accurate short read alignment with The contamination test script including documentation is pub- Burrows-Wheeler transform. Bioinformatics 25(14): 1754–1760. licly available on GitHub: https://github.com/kuhlwilm/ Locke DP, et al. 2011. Comparative and demographic analysis of orang- utan genomes. Nature 469(7331): 529–533. HuConTest (last accessed May 27, 2021). Files with the diag- Mallick S, et al. 2016. The Simons genome diversity project: 300 genomes nostic positions are publicly available on FigShare (doi: from 142 diverse populations. Nature 538(7624): 201–206. 10.6084/m9.figshare.14237834). De Manuel M, et al. 2016. Chimpanzee genomic diversity reveals ancient admixture with bonobos. Science (80-) 354(6311): 477–481. McKenna A, et al. 2010. The Genome Analysis Toolkit: MapReduce frame- Data Availability work for analyzing next-generation DNA sequencing data. Genome Res. 20(9): 1297–1303. There are no new data associated with this article. Meyer M, et al. 2016. Nuclear DNA sequences from the Middle Pleistocene Sima de los Huesos hominins. Nature 531(7595): 504–515. Acknowledgments Nater A, et al. 2017. Morphometric, behavioral, and genomic evidence for a new orangutan species. Curr Biol. 27(22): 3487–3498.e10. Wethank Marc de Manuel for help with data preparation. Peyregne S, Prüfer K. 2020. Present-day DNA contamination in ancient M.K. is supported by “la Caixa” Foundation (ID DNA datasets. BioEssays 42(9): 2000081. 100010434), fellowship code LCF/BQ/PR19/11700002. Prado-Martinez J, et al. 2013. Great ape genetic diversity and population history. Nature 499(7459): 471–475. M.A.E. is supported by an FPI (Formacio n de Personal Prüfer K, et al. 2012. The bonobo genome compared with the chimpanzee Investigador) PRE2018-083966 from Ministerio de Ciencia, and human genomes. Nature 486(7404): 527–531. n. T.M.-B is supported by funding Universidades e Investigacio Prüfer K, et al. 2014. The complete genome sequence of a Neanderthal from the European Research Council (ERC) under the from the Altai Mountains. Nature 505(7481): 43–49. European Union’s Horizon 2020 research and innovation pro- R Core Team. 2015. R: a language and environment for statistical com- puting. http://www.r-project.org/. Accessed May 27, 2021. gramme (grant agreement No. 864203), BFU2017-86471-P Schubert M, Lindgreen S, Orlando L. 2016. AdapterRemoval v2: rapid adapter (MINECO/FEDER, UE), “Unidad de Excelencia Marıa de trimming, identification, and read merging. BMC Res Notes. 9:88. Maeztu”, funded by the AEI (CEX2018-000792-M), van der Valk T, Dıez-del-Molino D, Marques-Bonet T, Guschanski K, Dal en Howard Hughes International Early Career, Obra Social “La L. 2019. Historical genomes reveal the genomic consequences of re- Caixa” and Secretaria d’Universitats i Recerca and CERCA cent population decline in eastern gorillas. Curr Biol. 29(1): 165–170.e6. Programme del Departament d’Economia i Coneixement de White LC, et al. 2019. A roadmap for high-throughput sequencing studies la Generalitat de Catalunya (GRC 2017 SGR 880). of wild animal populations using noninvasive samples and hybridiza- tion capture. Mol Ecol Resour. 19(3): 609–622. Xue Y, et al. 2015. Mountain gorilla genomes reveal the impact of long- References term population decline and inbreeding. Science 348(6231): 242–245. Auton A, 1000 Genomes Project Consortium, et al. 2015. A global refer- ence for human genetic variation. Nature 526(7571): 68–74. Associate editor: David Enard Genome Biol. Evol. 13(6): doi:10.1093/gbe/evab117 Advance Access publication 26 May 2021 5