Bio Phasing Imputation Imputation Qc

Version Compatibility

Reference examples tested with: bcftools 1.19+, cyvcf2 0.31+, pandas 2.2+, numpy 1.26+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show cyvcf2 pandas numpy then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

The imputation quality field is named differently by each engine: DR2 (Beagle 5.x; AR2 is legacy 4.x), R2 (Minimac4; ER2/EmpRsq for typed sites), INFO (IMPUTE5 and GLIMPSE). bcftools +fill-tags computes AF/MAF/HWE but CANNOT produce an imputation Rsq - the quality number comes from the imputer only. Record the engine, panel, and build, because a quality number is only comparable within the same engine and panel.

Imputation QC -- Filtering on a Self-Estimated Quality, Validating With Masked Truth

"Filter my imputed genotypes by quality and check the accuracy" -> Filter on the imputer's per-variant quality field with a MAF floor, and validate true accuracy by masking - because the routine quality score is the model grading its own posterior, it is confounded with allele frequency, and a flat cutoff silently deletes the rare variants of greatest interest.

• CLI: bcftools view -e 'INFO/DR2<0.3 || INFO/AF<0.01 || INFO/AF>0.99' imputed.vcf.gz (DR2 for Beagle; R2 for Minimac; INFO for GLIMPSE/IMPUTE; Minimac also emits a MAF tag, Beagle only AF)

Scope: QC of already-produced phasing/imputation output - what the numbers mean, which to trust, how to filter, and how poor QC propagates into false GWAS hits. Running imputation -> genotype-imputation. Phasing -> haplotype-phasing. Panel ancestry (which the metric cannot detect) -> reference-panels. The GWAS test on the filtered dosages -> population-genetics/association-testing. VCF field-parsing mechanics -> variant-calling/vcf-statistics. Read-backed phasing switch QC -> long-read-sequencing/haplotype-phasing.

The Single Most Important Modern Insight -- The Quality Score Is the Model Grading Its Own Posterior, Confounded With MAF, So a Flat Cutoff Is a Hidden Rare-Variant Filter

Every routine imputation quality score (IMPUTE INFO, Minimac R2/Rsq, Beagle DR2) estimates the same quantity - the squared correlation r2 between the imputed dosage and the unobserved true genotype - from the posterior spread, without ever seeing truth (Marchini & Howie 2010 Nat Rev Genet 11:499). It is a self-report of confidence, not a measured accuracy. Three facts organize all of QC:

1. The metric is structurally confounded with MAF. The estimator is the fraction of the HWE-expected dosage variance 2p(1-p) the imputed dosages recover; the denominator collapses toward zero as the allele gets rarer, so the estimate is noisier and systematically lower for rare variants. A single flat INFO/R2 >= 0.3 cutoff therefore deletes rare variants far more aggressively than common ones - a hidden MAF filter applied without anyone deciding to. Filter MAF-stratified, or report accuracy per MAF bin.
2. The only true accuracy comes from masking. Hide genotypes that are actually known, re-impute them, and compute the real squared correlation against the held-out truth, binned by MAF (the dosage-r2-by-MAF curve). That is the gold standard; the per-variant R2 is a convenient proxy never validated against truth for that variant. Minimac's EmpRsq (masked, measured) vs Rsq (self-estimated) is the only place the two meet, and a gap between them is a QC alarm (panel/strand/ancestry mismatch). Concordance is disqualified for rare variants: a do-nothing imputer that always calls the major homozygote scores ~98% concordance at MAF 1% (0.99^2) while carrying zero information.
3. The differential-imputation confound manufactures false hits, and every per-group QC metric passes. Because imputation quality is a function of the data and panel, imputing cases and controls (or batches, or ancestries) separately makes the imputation error differ between groups, and a case-control test cannot distinguish that artifactual allele-frequency difference from real association. The fix is structural (impute together / harmonize), not a filter; the absence of a QC red flag never clears it.

The Metrics, Precisely

All three estimate r2 (imputed dosage vs true genotype) from the posterior, never from observed truth. The variance-ratio intuition: estimated r2 = Var(imputed dosage) / [2p(1-p)]. Under perfect information the dosages equal the true 0/1/2 genotypes and recover the full HWE variance (ratio 1); under no information every dosage collapses to the mean 2p and the variance goes to 0 (ratio 0). The fraction of HWE variance recovered IS the estimate, which is also why it is noisier and lower at low MAF. That Var(dosage)/2p(1-p) form is specifically the Minimac/Beagle estimator (the spread of the point dosages across samples); IMPUTE INFO targets the same r2 but computes it differently, averaging each sample's within-individual posterior variance - the mechanistic reason the three numbers are not interchangeable across engines.

The provenance matters: the field name tells the tool, the numbers are NOT comparable across engines (a Beagle DR2 of 0.8 is not a Minimac R2 of 0.8), and bcftools +fill-tags cannot produce any of them.

The INFO x MAF Interaction and Thresholds

At low MAF there is little dosage variance to predict and few panel copies of the rare haplotype to anchor the estimate, so R2 is intrinsically noisier and downward-biased there - a property of the construction, not a bug. The standard cutoffs and their status:

True Accuracy by Masking

The accepted accuracy curve is aggregate r2 (squared Pearson correlation between imputed dosage and the masked-then-revealed true genotype) binned by MAF; it decreases monotonically as MAF falls. Workflow: mask typed genotypes (array sites or a held-out set), re-impute from the panel, compute squared correlation vs the withheld truth, bin by MAF. Leave-one-out (one variant at a time) is the per-variant version and is what produces Minimac's EmpRsq for typed sites. Never use concordance as the rare-variant accuracy metric (Ramnarine 2015 PLoS One 10:e0137601); it is dominated by the major-homozygote class and inflates rare-variant accuracy.

MAF-stratified quality summary

Goal: Reveal the hidden-MAF-filter effect by reporting imputation quality per MAF bin instead of as one global mean, so the rare-variant tail a flat cutoff would delete is visible.

Approach: Parse the imputed VCF for the engine's quality field and allele frequency with cyvcf2, derive MAF, bin it, and report mean quality and the fraction passing a candidate cutoff per bin.

import numpy as np
import pandas as pd
from cyvcf2 import VCF

def quality_by_maf(vcf_path, qual_key='DR2', cutoff=0.3):
    rows = []
    for v in VCF(vcf_path):
        q = v.INFO.get(qual_key)
        af = v.INFO.get('AF')
        if q is None or af is None:
            continue
        af = af[0] if isinstance(af, tuple) else af
        rows.append((min(af, 1 - af), float(q)))
    df = pd.DataFrame(rows, columns=['maf', 'qual'])
    bins = [0, 0.001, 0.01, 0.05, 0.5]   # rare-to-common; the rare bins are where a flat cutoff bites
    df['maf_bin'] = pd.cut(df['maf'], bins=bins)
    summary = df.groupby('maf_bin', observed=True).agg(n=('qual', 'size'), mean_qual=('qual', 'mean'), frac_pass=('qual', lambda q: (q >= cutoff).mean()))
    return summary

quality_by_maf('imputed.vcf.gz', qual_key='DR2', cutoff=0.3)

Phasing QC -- Switch Error Rate

Switch error rate (SER) = switch errors / opportunities, an opportunity being each consecutive heterozygous-site pair, scored against trio/duo/benchmark truth. A flip error is two switches one site apart (an isolated mis-assignment, not a long-range switch); Hamming distance counts overall haplotype differences and inflates on block swaps. Trio-based SER is biased upward by genotype error, which is why SHAPEIT5's switch reports SER and a genotyping-error rate jointly. Magnitudes are always MAC- and N-stratified (sub-0.5% common-variant SER for modern tools; single-digit-percent and rising as MAC approaches 1). Tools: whatshap compare, vcftools --diff-switch-error, SHAPEIT5 switch.