Bayesian PGx allele-frequency estimation

Problem

ACMG variant classification assigns benign evidence using hard allele-frequency thresholds: BA1 at 5% and BS1 at 1%. A single point estimate is compared against the threshold and a binary call follows.

When a sample size is small relative to the distance from the threshold, that point estimate gives false confidence. A frequency near the cutoff can flip classification with a handful of observed alleles, yet the binary call hides how uncertain the crossing is.

This bites hardest for populations underrepresented in gnomAD, where the reference allele numbers driving the estimate are thin. The project targets 9 pharmacogenomic loci across 7 African populations, drawing allele counts from 661 whole-genome samples in the AfriPharmaGen Catalog.

How it works

The pipeline models each locus with a two-level Beta-Binomial hierarchy fitted by MCMC in PyMC. A continental frequency is drawn from a Beta prior; subpopulation frequencies are drawn from a Beta centered on it with a Gamma-distributed concentration controlling how tightly populations cluster; observed alt-allele counts enter through a Binomial likelihood.

gnomAD v4 African frequencies and allele numbers become informative priors, downweighted by a factor of 10 so the prior informs without dominating the observed data. Loci absent from gnomAD fall back to a weakly informative prior.

For each of the 63 locus-population pairs the posterior yields the median, a 95% credible interval, and the two classification probabilities P(frequency > BA1) and P(frequency > BS1). These map to zones (confident BA1, confident BS1, below both, uncertain), replacing the binary threshold call with a graded one. A second entry point renders six publication figures in PNG, PDF, and SVG.

gnomAD v4 priors (downweighted) --\
                                  v
AfriPharmaGen allele counts --> Beta-Binomial hierarchy (PyMC MCMC)
                                  v
            Posterior per locus-population pair
         P(>BA1), P(>BS1)  |  Wilson CI comparison  |  threshold sensitivity
                                  v
            Zone classification --> disagreement set --> CPIC actionability

Hard parts

• Encoding gnomAD reference data as a prior without letting it swamp small local samples: the effective allele number is divided by a downweight factor, and missing loci default to an uninformative prior.
• Sharing information across related populations through the hierarchical layer, so a sparsely sampled population is pulled toward the continental estimate.
• Turning a binary threshold decision into a calibrated one by computing posterior tail probabilities against both BA1 and BS1 and defining explicit confident, below-both, and uncertain zones.
• Checking the Bayesian posterior adds information over a frequentist baseline: Wilson score intervals are computed for every pair and compared for narrower width, shifted median, and threshold straddle.
• Sweeping the classification confidence from 0.80 to 0.99 to see how the disagreement count moves, with posterior predictive checks and R-hat / effective-sample-size convergence monitoring.

Results

• Applied to 9 loci across 7 African populations, producing 63 locus-population classifications from 661 whole-genome samples.
• 19 of 63 classifications (30.2%) disagree between the point-estimate and Bayesian approaches; the Bayesian method labels these as uncertain.
• 10 of the 19 disagreements (53%) involve alleles with CPIC Level A evidence and confirmed functional impact, including CYP2C9*8 (warfarin), CYP2B6*18 (efavirenz), and CYP2D6*29 (opioids).
• CYP2C19*9 is flagged uncertain for BS1 across all 7 populations; CYP2C9*8 is uncertain for BA1 in 6 of 7.
• The full pipeline runs in about 20 seconds for all 9 loci at 4 chains and 7000 draws each.

Artifacts

• A single-command analysis pipeline: hierarchical fit, classification, posterior predictive checks, Wilson comparison, threshold sensitivity, and CPIC actionability stratification.
• A figure generator that regenerates all six figures in PNG, PDF, and SVG.
• Four JSON result files: Bayesian classifications, Wilson comparison, threshold sensitivity, and clinical actionability.
• Public GitHub repository under MIT license with a companion manuscript archived on Zenodo.