Genotype-level uncertainty and population-level bias in low-coverage imputation: an empirical framework for non-model organisms
Genotype-level uncertainty and population-level bias in low-coverage imputation: an empirical framework for non-model organisms
Abstract
Genotype imputation from low-coverage sequencing (lcWGS) is routine, but the calibration of its uncertainty estimates remains under-investigated in non-model species. We used 36 baboons at 1x coverage to evaluate how population-level quality control (QC) impacts imputation accuracy and the site frequency spectrum (SFS). While GLIMPSE2 and QUILT2 provide high baseline accuracy ( for MAF > 0.02), we show that automated filter optimization increases concordance to . However, we find that GLIMPSE2 is systematically overconfident: a filter retains genotypes that are 97.9% accurate in aggregate, with the lowest-confidence retained genotypes (GP 0.99–0.999) accurate only 92.4% of the time. This miscalibration, combined with the preferential removal of rare variants during QC, distorts the SFS in ways that could bias downstream demographic and selection analyses. To address the technical barriers in non-model genomics, we release a fully automated Snakemake pipeline and a public validation dataset for benchmarking and optimizing post-imputation QC.
Introduction
Genome-wide genotype data are a cornerstone of modern conservation and evolutionary genetics, enabling high-resolution analyses ranging from genome-wide association studies (GWAS) and eQTL mapping to kinship inference and demographic reconstruction. However, generating high-coverage sequencing data for the large sample sizes required to study natural populations remains prohibitively expensive. For a typical mammalian genome, sequencing to 1x coverage costs approximately 540 for 30x coverage (Watowich et al., 2023) — a 30-fold difference that makes population-scale genotyping feasible for non-model species only through low-coverage whole-genome sequencing (lcWGS) combined with genotype imputation.
A transformative opportunity for non-model organism research lies in the increasing availability of public high-coverage datasets. By leveraging existing reference panels from species with similar ancestry, researchers can significantly improve the quality of genotype data from wild populations without the cost of de novo high-coverage sequencing for every study individual. Reference-panel-based tools such as GLIMPSE2 (Rubinacci et al., 2023) and QUILT2 (Davies et al., 2021) have already demonstrated success in diverse systems, including corals (Fuller et al., 2020), Darwin's finches (Enbody et al., 2023), and baboons (Lin et al., 2024).
However, translating these tools into accessible workflows for field biologists and conservation geneticists remains a challenge. Setting up imputation pipelines requires navigating a complex ecosystem of software (e.g., bcftools, GLIMPSE2, Snakemake) and correctly parameterizing concordance evaluations. For researchers without dedicated bioinformatics support, this technical barrier sometimes leads to outsourcing imputation to commercial services (e.g., Fuller et al., 2020), which can function as a black box with limited visibility into the filtering and QC decisions applied to the data.
Furthermore, even when imputation is successfully implemented, two critical methodological nuances remain underexplored. First, we examine how well-calibrated the uncertainty estimates (genotype posterior probabilities, GP) produced by these tools are. Second, we investigate how post-imputation filtering reshapes the allele frequency spectrum (SFS). Because imputation errors are typically enriched at rare variants, aggressive quality control can inadvertently distort the SFS, potentially biasing downstream neutrality tests and demographic inferences. These population-level patterns — including Hardy-Weinberg Equilibrium (HWE) deviations and variant-level missingness — are especially vital to characterize in non-model systems where reference panels may be small or population structures complex.
Here, we address both the scientific and practical gaps by presenting an empirical investigation of genotype-level uncertainty and population-level bias in low-coverage imputation. We introduce a dual-method pipeline (GLIMPSE2 + QUILT2) that utilizes public reference data and automated filter optimization. By characterizing GP calibration and quantifying SFS distortion under different filtering regimes, we provide a framework for maximizing the utility of lcWGS data. The complete pipeline is released as an executable Snakemake workflow with a public test dataset, designed to make high-quality genotype imputation accessible and reproducible for the broader evolutionary biology community.
Results
Both imputation methods achieve comparable baseline accuracy
We imputed 36 baboon samples (Papio cynocephalus P. anubis hybrids) from the Amboseli ecosystem at 1x coverage using a phased reference panel of 250 non-target baboons (182 yellow, 68 anubis; ~37x) on a 5 Mb region of chromosome 20 (NC_044995.1:1-5,000,000). The target population is naturally admixed between two baboon species, and the reference panel was constructed from both parental species — a design that tests imputation under the challenging conditions typical of non-model systems where reference and target populations do not perfectly overlap in ancestry.
After merging per-sample imputed genotypes, monomorphic sites (variants from the reference panel catalog that are not polymorphic in the 36 target samples) were removed, yielding ~47,000 (GLIMPSE2) and ~47,000 (QUILT2) segregating variants for evaluation. Both methods achieved high variant-level concordance on these polymorphic sites (Table 1), with dosage increasing from ~0.85 for rare variants (MAF{\text{REF}} < 0.02) to ~0.94 for common variants (MAF{\text{REF}} 0.40–0.50). The two methods performed comparably across all MAF bins, with no consistent advantage for either.
Table 1. Pre-filter variant-level dosage by MAF bin (5 Mb region, 1x, polymorphic sites only).
| MAF bin | N variants | GLIMPSE2 | QUILT2 |
|---|---|---|---|
| 0–0.02 | 5,791 | 0.842 | 0.864 |
| 0.02–0.05 | 8,028 | 0.911 | 0.908 |
| 0.05–0.10 | 9,278 | 0.915 | 0.917 |
| 0.10–0.20 | 10,296 | 0.926 | 0.932 |
| 0.20–0.30 | 5,941 | 0.929 | 0.938 |
| 0.30–0.40 | 4,316 | 0.921 | 0.926 |
| 0.40–0.50 | 3,423 | 0.937 | 0.944 |
| Total | 47,073 |
Genotype posterior probabilities are systematically overconfident
Before evaluating GP-based filtering, we assessed the calibration of GLIMPSE2's posterior probabilities by binning genotypes by their maximum GP value and computing the fraction that matched the high-coverage truth (Table 2).
Table 2. GP calibration: observed accuracy vs. stated confidence (GLIMPSE2, 10 samples, ~188K genotypes).
| GP bin | N genotypes | Observed accuracy | Expected accuracy |
|---|---|---|---|
| 0.50–0.70 | 3,191 | 0.592 | ~0.60 |
| 0.70–0.80 | 2,087 | 0.672 | ~0.75 |
| 0.80–0.90 | 3,116 | 0.760 | ~0.85 |
| 0.90–0.95 | 3,017 | 0.814 | ~0.93 |
| 0.95–0.99 | 6,235 | 0.853 | ~0.97 |
| 0.99–0.999 | 9,411 | 0.924 | ~0.99 |
| 0.999–1.0 | 161,053 | 0.982 | ~1.00 |
GP values are systematically overconfident across all bins. A GP 0.99 filter retains genotypes from both the 0.99–0.999 bin (9,411 genotypes, 92.4% accurate) and the 0.999–1.0 bin (161,053 genotypes, 98.2% accurate). The aggregate accuracy of the retained set is (9,411 0.924 + 161,053 0.982) / (9,411 + 161,053) = 97.9% — not 99%. The 92.4% figure represents the marginal accuracy of the lowest-confidence retained genotypes (the 0.99–0.999 slice), while the vast majority of retained calls fall in the higher-confidence bin. Although the aggregate miscalibration is modest (~1 percentage point), it reflects a systematic pattern: every GP bin overestimates accuracy, and the lowest-confidence retained genotypes are correct substantially less often than their stated posterior would imply. The miscalibration is consistent with the Li & Stephens model's tendency to overestimate haplotype matching confidence when the reference panel is small or genetically distant from the target population.
This finding has implications beyond our pipeline. Any study that interprets GP = 0.9999 as "99.99% correct" is overestimating genotype quality. The actual quality depends on reference panel size, genetic distance, and local haplotype complexity — factors that GP does not fully capture.
Sequential filter selection via a heuristic score function
We applied four filters sequentially (MAF HWE GP FM), testing a gradient of thresholds at each step. At each step, we selected the threshold maximizing a heuristic score , where is a user-specified parameter that controls the relative penalty for variant loss.
This score function does not claim to optimize any specific downstream statistic. It is a tuning knob that allows researchers to explore the Pareto frontier of accuracy versus retention on their own data. At low (e.g., 0.01), the score is nearly indifferent to variant loss and selects for maximum genotype accuracy. At high (e.g., 1.0), accuracy and retention contribute equally, preserving more of the original variant set. The pipeline exposes as a configurable parameter (filter_alpha in the config file), and we report results at throughout this paper. The table below illustrates how the FM threshold — the final and most impactful filter step — shifts with on our test data:
| FM selected | Variants retained | Retention | Post-filter | |
|---|---|---|---|---|
| 0.01 | 0.10 | 31,893 | 58.7% | 0.981 |
| 0.05 | 0.20 | 40,250 | 74.1% | 0.974 |
| 0.10 | 0.50 | 45,429 | 83.6% | 0.965 |
| 0.50 | 1.00 | 46,546 | 85.7% | 0.962 |
| 1.00 | 1.00 | 46,546 | 85.7% | 0.962 |
Prior to filter optimization, monomorphic sites (AC = 0 or AC = 2N in the imputed population) are removed from the merged VCF. This ensures that the variant set entering the optimizer contains only sites polymorphic in the target population, preventing the reference panel's variant catalog from inflating the apparent variant count or distorting the SFS baseline.
Table 3. Filter optimization results (GLIMPSE2, 54,333 starting variants; QUILT2, 53,353 starting variants).
| Step | Selected (GLIMPSE2) | Selected (QUILT2) | Post-filter variants |
|---|---|---|---|
| MAF | 0.01 | 0.01 | 54,333 / 53,353 |
| HWE | 0.001 | 0.01 | ~54,000 / ~53,000 |
| GP | 0.99 | 0.99 | ~46,500 / ~47,300 |
| FM | 0.2 | 0.2 | 40,250 (74.1%) / 42,943 (80.5%) |
| Post-filter | 0.974 | 0.968 |
The GP filter produced the largest concordance gain for both methods. The GP step also removed variants that became monomorphic after low-confidence genotypes were set to missing — these are singletons or doubletons where the sole carrier had an uncertain genotype call. The improvement from GP filtering is mathematically expected when removing the least confident genotypes. The more informative question is whether the removal process introduces bias in the allele frequency spectrum, which we examine in the next section.
Both methods converged on similar optimal thresholds (MAF 0.01, GP 0.99, FM 0.2), indicating that these filters target shared imputation error modes. GLIMPSE2 achieved slightly higher post-filter (0.974 vs. 0.968) but retained fewer variants (74.1% vs. 80.5%), reflecting its more aggressive GP filtering. HWE had minimal effect at N = 36, consistent with limited statistical power at small sample sizes.
GP + FM filtering distorts the allele frequency spectrum
To assess whether our filtering framework introduces frequency-dependent bias, we compared the site frequency spectrum (SFS) across three datasets: the high-coverage truth genotypes, the raw imputed genotypes, and the filtered output (Table 4).
Table 4. Allele frequency spectrum comparison (proportion of polymorphic variants per MAF bin; GLIMPSE2). With N = 36 diploid individuals, the minimum possible MAF is 1/72 0.014.
| MAF bin | Truth (58,234) | Imputed (54,333) | Filtered (40,250) |
|---|---|---|---|
| 0.01–0.02 | 12.1% | 13.8% | 11.5% |
| 0.02–0.05 | 13.8% | 13.8% | 13.8% |
| 0.05–0.10 | 19.7% | 20.0% | 19.2% |
| 0.10–0.20 | 19.5% | 19.4% | 21.9% |
| 0.20–0.50 | 34.9% | 32.9% | 33.6% |
Monomorphic sites are removed at two stages: during the merge step (variants from the reference panel catalog not segregating in the target population) and after GP filtering (variants whose sole carrier had their genotype set to missing). The resulting imputed SFS closely matches truth across all bins. Post-filter, the lowest bin (singletons/doubletons, MAF 0.01–0.02) decreases from 13.8% to 11.5% — a 16.7% relative depletion — as GP + FM filtering preferentially removes low-confidence genotypes at rare variants. The 0.10–0.20 bin shifts from 19.4% to 21.9%. Absolute shifts are ~2 percentage points per bin, but the relative depletion of rare variants is substantial and its downstream impact depends on the analysis. QUILT2 shows a similar pattern.
Even on polymorphic sites, the GP + FM filters preferentially remove genotypes at rare variants (which have lower posterior confidence), subtly distorting the frequency spectrum toward common alleles. The magnitude of this distortion depends on : at (accuracy-favoring), more rare-variant genotypes are removed; at (retention-favoring), the SFS is better preserved at the cost of lower per-genotype accuracy.
This tradeoff has concrete downstream consequences. For eQTL mapping and GWAS, the 74–81% variant retention at represents a favorable tradeoff — high-quality genotypes yield more statistical power than noisy markers, and the filtered variant density remains sufficient for genome-wide analysis. For demographic inference (Stairway Plot, ), the preferential removal of rare variants could bias effective population size estimates. For neutrality tests (Tajima's ), the spectrum shift toward common alleles would inflate the statistic. For selection scans (iHS), loss of rare-variant haplotype signal reduces detection power. The parameter allows researchers to generate callsets at different points along this tradeoff and evaluate them against their specific downstream requirements, rather than accepting a single filter configuration.
Genome-wide concordance and Mendelian error validation
The 5 Mb results above provide a minimal reproducible example; here we report genome-wide validation on the full dataset, which is too large (~2.5 TB of BAMs across 532 samples and 20 autosomes) for public distribution but whose results confirm the patterns observed in the subset.
When applied to all 20 autosomes of the 36 high-coverage validation samples (coverage = 26.9 8.1 s.d.; range: 20.1–44.4x), the pipeline at 1x coverage using default GLIMPSE2 parameters achieved variant-level concordance (across samples, by MAF bin) of > 0.85 for common variants (MAF{\text{REF}} > 0.05). For the genome-wide application, the downstream target was eQTL mapping, which demands high per-genotype accuracy but tolerates aggressive variant loss — rare variants contribute little to cis-eQTL discovery, and genotype errors at common variants directly reduce statistical power. The optimizer accordingly selected stringent thresholds (GEM = 1, MAF{\text{REF}} \geq 0.05, HWE 0.01, GP 0.9999, FM 0.2), retaining 248,623 biallelic SNPs (39.3% of variable imputed sites on chromosome 20). This aggressive filtering is feasible in baboons because their nucleotide diversity (~0.25%) is roughly 2.5× that of humans, so even after discarding 60% of sites the remaining variant density (~4 SNPs/kb) exceeds what human studies achieve before filtering. Post-filter variant-level concordance exceeded > 0.96 for all MAF bins above 0.05 in the target population. Critically, the same sequential filter strategy — MAF, HWE, GP, FM applied in that order — produced concordance gains at every step in both the 5 Mb subset and the genome-wide analysis, confirming that the approach generalizes beyond the reproducible subset we provide.
As an independent validation that does not require high-coverage ground truth, we assessed Mendelian inheritance consistency using the established Amboseli baboon pedigree (Alberts & Altmann, 2012). The optimized pipeline was applied to all autosomes for all 532 Amboseli samples (36 high-coverage + 496 low-coverage; mean coverage = 5.8 7.2 s.d.; range: 1.2–44.4x). We then checked 259,074 imputed variants on chromosome 20 for Mendelian consistency across 138 full parent-offspring trios (238 unique individuals). The Mendelian error rate was extremely low at the trio level (mean = 0.00094 0.00055 s.d.; range: 0.00017–0.00474), translating to ~1 Mendelian violation per 1,000 variants. Trios with lower mean coverage had slightly higher error rates (Pearson's = −0.33, = 6.47 × 10). At the site level, each variant exhibited a mean of 0.98 ( 0.44 s.d.) errors across the 138 trios, indicating that Mendelian inconsistencies are not concentrated at specific loci but are distributed across the genome at very low frequency.
These genome-wide results confirm that the patterns observed in the 5 Mb reproducible subset — high baseline accuracy, substantial improvement with post-imputation filtering, and robustness across the allele frequency spectrum — generalize to the full genome. The Mendelian error analysis further validates imputation quality across samples of all coverage levels, not just the high-coverage validation subset.
Discussion
Reference panel design for admixed non-model systems. A key challenge in non-model organism imputation is constructing an appropriate reference panel. Our test system — naturally admixed yellow anubis baboons — uses a reference panel drawn from both parental species but explicitly excluding the target population. This design tests imputation under realistic conditions for conservation and evolutionary genetics, where reference panels are often small, drawn from related but non-identical populations, and must capture haplotype diversity across divergent ancestral lineages. The high baseline accuracy achieved here ( for common variants at 1x) demonstrates that reference-based imputation is feasible in admixed non-model systems with modest panel sizes (N = 250), extending previous demonstrations in baboons (Lin et al., 2024).
GP calibration and the meaning of filtering thresholds. Our calibration analysis reveals that GLIMPSE2's GP values are systematically overconfident. A GP 0.99 filter retains genotypes that are 97.9% accurate in aggregate — not 99% — with the lowest-confidence retained genotypes (the 0.99–0.999 bin) accurate only 92.4% of the time. This miscalibration does not invalidate GP-based filtering (it still successfully stratifies genotypes by quality) but it changes its interpretation. A GP 0.99 filter does not retain "99% confident" genotypes; it retains genotypes in the top ~91% of the posterior distribution, the vast majority of which (161K out of 170K) fall in the 0.999+ bin. Users should treat GP thresholds as empirical quality tiers, not probability guarantees. Future work should assess whether GP calibration varies with reference panel size, genetic distance, and local recombination rate.
The parameter is an exploratory heuristic, not a validated optimization. Our score function provides a systematic way to explore the Pareto frontier of accuracy versus retention, but it does not directly optimize any downstream biological statistic. We do not claim that a specific value maximizes GWAS power or minimizes demographic inference bias — such claims would require empirical power analyses that are beyond the scope of this work. Instead, serves as a knob that researchers can turn to generate multiple filtered callsets and evaluate them against their own downstream QC criteria. Our empirical table (Table 3) shows how the accuracy-retention tradeoff shifts across values from 0.01 to 1.0, enabling informed selection.
Population-level signals vs. individual-level evaluation. HWE, GP evaluated across samples, and FM are all population-level metrics — they require multiple samples and capture patterns invisible in single-individual concordance. While HWE filtering had minimal impact in this 36-sample study, its utility is expected to scale non-linearly with sample size, becoming a critical diagnostic for systematic mapping errors in larger cohorts () where it can detect reference allele bias, paralogous mapping, and batch effects. The automated framework accommodates this by testing HWE thresholds rather than hardcoding them — at small N the optimizer correctly selects no HWE filter, while at larger N it can exploit the additional signal.
Method comparison and tool selection for non-model organisms. A practical question for researchers adopting imputation in non-model systems is which tool to use. We compared GLIMPSE2 and QUILT2 — the two leading BAM-native reference-panel-based methods — under identical conditions, providing benchmarking context that is otherwise unavailable for non-human species. Both methods operated directly on sequencing reads, avoiding the information loss inherent in pre-calling genotype likelihoods (as required by Beagle 5.4, which we excluded because this intermediate step reduces accuracy at low coverage; Rubinacci et al., 2023). The two methods achieved comparable baseline accuracy and selected nearly identical optimal GP and FM thresholds after independent filter optimization, indicating that post-imputation error modes are shared across methods rather than tool-specific. For non-model organism researchers, this means the choice between GLIMPSE2 and QUILT2 is less consequential than the choice of reference panel and filtering strategy. Our pipeline runs both methods by default and selects the method that maximizes the heuristic score at the user-specified , or can be configured to run a single method when resources are limited.
On the circularity of filter-driven improvement. Removing low-confidence genotypes is expected to raise concordance; the practical questions are where to set the thresholds and whether the filters are redundant. Our sequential framework shows they are not: GP operates at the genotype level (setting individual uncertain calls to missing), while FM operates at the site level (removing variants where too many samples lost their genotype after GP filtering). These two filters attack poorly imputed genotypes from complementary angles — GP identifies unreliable individual calls, and FM identifies sites where unreliability is pervasive across the population. HWE adds a third angle by flagging sites with systematic deviations that may indicate mapping artifacts rather than true variation. The automated threshold selection removes the subjectivity of choosing cutoffs manually.
Scope and generalizability. The reproducible subset we provide covers 36 samples on a 5 Mb region of a single chromosome — deliberately small to enable rapid validation (~8 minutes). However, the genome-wide results on all 20 autosomes confirm that the same filter strategy and thresholds generalize beyond this window: variant-level concordance improves from > 0.85 to > 0.96 for common variants, and Mendelian error rates across 532 samples and 138 trios remain below 0.1%. We provide the 5 Mb dataset as a minimal reproducible example; the pipeline is designed for whole-genome application. Filter threshold recommendations should nonetheless be revalidated for different sample sizes, coverage levels, or species. The HWE filter, in particular, is expected to become much more powerful at N > 100. Centromeric and other repetitive regions, which we did not test, may behave differently. The pipeline currently requires a pre-phased reference panel; integrating upstream phasing (e.g., SHAPEIT5; Hofmeister et al., 2023) would extend applicability. Finally, our method comparison (GLIMPSE2 vs. QUILT2) is based on a single genomic window — the difference of 0.006 is not statistically significant, and both methods should be considered equivalent in the absence of multi-region validation.
Accessibility. Our pipeline automates the complete workflow from tool installation through dual-method imputation, concordance evaluation, filter optimization, and method comparison. The Snakemake framework ensures reproducibility; the preflight checks diagnose setup problems before computation begins; and the 148 MB test dataset at https://github.com/lycium0605/baboon_imputation_repo enables validation in ~8 minutes. This is designed to reduce the barrier for researchers in non-model systems who need imputed genotypes but lack dedicated bioinformatics support.
Methods
Test data
36 high-coverage baboon samples (Papio cynocephalus P. anubis) downsampled to 1x on a 5 Mb region of chromosome 20 (NC_044995.1:1-5,000,000). Reference panel: 250 non-target baboons (182 yellow, 68 anubis; ~37x), phased with SHAPEIT5 (Hofmeister et al., 2023). Data: https://github.com/lycium0605/baboon_imputation_repo.
Genome-wide validation dataset
For genome-wide validation, we imputed 532 Amboseli baboons comprising 36 high-coverage samples (coverage = 26.9 8.1 s.d.; range: 20.1–44.4x) and 496 low-coverage samples, for a combined mean coverage of 5.8 7.2 s.d. (range: 1.2–44.4x). Imputation was performed across all 20 autosomes using the same GLIMPSE2 pipeline described above. The optimizer selected stringent thresholds appropriate for eQTL mapping: GEM = 1, MAF 0.05, HWE 0.01, GP 0.9999, FM 0.2.
To assess imputation accuracy independently of high-coverage ground truth, we evaluated Mendelian inheritance consistency using 138 full parent-offspring trios (238 unique individuals) identified from long-term field observations and microsatellite genotyping (Alberts & Altmann, 2012). Mendelian error was assessed on 259,074 imputed variants on chromosome 20 using bcftools +mendelian2: trio-level error rates were computed via count mode (-m c), and site-level error distributions were computed via annotate mode (-mode a).
Imputation
GLIMPSE2 v2.0.1: chunking, binary reference splitting, per-sample BAM-native imputation, ligation. QUILT2 v2.0.4: IMPUTE2-format reference, simultaneous multi-sample imputation. Both used the same reference panel and genetic map.
Concordance
GLIMPSE2_concordance with --min-val-gl 0.01, --min-val-dp 5, binned by reference panel allele frequency (--af-tag AF). For GP evaluation, --min-tar-gp computed over genotypes with max(GP) threshold. Weighted mean across all MAF bins was the optimization target.
GP calibration
Genotypes from 10 samples were binned by max(GP) value. For each bin, we computed the fraction matching the high-coverage truth genotype (best-guess from GP vs. diploid GT in truth). Calibration is defined as the agreement between observed accuracy and the GP bin midpoint.
Filter optimization
Sequential optimization (MAF HWE GP FM) with . MAF used reference panel allele frequencies. After GP set genotypes to missing (bcftools +setGT -e 'FMT/GP>={thresh}'), bcftools +fill-tags recomputed population statistics. HWE preceded GP because GP preferentially removes heterozygotes, biasing HWE if applied first.
SFS comparison
For each of truth, raw imputed, and filtered VCFs, we computed minor allele frequency from genotypes and binned into five MAF categories. The truth VCF contains variants polymorphic in the 36 high-coverage samples; the imputed VCF contains all variants from the reference panel catalog.
Data and code availability
Pipeline code (Snakefile, scripts, SKILL.md): https://github.com/lycium0605/genotype-imputation-skill. Test data (36 downsampled BAMs, reference panel, validation VCF): https://github.com/lycium0605/baboon_imputation_repo (187 MB). Tools: GLIMPSE2 (static binaries), QUILT2 (R package), bcftools (compiled with plugins). Automated installation: scripts/install_tools.sh.
References
Alberts, S. C., & Altmann, J. (2012). The Amboseli Baboon Research Project: 40 years of continuity and change. In P. M. Kappeler & D. P. Watts (Eds.), Long-Term Field Studies of Primates (pp. 261–287). Springer. https://doi.org/10.1007/978-3-642-22514-7_12
Davies, R. W., Flint, J., Myers, S., & Mott, R. (2016). Rapid genotype imputation from sequence without reference panels. Nature Genetics, 48(8), 965–969. https://doi.org/10.1038/ng.3594
Davies, R. W., Kucka, M., Su, D., Shi, S., Flanagan, M., Cunniff, C. M., Chan, Y. F., & Myers, S. (2021). Rapid genotype imputation from sequence with reference panels. Nature Genetics, 53(7), 1104–1111. https://doi.org/10.1038/s41588-021-00877-0
Enbody, E. D., Sendell-Price, A. T., Sprehn, C. G., Rubin, C.-J., Visscher, P. M., Grant, B. R., Grant, P. R., & Andersson, L. (2023). Community-wide genome sequencing reveals 30 years of Darwin's finch evolution. Science, 381(6665), eadf6218. https://doi.org/10.1126/science.adf6218
Fuller, Z. L., Mocellin, V. J. L., Morris, L. A., Cantin, N., Shepherd, J., Sarre, L., Peng, J., Liao, Y., Pickrell, J., Andolfatto, P., Matz, M., Bay, L. K., & Przeworski, M. (2020). Population genetics of the coral Acropora millepora: Toward genomic prediction of bleaching. Science, 369(6501), eaba4674. https://doi.org/10.1126/science.aba4674
Hofmeister, R. J., Ribeiro, D. M., Rubinacci, S., & Delaneau, O. (2023). Accurate rare variant phasing of whole-genome and whole-exome sequencing data in the UK Biobank. Nature Genetics, 55(7), 1243–1249. https://doi.org/10.1038/s41588-023-01415-w
Lin, W., Wall, J. D., Li, G., Newman, D., Yang, Y., Abney, M., VandeBerg, J. L., Olivier, M., Gilad, Y., & Cox, L. A. (2024). Genetic regulatory effects in response to a high-cholesterol, high-fat diet in baboons. Cell Genomics, 4(3). https://doi.org/10.1016/j.xgen.2024.100509
Rubinacci, S., Hofmeister, R. J., Sousa da Mota, B., & Delaneau, O. (2023). Imputation of low-coverage sequencing data from 150,119 UK Biobank genomes. Nature Genetics, 55(7), 1088–1090. https://doi.org/10.1038/s41588-023-01438-3
Watowich, M. M., Chiou, K. L., Graves, B., Montague, M. J., Brent, L. J. N., Higham, J. P., Horvath, J. E., Lu, A., Martinez, M. I., Platt, M. L., Schneider-Crease, I. A., Lea, A. J., & Snyder-Mackler, N. (2023). Best practices for genotype imputation from low-coverage sequencing data in natural populations. Molecular Ecology Resources, 25(5), e13854. https://doi.org/10.1111/1755-0998.13854
Reproducibility: Skill File
Use this skill file to reproduce the research with an AI agent.
---
name: genotype-imputation
description: Dual-method genotype imputation pipeline (GLIMPSE2 + QUILT2) for low-coverage WGS data via Snakemake. Runs both methods, optimizes post-imputation filters independently, and selects the best method+filter combination. Use when the user mentions genotype imputation, GLIMPSE2, QUILT2, low-coverage WGS variant calling, imputation accuracy, or filtering imputed variants.
---
# Genotype Imputation from Low-Coverage Sequencing (GLIMPSE2 + QUILT2)
Dual-method pipeline: imputes with both GLIMPSE2 and QUILT2, optimizes filters for each, and selects the best. Works for any diploid organism with a phased reference panel. Validated on baboon data (1x, 250-sample panel). Repo: https://github.com/lycium0605/baboon_imputation_repo
## When the user asks about this pipeline, Claude should:
- **Always confirm** the three required inputs: (1) phased reference panel BCF, (2) sites-only VCF, (3) sorted indexed BAMs
- **Always ask** what organism/chromosome naming convention they use — affects `chrom` and `region` in config
- **Flag immediately** if reference panel N < 50 haplotypes — rare variant imputation will be poor
- **Recommend preflight checks first**: `check_software` → `check_data` → dry run
- **Distinguish** variant-level R² (per MAF bin, across samples) from within-sample R² (across variants) — users frequently conflate these
- **Ask about downstream analysis** — if frequency-sensitive (Tajima's D, SFS, iHS), recommend omitting MAF filter
## Required Tools
| Tool | Version | Purpose |
|------|---------|---------|
| snakemake | >= 7.0 | Workflow engine |
| bcftools | >= 1.17 | VCF manipulation + plugins (fill-tags, setGT) |
| GLIMPSE2 | v2.0.1 | chunk, split_reference, phase, ligate, concordance |
| R + QUILT | >= 2.0.4 | QUILT2 imputation (`install.packages("QUILT")` or conda) |
Quick install for GLIMPSE2 + bcftools: `bash scripts/install_tools.sh`. QUILT2 requires R with the QUILT package installed separately.
**Critical env vars after install:**
```bash
export PATH="$(pwd)/tools:$(pwd)/tools/bin:${PATH}"
export BCFTOOLS_PLUGINS="$(pwd)/tools/libexec/bcftools"
```
## Pipeline Structure
Both methods share preflight checks, chunking, and reference panel AF computation. Then:
- **GLIMPSE2 path:** split reference → per-sample impute+ligate → merge
- **QUILT2 path:** convert ref to IMPUTE2 format → convert genetic map → impute all samples at once → merge
Each method's merged BCF goes through independent concordance evaluation and serial filter optimization. A final comparison step selects the best method+filter combination.
To run only one method, set `methods: ["glimpse2"]` or `methods: ["quilt2"]` in config.
## Execution Order
```bash
snakemake check_software --configfile config.yaml -j 1
snakemake check_data --configfile config.yaml -j 1
snakemake --configfile config.yaml -j 4 -n # dry run
snakemake --configfile config.yaml -j 4 # execute
cat results/pipeline_summary.txt # review comparison
```
## Filter Optimization Logic
Monomorphic sites (AC=0 or AC=2N in the imputed samples) are removed during merge, so only variants polymorphic in the target population enter the filter optimization.
Each method is filtered independently. Filters applied sequentially — each step selects threshold maximizing `score = R² × retention^alpha`, where `alpha` is configurable via `filter_alpha` in config:
| alpha | Behavior | Use case |
|-------|----------|----------|
| 0.01 | Almost pure accuracy | eQTL, fine-mapping |
| 0.05 | Strongly favor accuracy (default) | GWAS, association |
| 0.5 | Balanced | Frequency-sensitive analyses |
| 1.0 | Equal weight accuracy and retention | Preserve SFS for demography |
1. **MAF_REF** — reference panel allele frequency, not imputed-sample MAF
2. **HWE** — before GP, because GP preferentially removes heterozygotes which biases HWE
3. **GP** — genotype-level filter (`bcftools +setGT -e 'FMT/GP>={thresh}'`); evaluated with `--min-tar-gp` for genotype-level R²
- **GP >= 0.99** (typical default): retains more genotypes but includes a marginal tier (0.99-0.999) with ~92.4% accuracy
- **GP >= 0.999**: more stringent, retains only genotypes with >98% observed accuracy, but removes more variants
- Users working with rare variants or frequency-sensitive analyses may prefer GP >= 0.999 for higher per-genotype quality
4. **FM** — variant-level missingness after GP filtering
After GP, `bcftools +fill-tags` recomputes AF, MAF, HWE, F_MISSING from surviving genotypes.
**Outputs per method:**
- `results/{method}/filtered/best_params.json` — selected thresholds
- `results/{method}/filtered/optimization_report.tsv` — full gradient table
- `results/{method}/filtered/filtered_*.vcf.gz` — final filtered VCF
**Comparison output:**
- `results/pipeline_summary.txt` — side-by-side pre/post-filter R² for both methods
- `results/best_method.json` — winning method + parameters
## Alpha Sensitivity Analysis
Users can evaluate different `filter_alpha` values to find the right accuracy/retention tradeoff for their downstream application. The scoring function `score = R² × retention^alpha` shifts behavior from near-pure accuracy (low alpha) to balanced accuracy+retention (high alpha). The pipeline's `filter_alpha` config parameter controls this.
Empirical results from baboon validation data (GLIMPSE2, 1x, 5 Mb test region, MAF>=0.01, HWE>=0.001, GP>=0.99 fixed; FM threshold selected by alpha):
| alpha | FM selected | Variants retained | Retention | Post-filter R² | Recommended use |
|-------|------------|-------------------|-----------|----------------|-----------------|
| 0.01 | ≤ 0.10 | 31,893 | 58.7% | 0.981 | eQTL, fine-mapping |
| 0.05 | ≤ 0.20 | 40,250 | 74.1% | 0.974 | GWAS, association (default) |
| 0.10 | ≤ 0.50 | 45,429 | 83.6% | 0.965 | General-purpose |
| 0.50 | ≤ 1.00 | 46,546 | 85.7% | 0.962 | SFS-aware analyses |
| 1.00 | ≤ 1.00 | 46,546 | 85.7% | 0.962 | Demography, neutrality tests |
GP threshold also interacts with alpha. At 1x coverage on full chromosome 20:
| GP threshold | Genotype retention | Concordance | Dosage R² |
|-------------|-------------------|-------------|-----------|
| none | 100.0% | 95.94% | 0.841 |
| ≥ 0.9 | 95.0% | 97.62% | 0.900 |
| ≥ 0.99 | 89.4% | 98.43% | 0.930 |
| ≥ 0.999 | 84.2% | 98.79% | 0.944 |
| ≥ 0.9999 | 80.8% | 98.93% | 0.951 |
Users should run the pipeline with different alpha values on their own data. Actual values depend on coverage, panel size, and population divergence from the reference.
## Key Diagnostic Decisions
**GLIMPSE2 vs QUILT2 — which is better?**
- Usually GLIMPSE2 wins on accuracy (BAM-native HMM). QUILT2 is competitive and sometimes better for rare variants or when reference panel is large. The pipeline tests both and picks the winner automatically.
**R² lower than expected?**
- Check panel–target overlap: `bcftools isec`
- Confirm panel N > 50 haplotypes
- Check BAM coverage distribution
**GP optimization R² drops?**
- N < 30 → per-variant concordance unreliable. Optimizer correctly picks GP=0.
**No genetic map?**
- `awk -v chr="CHR" '$1==chr {print chr, ".", $2/1e6, $2}' genome.fa.fai > uniform_map.txt`
**QUILT2 fails or is slow?**
- Ensure R + QUILT package installed: `Rscript -e 'library(QUILT)'`
- QUILT2 imputes all samples at once (not per-sample) — memory scales with N×variants
- If QUILT2 is not needed, set `methods: ["glimpse2"]` to skip it
## Common Issues
| Symptom | Cause | Fix |
|---------|-------|-----|
| `check_software` fails on plugin | `BCFTOOLS_PLUGINS` not set | `export BCFTOOLS_PLUGINS=$(pwd)/tools/libexec/bcftools` |
| `check_software` fails on QUILT | R/QUILT not installed | `conda install -c bioconda r-quilt` or `install.packages("QUILT")` |
| GLIMPSE2_phase segfault | Chunk too large for RAM | Reduce chunk size |
| Concordance "No variant found" | Stale `.csi` index | `bcftools index -f <file>` |
| Sample name mismatch | BAM suffix not stripped | Check `bam_suffix` in config |
| GP optimization R² drops | N < 30 samples | Skip GP or use conservative threshold |
| QUILT2 sample names differ | QUILT uses BAM filenames | merge rule strips suffix; verify with `bcftools query -l` |
## References
- `references/pipeline-details.md` — config reference, expected results, runtime
- `references/filter-parameters.md` — filter rationale, gradient defaults, caveats
- `scripts/run_quilt.R` — QUILT2 R wrapper script
- Rubinacci et al. (2023) *Nature Genetics* 55:1088–1096
- Hofmeister et al. (2023) *Nature Genetics* 55:1243–1249
Discussion (0)
to join the discussion.
No comments yet. Be the first to discuss this paper.