How I built a dual-condition genotype imputation benchmark on 1000 Genomes data, and what I found: Beagle's DR2 calibration degrades 8-fold under reference-panel ancestry mismatch, even after abstaining on half the genome.
- Beagle's DR2 confidence score is well calibrated when the reference panel matches the target ancestry (ECE 0.026) and badly miscalibrated when it does not (ECE 0.205). Same people, same chip, same software; only the panel changed.
- A "trust DR2 >= 0.8" rule tuned on the matched case abstains on 54% of the genome under mismatch and the surviving variants still deliver r2 of 0.67 against a promised 0.8. Abstention makes the failure quieter, not smaller.
- Global confidence thresholds on imputed data quietly hand minority-ancestry cohorts worse genotypes than the threshold implies. Fixing it needs ancestry-aware calibration, not stricter global filters.
- Everything is public: github.com/gbadedata/imputation-calibration-benchmark: 63 tests, exact accounting, results tables generated from the machine-readable report.
Genotyping arrays are cheap because they read about 700,000 positions and let software guess the rest. The guessing is called imputation: match each person's chip data against a reference panel of fully sequenced genomes, infer the unread positions, and attach a confidence score to every guess. Beagle calls its score DR2.
Downstream systems trust DR2. Recall-by-genotype studies recruit real people based on imputed variants. Polygenic risk scores are computed on imputed dosages. Filters like "keep variants with DR2 above 0.8" appear in methods sections everywhere, applied globally, as if the score means the same thing for everyone.
Here is the problem: most reference panels are European-heavy, and most confidence thresholds are tuned on European-majority data. If DR2 is calibrated for one ancestry and quietly miscalibrated for another, every genotype-conditioned decision inherits that bias. Silently, because a confidence score that overpromises does not announce itself.
So I measured it. Here is what I built, what broke, and what the numbers mean.
The design
1000 Genomes Phase 3 provides fully sequenced truth genotypes for 2,504 people. That makes a controlled experiment possible: simulate a genotyping chip by masking each target's genome down to an array-density site set, impute the masked positions, then grade every imputed genotype against the withheld truth.
The experiment has two arms that differ in exactly one variable:
| Condition A (control) | Condition B (mismatch) | |
|---|---|---|
| Targets | 100 African-ancestry samples | the same 100 samples |
| Reference panel | 533 AFR samples | 498 EUR samples |
| Chip data | identical 17,698 QC-passed sites | identical |
| Software, map, settings | identical | identical |
Same people, same chip, same Beagle. Only the panel changes. Any difference in how honest DR2 is between the two arms is attributable to ancestry mismatch.
The design in one picture. The withheld truth is the answer key; Beagle never sees it.
The analysis then asks two questions. First, calibration: within each DR2 bin, does the promised accuracy match the empirical squared correlation between imputed dosage and truth? Second, abstention transfer: choose a trust threshold on Condition A using a stated rule, apply that identical threshold to Condition B, and report what actually happens.
Every number carries a bootstrap confidence interval. The threshold rule is written down before the results exist. Whatever comes out, ships.
Three things that broke, one hypothesis that died
Real engineering does not go according to plan, and the failures are usually more instructive than the successes.
1. The 1000 Genomes v5b chr20 release carries no rsIDs
The original design intersected array sites with the HapMap3 SNP list by rsID, the standard move. Then the site partition returned zero array sites, with all 1,739,315 records classified missing-ID. The raw file confirmed it: every ID column entry is a dot. The v5b re-release stripped rsIDs entirely, something no tutorial mentions because every tutorial references the older files.
Bonus trap: the widely copied filename for this file (v5a) returns 404. Only listing the live server directory revealed the truth. Fix: deterministic thinning of common variants to 20,000 chr:pos-keyed sites, a fallback authorized in the build spec before it was needed. Lesson: upstream re-releases change content, not just filenames, and ID-dependent designs must verify ID presence before building on them.
2. plink2's BGEN export is ref-last, and the wrong import mode mirrors your genome
I validate every format conversion with a content-signature round trip: VCF to BGEN to VCF, compared by per-site ALT allele count. The first run failed with a perfect mirror: every count c became 2N minus c across all 13,367 sites. That signature diagnoses an allele-order flip exactly: the export wrote alleles ref-last, the import assumed ref-first, and every genotype was read through a mirror. One flag fixed it, and the verifier now recognises the complement pattern and names the fix in its own error message:
# Diagnose the known allele-order failure: every count complemented to 2N - c
if all(
a.pos == b.pos and b.alt_count == 2 * a.n_samples - a.alt_count
for a, b in zip(original, round_tripped, strict=True)
):
raise RoundTripError(
"ALLELE-ORDER FLIP: every alt count is complemented (c -> 2N - c). "
"The BGEN import REF/ALT mode is wrong - plink2 a7.4 bgen-1.2 export is "
"ref-LAST; import must use 'ref-last'."
)
Lesson: "the command exited 0" is not verification. Formats need content-level checks, and a failed check with a diagnostic signature is worth more than a lucky pass.
3. My dropped-sites hypothesis was falsified
I expected the EUR panel to force Beagle to drop AFR-common sites that are monomorphic in Europeans, predicting a coverage gap between conditions. Wrong: coverage was 100 percent in both arms. Beagle imputes every reference marker, including panel-monomorphic ones; it just imputes them badly. The validation gate I had built asymmetric coverage floors for turned out unnecessary, and the mismatch effect showed up entirely where the calibration analysis looks. I kept the falsified prediction in the repo history. Being wrong in a measurable way is the point of building instruments.
And one silent echo worth recording: chip QC on the 100 AFR targets removed 2,280 of 20,000 array sites (11.4 percent) for low within-cohort MAF: sites common globally but rare in this cohort. At evaluation, another 19,093 sites (7.8 percent) were monomorphic within the targets and excluded with a count. The frequency divergence between "common in the world" and "common in this cohort" appears twice before imputation is even graded. It is the same story the headline result tells, visible earlier and cheaper.
The numbers
| Metric | A: matched (AFR panel) | B: mismatched (EUR panel) |
|---|---|---|
| Expected calibration error (ECE) | 0.026 | 0.205 |
| Fraction abstained (shared DR2 >= 0.8 rule) | 12.9% | 54.1% |
| Mean empirical r2 among retained variants | 0.940 | 0.671 |
| Evaluation variants analysed | 226,943 | 226,943 |
| Evaluation-site coverage | 100% | 100% |
| Wall-clock, full chr20 imputation | 51 s | 42 s |
Pipeline scale for context: 1,739,315 raw biallelic SNPs, 20,000 simulated array sites, 17,698 after chip QC, 33 samples removed by KING kinship filtering across all AFR and EUR candidates before the target draw. 63 automated tests, exact accounting at every stage: every excluded site and sample is counted, never silently dropped.
What the numbers actually mean
Figure 1. The dotted diagonal is perfect honesty (promised accuracy = delivered accuracy). Vertical distance below it IS the overconfidence. The dash-dot line is the abstention threshold chosen on A and transferred to B; small multiples repeat the comparison within MAF strata.
ECE 0.026 versus 0.205 is an 8-fold calibration collapse. On the matched panel, DR2's promise and the delivered accuracy differ by about 3 points on average. Under mismatch, the average shortfall is 20 points. Same software, same people, same sites.
The abstention result is the operational punchline. DR2 is not oblivious to mismatch: under the identical trust rule, Condition B abstains on 54 percent of variants versus 13 percent for A. Relative ranking partially survives. But absolute calibration does not: the variants B still vouches for deliver r2 of 0.671 against a promised floor of 0.8. Read that again, because it is the finding: the confidence score overpromises even after it has already thrown away half the genome. An abstention policy calibrated on the majority-ancestry case does not protect the minority-ancestry case. It just fails more quietly.
Figure 2. What accuracy each condition can buy at each abstention price. Condition B cannot reach Condition A's retained accuracy at any threshold on this sweep; the ceiling itself is lower.
Why this happens is a modelling point, not a data-volume point. DR2 is Beagle's model-internal accuracy estimate, derived from posterior probabilities under its haplotype model. Those posteriors are only as honest as the assumption that the panel spans the target haplotypes. A EUR panel lacks AFR haplotype backbones and carries different LD structure, so the model confidently matches targets to the wrong haplotype copies. More data from the same panel would not fix it. Representativeness would.
Design decisions worth understanding
Why thin the array from common variants only? Real chips are deliberately common-biased, and an "array" of rare singletons would give the phasing model nothing to work with. The thinning pool is restricted to MAF >= 0.01 before deterministic even-spaced selection, reproducible from the repo with no seed.
Why chip QC on the targets and not the panels? Chip QC belongs to the chip cohort; that is what a production workflow QCs. Computing HWE across mixed AFR and EUR panels would trigger spurious Wahlund-effect failures, and MAF-filtering the array by EUR panel frequencies would delete exactly the AFR-common sites where mismatch acts, biasing the experiment toward its own null. The original spec had this wrong; the redesign is documented in the repo with the reasoning.
Why choose the abstention threshold with an "every bin above must qualify" rule? The naive rule (smallest bin whose CI lower bound clears the floor) steps over non-monotone dips. If bin 0.6 qualifies but bin 0.7 dips below the floor, naive selection trusts a region it should not. There is a unit test that constructs exactly this trap and proves the rule refuses it. The whole rule is eight lines:
def select_threshold(calib_a: CalibrationResult, accuracy_floor: float) -> float | None:
"""Smallest bin lower-edge with all bins at/above clearing the CI-lower floor."""
bins = calib_a.bins
for i, b in enumerate(bins):
if all(bb.r2_ci_low >= accuracy_floor for bb in bins[i:]):
return b.lo
return None # nothing qualifies: abstain on everything, reported as such
Why bootstrap everything? Because a headline number without an interval is noise wearing a suit. Every bin statistic and every retained-accuracy figure carries a 1,000-rep percentile CI, seeded and deterministic.
The limitations, stated plainly
Within-cohort circularity, and it cuts both ways. Condition A's panel and targets come from the same cohort: same platform, same joint calling, shared LD structure even after kinship filtering. A's absolute numbers are a favourable case, not a real-world estimate, and part of the A-versus-B contrast may reflect cohort effects beyond ancestry per se. The clean follow-up is an external AFR panel for Condition A. The claim this design supports is the contrast, run under controlled conditions, not the absolute magnitudes.
Chromosome 20 is a model system, sized for a 32 GB laptop. Per-variant confidence calibration is measurable on one chromosome. Nothing here is a polygenic score result, and the repo says so in its first screen.
The simulated array mimics density and common-variant bias, not any commercial chip's site-selection strategy.
What this looks like in a real genomics organisation
Any organisation that selects people by imputed genotype (recall-by-genotype studies, PRS-tiered recruitment, variant-conditioned trial eligibility) and applies one global DR2 threshold is handing minority-ancestry participants quietly worse data than the threshold implies. The mitigation is not a stricter global threshold; Condition B already abstains on half the genome and still underdelivers. It is ancestry-aware calibration: per-ancestry thresholds, or empirical recalibration of DR2 against a held-out truth subset per cohort. This benchmark measures the size of the problem on a controlled system, and the machinery to run the measurement (the accounting, the round-trip checks, the bootstrap harness) is exactly what the production version would need.
Methods at a glance
| Data | 1000 Genomes Phase 3, chr20, release 20130502 (v5b) |
| Targets / panels | 100 AFR / 533 AFR (A) / 498 EUR (B), KING-filtered at 0.0884 |
| Array simulation | 20,000 common-variant sites, deterministic thinning; 17,698 after chip QC |
| Imputation | Beagle 5.5 (27Feb25.75f), GRCh37 plink map, default seed |
| Accuracy metric | per-variant r2, imputed dosage vs withheld truth, 226,943 variants |
| Uncertainty | 1,000-rep seeded percentile bootstrap on every reported statistic |
| Tooling | bcftools, plink2 alpha 7.4, Python 3.12, Nextflow DSL2 wrap, 63 pytest tests |
| Hardware | one laptop: WSL2, 32 GB RAM, ~1 minute per full-chr20 imputation |
Where this sits in the literature
The unequal performance of genomic tools across ancestries is well established: Martin and colleagues showed in 2019 that polygenic scores trained on European cohorts lose substantial accuracy when applied to African-ancestry individuals, and the field's answer has been bigger, more diverse panels, with the TOPMed program and its imputation server as the flagship. This benchmark asks a narrower and less-examined question: not whether mismatch reduces accuracy (it does, and that is known), but whether the confidence score remains honest about the reduction. The distinction matters operationally, because pipelines rarely consume accuracy directly; they consume DR2 thresholds. A score that degraded but stayed calibrated would be self-correcting under any threshold. The measured result is that it does not stay calibrated, which is precisely the failure mode global thresholds cannot see.
Reproduce it
git clone https://github.com/gbadedata/imputation-calibration-benchmark
cd imputation-calibration-benchmark
python3 -m venv .venv && source .venv/bin/activate
pip install -e ".[dev]"
ruff check . && pytest # 63 tests, no data needed
# place raw inputs per data/PROVENANCE.md (URLs, checksums, versions), then:
bash scripts/phase1.sh # derivation, exact site accounting
bash scripts/phase2.sh # kinship, chip QC, BGEN round trip
bash scripts/phase3.sh # dual-condition Beagle, ~1 min each
bash scripts/phase4.sh # calibration, abstention, figures, report
Total compute is a laptop-hour; the slowest step is downloading the 312 MB chr20 VCF.
Repository
Everything in this post is public, with the full verification log including both upstream surprises:
github.com/gbadedata/imputation-calibration-benchmark
The README leads with the limitations, the results tables are generated from the machine-readable report rather than hand-typed, and the figures (the design schematic, the calibration curves, and the abstention trade-off) are rendered inline.
What is next
The obvious follow-up is replacing Condition A's within-cohort panel with an external ancestry-matched reference, which converts the favourable-case control into a realistic one. Beyond that, this calibration-and-abstention lens extends naturally to somatic variant calling at low allele fractions, where the cost of a confidence score that overpromises is measured in missed cancer variants rather than miscounted genotypes.
Working on genotype imputation, statistical genetics, or calibration in genomics pipelines? Connect on GitHub or LinkedIn.
Top comments (0)