bio-clinical-databases-gnomad-frequencies

## Version Compatibility Reference examples tested with: requests 2.31+, hail 0.2.130+, pandas 2.2+, myvariant 1.0+. Current gnomAD release is **v4.1 (May 2024)**; v4.1 fixed the v4.0 AN under-counting issue that inflated rare-variant AF estimates by 5-10%. Before using code patterns, verify installed versions match. If versions differ: - Python: `pip show <package>` then `help(module.function)` to check signatures - Hail: `hl.version()`; pin to >=0.2.130 for v4 schema If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. The gnomAD browser GraphQL API at `https://gnomad.broadinstitute.org/api` is the supported public endpoint; Hail Tables on Google Cloud Storage at `gs://gcp-public-data--gnomad/` are the supported bulk access. # gnomAD Frequency Queries and Constraint **'How rare is this variant in the general population?'** -> Pull allele frequency, grpmax FAF95 (the ACMG-grade frequency), LOEUF gene-level constraint, structural variant catalog, mtDNA frequencies, and the appropriate dataset version per use case. - Python (single variant): GraphQL via `requests.post('https://gnomad.broadinstitute.org/api', json={'query': ..., 'variables': ...})` - Python (aggregator): `myvariant.MyVariantInfo().getvariant(hgvs, fields=['gnomad_exome', 'gnomad_genome'])` - Python (bulk): `hl.read_table('gs://gcp-public-data--gnomad/release/4.1/ht/exomes/gnomad.exomes.v4.1.sites.ht')` ## v2.1.1 / v3.1.2 / v4.x: When to Use Which This is the most consequential decision in any gnomAD query. The releases are **not interchangeable**; choice determines what can and cannot be said about a variant. | Release | Build | Samples | Use when | Fails when | |---------|-------|---------|----------|-----------| | **v2.1.1** | GRCh37 | 125,748 exomes + 15,708 genomes | Constraint metrics needed (LOEUF v2 most-validated); chrX/Y constraint required; GRCh37 native non-negotiable | GRCh38 native cohort; modern rare-variant FAF95 (use v4) | | **v3.1.2** | GRCh38 | 76,156 genomes (NO exomes) | Non-coding region rare variants on GRCh38; mtDNA frequencies | Exome variants needed (no exomes); 76k cohort smaller than v4 | | **v4.0/v4.1** | GRCh38 | 730,947 exomes + 76,215 genomes = **807,162 total** | Default for everything; rare-variant filtering, FAF95, gene queries | chrX/Y constraint (not released); cancer-cohort analysis (no TCGA in v4) | **Critical caveats:** - v4 genomes are the SAME 76,215 v3 samples reprocessed against GRCh38 with updated pipelines; not independent. - ~81% of v2 genomes are also in v3; joint v2+v3 meta-analysis must dedupe at sample ID. - v4 includes 416,555 UK Biobank exomes under a specific collaboration agreement; check use terms. - **v4 does NOT include TCGA**, so the `non_cancer` subset is unnecessary; the v4 subset is `non_ukb` (excludes UKB exomes for ancestry rebalancing). - Liftover v2 (GRCh37 -> GRCh38) is NOT equivalent to v4 native; variant representation differs at ~0.5-1% of sites due to assembly fixes. ## v4 Ancestry Groups: popmax -> grpmax Terminology v4 ancestry groups: **AFR, AMR, ASJ, EAS, FIN, MID, NFE, SAS, AMI, REMAINING**. The **MID (Middle Eastern) group was new in v4**; previously absorbed into "OTH". The **REMAINING** group (31,256 v4 samples) is individuals who did not cluster with any reference; they contribute to overall AF but not to grpmax. **Terminology shift:** gnomAD documentation and ACMG-facing narrative uses **grpmax** (genetic ancestry group max) -- replacing the older **popmax** ("population max") term -- to disambiguate genetic ancestry from self-reported race/ethnicity. The public GraphQL schema still exposes legacy field names containing `popmax` (e.g. `faf95.popmax`, `faf95.popmax_population`); these are the grpmax values under the modern terminology. Always check the schema version when writing queries; new browsers may rename these fields. `grpmax_faf95` is the operational ACMG field. It computes the maximum 95% lower-CI allele frequency, **excluding bottleneck groups** (AMI, ASJ, FIN, REMAINING) because pathogenic founder variants in those groups would otherwise falsely trigger BS1/BA1. MID is included in grpmax but is the smallest non-bottleneck group with highest per-allele variance. ## Filtering Allele Frequency (FAF95): The ACMG-Grade AF Whiffin 2017 *Genet Med* 19:1151 introduced FAF95 = Poisson lower bound of 95% CI for AF. By construction, AF > FAF95; FAF95 is the conservative frequency for ACMG application. **Max-credible-AF formula:** `(prevalence x heterogeneity x allelic-contribution) / (penetrance x 2)`. Plug in disease parameters to get the gene-specific BA1 / BS1 threshold; compare against `grpmax_faf95`. | Code | Threshold | Notes | |------|-----------|-------| | **BA1** | AF > 5% in any non-bottleneck group | ClinGen SVI default; VCEPs may override (Hearing Loss VCEP uses 0.5%) | | **BS1** | AF > gene-specific max-credible-AF | Computed per gene via Whiffin formula | | **PM2_Supporting** | Absent or ultra-rare in gnomAD | Downgraded from PM2_Moderate in SVI 2020 | Use `grpmax_faf95`, not raw AF, for BS1/BA1 application; this is the ClinGen-recommended approach. ## Constraint Metrics: pLI, LOEUF, missense Z Karczewski 2020 *Nature* 581:434 defined LOEUF as the upper bound of the 90% CI of observed/expected pLoF count per gene. LOEUF is **recommended over pLI** because it is continuous and accounts for gene size more rigorously. | Metric | What | Interpretation | |--------|------|----------------| | **LOEUF** | Upper bound of 90% CI of LoF observed/expected ratio | Lower = more LoF-intolerant; **first decile (LOEUF < 0.35 v2; < 0.6 v4) = strongly intolerant** | | **pLI** | Probability LoF intolerant | Still used; gnomAD team recommends LOEUF for ranking | | **Missense Z** | Z-score of observed-vs-expected missense | Z > 3.09 = 'constrained' set (~p < 0.001, one-tailed) | | **Missense O/E** | Observed/expected missense ratio | Continuous form of missense Z | **Critical version mismatch:** - v2.1.1 constraint metrics published 2020; v4 constraint published **March 2024** (4 months after v4 data release). - **v4 constraint is autosomes only; chrX and chrY constraint metrics in v4 are NOT released**. For X/Y constraint, fall back to v2.1.1. - **LOEUF first decile shifted v2 to v4**: v2 < 0.35; v4 < 0.6 (larger sample shifted the distribution). Gene rank in deciles is stable across versions but absolute thresholds are NOT interchangeable. ## Subsets: non_cancer, non_neuro, controls | Release | Subset | Removes | Use when | |---------|--------|---------|----------| | v2.1.1 | `non_cancer` | TCGA | Cancer-related variant analysis (avoids circularity) | | v2.1.1 | `non_neuro` | Psychiatric/neuro cohorts | Neuropsychiatric variant analysis | | v2.1.1 | `controls` | Cases with known disease (~60k samples) | Disease-association calibration | | v3.1.2 | `non_v2` | v2 overlapping samples | Independent of v2 | | v3.1.2 | `controls_and_biobanks` | Disease cases retained, biobanks emphasized | Population-level reference | | v4 | `non_ukb` | UK Biobank exomes | When EUR-skew of UKB problematic | | v4 | `non_neuro` | Deprecated | -- | | v4 | `non_cancer` | Unnecessary (no TCGA in v4) | -- | ## SV Catalog and CNV | Resource | Release | Samples | Coverage | |----------|---------|---------|----------| | gnomAD-SV v2 | Collins 2020 *Nature* 581:444 | 14,891 unrelated WGS | 433k SVs, GRCh37 | | gnomAD-SV v4 | Nov 2023 | 63,046 unrelated WGS | 1,199,117 high-confidence SVs, GRCh38 | | gnomAD-CNV v4 | Nov 2023 | 464,297 individuals (exome-derived gCNV) | Rare (AF < 1%) autosomal coding CNVs | gnomAD-CNV v4 is the resource that democratized exome-derived CNV background frequencies; previously only ExAC-CNV provided this at scale. ## mtDNA (Laricchia 2022 *Genome Res* 32:569) 10,850 unique mtDNA variants across 56,434 individuals (v3.1). Frequencies reported per nuclear-ancestry AND per mitochondrial-haplogroup. Heteroplasmy >=10% threshold; ~1/250 individuals carry pathogenic mtDNA variant at heteroplasmy >=10%. mtDNA inheritance is non-Mendelian; standard ACMG criteria do not apply directly; use MITOMAP and HmtVar in parallel. ## VEP Version Pinning Each gnomAD release pins to a VEP version: - **v4 uses VEP 105** with GENCODE 39 / Ensembl 105 transcripts - v2.1.1 uses VEP 85 A variant's consequence prediction can flip between v2 and v4 due to MANE Select adoption and transcript-set updates. Always pin VEP version when reproducing gnomAD annotations. ## Decision Tree by Query Scenario | Scenario | Recommended path | Why | |----------|------------------|-----| | Single variant AF lookup | GraphQL API or myvariant.info | Lowest latency; returns full per-ancestry breakdown | | ACMG BS1/BA1 application | `grpmax_faf95` from v4 | The ClinGen-recommended field | | Gene-level LoF constraint (autosomes) | LOEUF from v4 March 2024 release | Larger sample, more stable | | Gene-level LoF constraint (chrX/Y) | LOEUF from v2.1.1 | v4 X/Y constraint NOT released | | Bulk rare-variant filter (cohort-scale) | Hail Table on GCS | No rate limits; full schema | | SV frequency | gnomAD-SV v4 (WGS) or gnomAD-CNV v4 (exome) | Choose by data type | | mtDNA frequency | v3.1 mtDNA release (Laricchia 2022) | Only gnomAD release with mtDNA | | Cancer-variant analysis | v2.1.1 `non_cancer` subset OR v4 (no TCGA) | Avoid TCGA circularity in v2 | | Comparison across builds | Use canonical SPDI or CA ID, normalize first | Liftover != native | ## Single Variant Query (GraphQL) **Goal:** Retrieve exome + genome AF, grpmax, FAF95, and per-ancestry breakdown for one variant. **Approach:** Hit gnomAD's GraphQL API with explicit dataset version; parse the nested response. ```python import requests GNOMAD_API = 'https://gnomad.broadinstitute.org/api' def query_variant(chrom, pos, ref, alt, dataset='gnomad_r4'): '''Query gnomAD GraphQL for variant frequency + grpmax FAF95. dataset options: gnomad_r4 (v4.1, default), gnomad_r3, gnomad_r2_1 ''' query = ''' query VariantById($variantId: String!, $dataset: DatasetId!) { variant(variantId: $variantId, dataset: $dataset) { variant_id rsids exome { ac an af homozygote_count filters populations { id ac an } faf95 { popmax popmax_population } } genome { ac an af homozygote_count filters populations { id ac an } faf95 { popmax popmax_population } } } } ''' variant_id = f'{chrom}-{pos}-{ref}-{alt}' r = requests.post(GNOMAD_API, json={'query': query, 'variables': {'variantId': variant_id, 'dataset': dataset}}, timeout=30) r.raise_for_status() return r.json().get('data', {}).get('variant') def grpmax_faf95(payload): '''Extract the grpmax FAF95; the ACMG-grade frequency. Excludes bottleneck groups.''' exome = payload.get('exome') if payload else None if exome and exome.get('faf95'): return { 'faf95': exome['faf95'].get('popmax'), 'grpmax_ancestry': exome['faf95'].get('popmax_population'), 'source': 'exome' } genome = payload.get('genome') if payload else None if genome and genome.get('faf95'): return { 'faf95': genome['faf95'].get('popmax'), 'grpmax_ancestry': genome['faf95'].get('popmax_population'), 'source': 'genome' } return {'faf95': 0.0, 'grpmax_ancestry': None, 'source': 'absent'} ``` ## ACMG BS1/BA1 Application **Goal:** Apply Whiffin max-credible-AF framework to a candidate variant. **Approach:** Compute the gene-specific BS1 threshold from disease parameters, compare to `grpmax_faf95`. ```python def max_credible_af(prevalence, max_allelic_contribution=1.0, max_genetic_contribution=1.0, penetrance=1.0): '''Whiffin 2017 max-credible-AF formula. Args: prevalence: disease prevalence (e.g., 1/10000 = 1e-4) max_allelic_contribution: max contribution of single allele to disease in any case max_genetic_contribution: max contribution of this gene to disease in any case penetrance: probability that variant carriers develop disease Returns: max-credible per-allele frequency under dominant inheritance (use /2 for AR) ''' return (prevalence * max_genetic_contribution * max_allelic_contribution) / (penetrance * 2) def apply_bs1_ba1(grpmax_faf95_val, max_credible, ba1_threshold=0.05): '''Apply ClinGen SVI BS1/BA1 criteria. BA1 default 5% per ClinGen SVI; VCEP-specific overrides exist (Hearing Loss = 0.5%). BS1 = max-credible-AF specific to gene+disease. ''' if grpmax_faf95_val is None: return 'PM2_Supporting' # Absent or ultra-rare if grpmax_faf95_val > ba1_threshold: return 'BA1' if grpmax_faf95_val > max_credible: return 'BS1' return None # No criterion triggered; variant is consistent with rare-disease causation ``` ## Gene-Level Constraint (LOEUF) **Goal:** Retrieve gene constraint metrics with awareness of version mismatch for chrX/Y. **Approach:** Use v4 LOEUF for autosomes; fall back to v2.1.1 for chrX/Y. Report LOEUF decile, not raw value, to avoid cross-version comparison errors. ```python def query_gene_constraint(gene_symbol, dataset='gnomad_r4'): '''Pull gene constraint metrics. Note: v4 has no chrX/Y constraint; use v2 fallback.''' query = ''' query GeneById($symbol: String!) { gene(gene_symbol: $symbol, reference_genome: GRCh38) { gene_id symbol chrom gnomad_constraint { oe_lof oe_lof_lower oe_lof_upper oe_mis oe_mis_upper pli mis_z } } } ''' r = requests.post(GNOMAD_API,