Genotype imputation from low-coverage whole genome sequencing data in wild baboons

Abstract

Genotype imputation from low-coverage whole genome sequencing (lcWGS) enables cost-effective generation of genome-wide genotype data, but its application in wild, non-model populations remains underexplored. Here, we developed and validated a genotype imputation pipeline for lcWGS data from 532 wild baboons (Papio cynocephalus x P. anubis hybrids) in the Amboseli ecosystem, Kenya. Using SHAPEIT5 for reference panel phasing and GLIMPSE2 for imputation, with an external reference panel of 250 non-Amboseli baboons (182 yellow, 68 anubis; mean coverage ~37x), we achieved high baseline imputation accuracy ($R^2 = 0.948 \pm 0.009$) even at 1x coverage. We then developed a multi-parameter post-imputation filtering framework incorporating GEM mappability, minor allele frequency, Hardy-Weinberg equilibrium, genotype posterior probability, and fraction of missingness filters. This filtering improved concordance to $R^2 = 0.991 \pm 0.004$ at 1x while retaining 39.3% of imputed variants. Benchmarking against three alternative methods (Beagle 5.4, STITCH, QUILT2) on the same data confirmed GLIMPSE2 as the top performer (dosage $R^2 = 0.90$ vs. 0.43, 0.28, and 0.88, respectively, for common variants at 1x coverage). Mendelian error analysis across 138 pedigree trios yielded extremely low error rates (mean 0.00094 per variant). As a downstream application, pedigree reconstruction with sequoia recovered all 469 known parent-offspring links and discovered 87 new ones, predominantly paternal relationships. Our results establish the feasibility of reference-based imputation in a wild, admixed primate population and provide a generalizable quality-control framework for imputation in non-model systems.

Introduction

Genetic data are essential for reconstructing evolutionary history, interpreting phenotypic variation, and projecting future evolutionary trajectories. However, genome-wide data can be difficult to obtain for non-model systems, especially from natural populations. While small data sets of distributed unlinked markers (e.g., microsatellites or reduced representation sequencing) provide important information about kinship, demographic history, and natural selection, denser genome-wide genotypes are necessary for genome-wide association studies (GWAS), expression quantitative trait loci (eQTL) mapping, identity-by-descent (IBD) analysis, and runs of homozygosity (ROH) estimation.

Whole genome sequencing (WGS) data are now increasingly feasible to generate even for non-model species. WGS does not require prior investment in species-specific genetic markers and provides near-complete representation of the genome. However, generating high-coverage data for large numbers of samples remains costly. For example, sequencing a mammalian genome (~3 Gb) to 30x costs approximately $540 per sample, while sequencing to 1x costs only ~$18, albeit with high uncertainty and many missing variants.

To achieve both cost-efficiency and genotyping accuracy, researchers have increasingly combined low-coverage sequencing (LCS) with genotype imputation. Imputation leverages linkage disequilibrium between nearby variants to identify shared haplotypes between a reference panel and target data, using these to fill in missing genotypes. Early imputation methods (IMPUTE, Beagle, MaCH, Minimac) were developed primarily for human GWAS from SNP arrays. More recently, tools designed for sequencing reads have emerged, including both commercial (loimpute) and open-source methods (STITCH, GeneImp, GLIMPSE, QUILT). Among these, GLIMPSE2 has become the most commonly used tool, achieving high computational efficiency and accuracy, and is now embedded in community-curated pipelines (e.g., nf-core phaseimpute) and commercial sequencing workflows.

LCS-based imputation has begun featuring in studies of non-model systems, including corals, geladas, Darwin's finches, captive baboons, rhesus macaques, and hihi/stitchbirds. However, predicting imputation accuracy a priori for a non-model population is difficult, as success depends jointly on reference panel quality, target population haplotype complexity, and the extent of haplotype sharing between reference and target populations.

Here, we developed and optimized a pipeline for genotype imputation from lcWGS data from wild baboons in the Amboseli ecosystem. This population consists of hybrid yellow (Papio cynocephalus) and anubis (Papio anubis) baboons, with continuous monitoring since 1971 and rich longitudinal behavioral and demographic data. We leveraged all existing external genomic resources from non-Amboseli yellow and anubis baboon populations to construct a reference panel, developed a multi-parameter post-imputation quality control framework, benchmarked GLIMPSE2 against three alternative imputation methods, and demonstrated the utility of imputed genotypes for pedigree reconstruction.

Results

The default GLIMPSE pipeline produces highly accurate imputed genotypes

We performed imputation using SHAPEIT5 for reference panel phasing and GLIMPSE2 for genotype inference. The reference panel included 182 non-Amboseli yellow baboons and 68 non-Amboseli anubis baboons, all sequenced to high coverage (yellow: $37.4\times \pm 9.6$ s.d.; anubis: $36.6\times \pm 5.5$ s.d.).

We first assessed whether variants discoverable in the reference panel capture variants segregating in the Amboseli population by comparing biallelic SNPs on chromosome 20 (the smallest baboon autosome) between the reference panel and a validation panel of 36 high-coverage Amboseli sequences ($26.9\times \pm 8.1$ s.d.; range: 20.1x-44.4x). We discovered 1,466,388 SNPs on chromosome 20 in the reference panel and 661,966 in the validation panel, with 580,890 shared. Importantly, 90.3% of variants with MAF > 0.05 in the Amboseli data were shared with the reference panel, though only 81.1% of variants with MAF $\leq$ 0.05 were shared.

To evaluate imputation accuracy, we down-sampled the 36 high-coverage samples to 10x, 5x, 3x, 2x, and 1x, performed imputation using default GLIMPSE2 parameters, and estimated concordance using the Pearson correlation $R^2$ between imputed and ground-truth genotype dosages. Even at 1x coverage, we achieved a mean within-sample concordance of $0.948 \pm 0.009$ s.d. (range: 0.927-0.961). Common variants with MAF > 0.05 in the reference panel achieved > 0.85 concordance at 1x, > 0.875 at 2x, and > 0.9 at 3x.

Post-imputation quality control filters substantially improve accuracy

We developed a multi-parameter filtering framework targeting known factors influencing imputation accuracy:

1. GEM mappability (GEM = 1): Exclude variants in regions with ambiguous read mapping
2. Minor allele frequency (MAF$_{\text{REF}}$ $\geq$ 0.05): Remove rare variants that are difficult to impute accurately
3. Hardy-Weinberg equilibrium (HWE $p \geq 0.01$): Filter variants deviating from equilibrium expectations
4. Genotype posterior probability (GP $\geq$ 0.9999): Remove individual genotype calls with low GLIMPSE2 confidence
5. Fraction of missingness (FM $\leq$ 0.2): Exclude variants with excessive missing genotypes after GP filtering

Filters were applied serially, with each parameter optimized to balance accuracy improvement against variant retention at 1x coverage. The GP and FM filters were applied last because applying GP filtering earlier preferentially removes heterozygotes, biasing HWE estimation.

The optimized filtering combination retained 248,623 biallelic SNPs (39.3% of imputed sites on chromosome 20) while substantially improving concordance. Mean within-sample concordance at 1x coverage increased from 0.948 to $0.991 \pm 0.004$ s.d. (range: 0.980-0.998). Variant-level concordance increased from > 0.85 to > 0.96 for common variants across all MAF bins.

Benchmarking against alternative imputation methods

To contextualize the performance of GLIMPSE2, we benchmarked it against three alternative imputation methods on the same data (chromosome 20, 36 samples at 1x coverage):

GLIMPSE2 achieved the highest dosage $R^2$ across all MAF bins, followed by QUILT2. Beagle 5.4 and STITCH performed substantially worse, likely because Beagle was originally designed for array-based genotyping data and STITCH reconstructs haplotypes de novo from the target population without using the external reference panel, which is disadvantageous when the target sample size is modest (N=36) relative to the haplotype diversity.

Near-elimination of Mendelian errors in a large pedigree

Using the optimized pipeline, we produced imputed genotypes across all autosomes for 532 Amboseli individuals ($5.8\times \pm 7.2$ s.d.; range: 1.2x-44.4x). We checked consistency with Mendelian inheritance for 259,074 imputed variants on chromosome 20 across 138 full parent-offspring trios (238 unique individuals).

The Mendelian error rate was extremely low (mean $= 0.00094 \pm 0.00055$ s.d.; range: 0.00017-0.00474), translating to approximately 1 violation per 1,000 variants (~243 errors in 259,074 variants). Trios with lower mean coverage had slightly higher error rates (Pearson's $R = -0.33$, $p = 6.47 \times 10^{-5}$). At the site level, we observed a mean of $0.98 \pm 0.44$ errors per site across 138 trios.

Application: pedigree reconstruction recovers known links and discovers new ones

As a downstream application, we reconstructed pedigree links for 532 individuals using imputed genotypes at 2,505 LD-pruned, highly variable sites (MAF > 0.35, $r^2 < 0.05$) with the R package sequoia. Comparison to the independently established Amboseli pedigree showed that the imputed pedigree fully recovered all known parent-offspring links (174 father-offspring and 295 mother-offspring links).

Notably, the analysis also discovered 87 new parent-offspring links, predominantly father-offspring relationships (n = 83). Of these, 71 involved parents not in the resequencing data set who were identified through clustering of half- or full-sibling offspring. For example, one target individual (DIW) was identified as a half-sibling of six other genotyped individuals who had previously been assigned the same father (AYU), despite AYU not being in the resequencing data set. Because these half-siblings have known different mothers, the paternal half-sibling relationship implies that AYU must also be DIW's father. This demonstrates the potential for using imputed genotypes to investigate kinship and pedigree structure beyond traditional parent-offspring assignment.

Discussion

The combination of low-coverage resequencing and genotype imputation is a cost-effective method for producing high-quality genome-wide genotype data. Here, we demonstrate that genotype imputation is highly accurate for lcWGS data from wild baboons, validating its feasibility in an otherwise well-characterized long-term field study of wild primates.

Reference panel construction for non-model systems. One key limitation of imputation in non-model systems is that reference panels are often small and sometimes drawn from the target population itself. Previous studies have shown that reference panels can be augmented using data from non-target populations of the same species. Here, we show that imputation can also succeed in naturally occurring hybrids, using a reference panel based on individuals from both parental ancestries drawn from wild and captive non-target populations. For common variants, our modest reference panel (N=250) captured most genetic variation in the target population and imputed genotypes with high accuracy.

Multi-parameter quality control framework. Consistent with previous studies, our results demonstrate that imputation accuracy is influenced by several site-level characteristics. Variants in complex genomic regions, rare variants, sites deviating from HWE, and genotypes with low posterior probabilities all contribute to imputation errors. Our serial filtering framework addresses each of these factors systematically. However, we caution that filtering disproportionately removes rare variants and distorts the allele frequency spectrum. This has minimal impact for GWAS or eQTL mapping but can affect selection tests relying on the allele frequency spectrum (e.g., Tajima's D, Fay and Wu's H).

Method benchmarking. Our comparison of four imputation methods confirms that reference-panel-based methods (GLIMPSE2, QUILT2) substantially outperform methods designed for array data (Beagle 5.4) or de novo haplotype reconstruction (STITCH) in this context. GLIMPSE2 achieved the highest accuracy across all MAF bins, supporting its use as the default imputation tool for non-model systems with available external reference panels.

Limitations. Our Mendelian error analysis is inherently limited to pedigree-connected individuals, potentially slightly overestimating accuracy relative to unrelated individuals. Additionally, all samples in this study were generated from high-quality blood or tissue-derived DNA with relatively even read distributions. Imputation from non-invasive samples (e.g., fecal DNA) or ancient DNA may produce poorer results due to highly uneven read distributions, even at equivalent genome-wide coverage.

Methods

Study subjects

Subjects were 532 individually identified baboons (266F, 266M) from a long-term study of hybrid yellow and anubis baboons in Amboseli, Kenya. Blood or tissue samples were collected 1989-2023 through darting or opportunistic tissue collection from corpses.

Data generation

The resequencing dataset includes data from 442 previously analyzed individuals plus 90 newly sequenced individuals. DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen) and libraries were prepared with the NEBNext Ultra II FS DNA Library Prep Kit. Sequencing was performed on Illumina NovaSeq 6000 or NovaSeq X platforms. Reads were processed following GATK Best Practices: adapter marking, mapping to the Panubis1.0 genome (GCF_008728515.1) with bwa-mem, duplicate removal, and quality filtering (MAPQ > 10). The final dataset had mean coverage $5.8\times \pm 7.2$ s.d. (range: 1.2x-44.4x).

Reference panel construction

We constructed a reference panel from 250 baboons with putatively unadmixed yellow (N=182) or anubis (N=68) ancestry from captive (SNPRC) and wild east African populations. All were sequenced to high coverage (~37x). Variants were called using GATK, filtered by standard quality criteria, HWE ($p > 10^{-30}$), and genotype missingness (< 0.1). Phasing was performed with SHAPEIT5, incorporating pedigree information from 57 SNPRC individuals.

Genotype imputation

Imputation was performed with GLIMPSE2 following the author-recommended workflow: phasing/imputation per chunk with GLIMPSE2_phase, ligation with GLIMPSE2_ligate, and merging with bcftools. Concordance was estimated using GLIMPSE2_concordance comparing imputed dosages to ground-truth genotypes from high-coverage data.

Benchmarking

We benchmarked GLIMPSE2 against Beagle 5.4, STITCH, and QUILT2 using 36 validation samples downsampled to 1x on chromosome 20 (NC_044995.1). All methods used the same reference panel (where applicable) and ground-truth validation set. Concordance was assessed using GLIMPSE2_concordance for all methods.

Post-imputation filtering

We optimized five serial filters: GEM mappability (score = 1), MAF in reference panel ($\geq$ 0.05), HWE p-value ($\geq$ 0.01), genotype posterior probability ($\geq$ 0.9999), and fraction of missingness ($\leq$ 0.2). Each parameter was tuned to balance concordance improvement against variant retention at 1x coverage.

Mendelian error analysis

Mendelian inconsistencies were quantified using bcftools +mendelian2 across 138 parent-offspring trios (238 unique individuals) identified from field observations and microsatellite genotyping.

Pedigree reconstruction

Pedigree reconstruction was performed with the R package sequoia using 2,505 LD-pruned SNPs (MAF > 0.35, $r^2 < 0.05$, pruned with PLINK v1.9). Age gap distributions were generated from the existing pedigree and used as priors for reconstructing pedigree links.

Data and code availability

The imputation pipeline is implemented as an 11-stage Snakemake workflow with associated R analysis scripts. Sequencing data from previous studies are available under BioProject accessions PRJNA433868 and PRJEB49549. The benchmarking pipeline comparing GLIMPSE2, Beagle 5.4, STITCH, and QUILT2 is available as a reproducible Snakemake workflow.

References

Alberts SC, Altmann J (2012). The Amboseli Baboon Research Project: 40 Years of Continuity and Change. In: Long-Term Field Studies of Primates, Springer, pp. 261-287.

Buckley RM et al. (2022). Best practices for analyzing imputed genotypes from low-pass sequencing in dogs. Mammalian Genome 33:213-229.

Davies RW et al. (2016). Rapid genotype imputation from sequence without reference panels. Nature Genetics 48:965-969.

Davies RW et al. (2021). Rapid genotype imputation from sequence with reference panels. Nature Genetics 53:1104-1111.

Enbody ED et al. (2023). Community-wide genome sequencing reveals 30 years of Darwin's finch evolution. Science 381:eadf6218.

Freudiger A et al. (2025). Estimating realized relatedness in free-ranging macaques by inferring identity-by-descent segments. PNAS 122:e2401106122.

Fuller ZL et al. (2020). Population genetics of the coral Acropora millepora. Science 369:eaba4674.

Hofmeister RJ et al. (2023). Accurate rare variant phasing of whole-genome and whole-exome sequencing data in the UK Biobank. Nature Genetics 55:1243-1249.

Huisman J (2017). Pedigree reconstruction from SNP data. Molecular Ecology Resources 17:1009-1024.

Lin W et al. (2024). Genetic regulatory effects in response to a high-cholesterol, high-fat diet in baboons. Cell Genomics 4:100509.

Robinson JA et al. (2019). Analysis of 100 high-coverage genomes from a pedigreed captive baboon colony. Genome Research 29:848-856.

Rubinacci S et al. (2021). Efficient phasing and imputation of low-coverage sequencing data using large reference panels. Nature Genetics 53:120-126.

Rubinacci S et al. (2023). Imputation of low-coverage sequencing data from 150,119 UK Biobank genomes. Nature Genetics 55:1088-1090.

Sørensen EF et al. (2023). Genome-wide coancestry reveals details of ancient and recent male-driven reticulation in baboons. Science 380:eabn8153.

Stahl K et al. (2021). Assessment of Imputation Quality: Comparison of Phasing and Imputation Algorithms in Real Data. Frontiers in Genetics 12:724037.

Tan HZ et al. (2025). High Imputation Accuracy Can Be Achieved Using a Small Reference Panel. Molecular Ecology Resources 25:e70024.

Vi T et al. (2025). Assessing Genotype Imputation Methods for Low-Coverage Sequencing Data. Molecular Ecology Resources 25:e70049.

Vilgalys TP et al. (2022). Selection against admixture and gene regulatory divergence in a long-term primate field study. Science 377:635-641.

Watowich MM et al. (2023). Best practices for genotype imputation from low-coverage sequencing data in natural populations. Molecular Ecology Resources 25:e13854.