e28d9673130c2d129d7f7684897fb809b5a17325 lrnassar Thu Aug 13 16:38:07 2026 -0700 Read TP53 BA1/BS1 from source gnomAD VCF for exact per-ancestry faf95. refs #37399 Replace the raw AF_grpmax approximation with the CSpec-exact metric: the max faf95 across non-founder continental ancestry groups (afr/amr/eas/nfe/sas), read from the source gnomAD v4.1 sites VCF via tabix. The /gbdb bigBed only carries the overall faf95 and the raw grpmax AF, not per-group faf95. faf95 is a Poisson 95% CI lower bound, so it discounts small groups (subsuming the CSpec >=2000 allele rule) and removes the raw-AF over-calls. Founder groups (asj/fin/mid/remaining) are excluded by omission from the max. Effect vs the AF_grpmax approximation: AF-track BS1 75->53; Provisional over-calls eliminated (Y107H BA1->BS1, now matches the EvRepo final call, and no we-only BA1/BS1 remain); T312S drops to no AF code because its exact faf95 0.000299 is just under the 0.0003 BS1 threshold. Also repurposes the AF bed's chipNote column to fafGroup (the ancestry group giving the faf95). diff --git src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py index 090f856c9b8..96b8d49dbc3 100644 --- src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py +++ src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py @@ -1,213 +1,238 @@ #!/usr/bin/env python3 """ TP53 VCEP Allele Frequencies (BA1/BS1/PM2) track generator. -Pulls gnomAD v4.1 exome variants at the TP53 locus and classifies them -per CSpec GN009 v2.4.0 thresholds: +Reads gnomAD v4.1 exome variants at the TP53 locus from the source sites VCF +(via tabix) and classifies them per CSpec GN009 v2.4.0 thresholds: - BA1 non-founder grpmax AF >= 0.001 stand-alone B - BS1 0.0003 <= non-founder grpmax AF < 0.001 -4 pts + BA1 non-founder ancestry-group faf95 >= 0.001 stand-alone B + BS1 0.0003 <= non-founder ancestry-group faf95 < 0.001 -4 pts PM2_Supporting AF < 0.00003 global AND grpmax AF < 0.00004 +1 pt -BA1/BS1 use the max AF of a single non-founder continental ancestry group -(AF_grpmax, col 27), per the CSpec's continental-subpopulation FAF rule, when -that group has >=2000 alleles tested (AN_grpmax, col 26). Global AF (col 15) -and faf95 (col 16) are shown in the mouseover. CHIP note (col 29) surfaced too. - -Founder-effect ancestry groups (AJ/FIN/AMI/MID/Remaining) are EXCLUDED per -CSpec: if the top group is a founder group we do not apply BA1/BS1, and because -gnomAD v4.1 exposes only the single top group we cannot recover the next-highest -non-founder group in that case (a documented limitation). +BA1/BS1 use the CSpec's continental-subpopulation filtering allele frequency: +the maximum faf95 across the non-founder ancestry groups (afr/amr/eas/nfe/sas). +The /gbdb gnomAD bigBed only carries the overall faf95 and the raw grpmax AF, so +we read the source VCF, which has per-group faf95 (faf95_afr, faf95_amr, ...). +faf95 is a Poisson 95% CI lower bound, so it already discounts small groups (the +CSpec's >=2000-allele requirement is subsumed). Founder-effect groups +(AJ/FIN/MID/Remaining) are excluded by simply not including them in the max. +PM2_Supporting uses global AF plus AF_grpmax as a per-ancestry proxy. """ import argparse import os import subprocess import sys sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import tp53FuncLib as lib DEFAULT_OUTDIR = "/hive/users/lrnassar/claude/RM37399/afFrequencies" -GNOMAD_BB = "/gbdb/hg38/gnomAD/v4.1/exomes/exomes.bb" +# Source gnomAD v4.1 exomes sites VCF (tabix-indexed, chr17). Read directly +# because the /gbdb bigBed does not carry per-ancestry faf95. +GNOMAD_VCF = "/hive/data/outside/gnomAD.4/v4.1/exomes/gnomad.exomes.v4.1.sites.chr17.vcf.bgz" +TABIX = "/cluster/bin/x86_64/tabix" # CSpec GN009 v2.4.0 TP53 thresholds BA1_FAF = 0.001 BS1_FAF_LOW = 0.0003 BS1_FAF_HIGH = 0.001 PM2_AF_GLOBAL_MAX = 0.00003 PM2_AF_GRPMAX_MAX = 0.00004 -# CSpec GN009 excludes these founder-effect ancestry groups from the -# frequency-based codes; BA1/BS1 use the max AF of a *non-founder* continental -# ancestry group. gnomAD v4.1 exposes only the single top group (AF_grpmax) with -# no per-group filtering AF, so we threshold on AF_grpmax when that group is -# non-founder and >=2000 alleles were tested (a documented approximation of the -# continental-subpopulation FAF the CSpec specifies). -FOUNDER_GROUPS = {'Ashkenazi Jewish', 'Finnish', 'Amish', - 'Middle Eastern', 'Remaining'} -GRPMAX_MIN_AN = 2000 +# Non-founder continental ancestry groups whose faf95 counts toward BA1/BS1 +# (gnomAD abbrev -> display name). The founder-effect groups the CSpec excludes +# (asj/fin/mid/remaining) are simply left out. +NON_FOUNDER_FAF = { + 'afr': 'African/African American', + 'amr': 'Admixed American', + 'eas': 'East Asian', + 'nfe': 'European (Non-Finnish)', + 'sas': 'South Asian', +} +# For displaying the AF_grpmax population (may be a founder group). +GRPMAX_POP_NAMES = { + 'afr': 'African/African American', 'amr': 'Admixed American', + 'asj': 'Ashkenazi Jewish', 'eas': 'East Asian', 'fin': 'Finnish', + 'mid': 'Middle Eastern', 'nfe': 'European (Non-Finnish)', + 'sas': 'South Asian', 'remaining': 'Remaining', +} COLORS = { 'BA1': '2,82,66', # dark teal (stand-alone B) 'BS1': '35,159,134', # teal 'PM2_Supporting': '138,111,158', # purple } POINTS = { 'BA1': 'stand-alone B', 'BS1': '-4 pts', 'PM2_Supporting': '+1 pt', } RULES = { - 'BA1': 'gnomAD v4.1 non-founder ancestry-group AF >= 0.001 (0.1%)', - 'BS1': 'gnomAD v4.1 non-founder ancestry-group AF in [0.0003, 0.001)', + 'BA1': 'gnomAD v4.1 non-founder ancestry-group faf95 >= 0.001 (0.1%)', + 'BS1': 'gnomAD v4.1 non-founder ancestry-group faf95 in [0.0003, 0.001)', 'PM2_Supporting': 'gnomAD v4.1 AF < 3e-5 global AND grpmax < 4e-5', } AUTOSQL = """table TP53AF "TP53 VCEP ACMG allele frequency classifications from gnomAD v4.1 exomes" ( string chrom; "Reference sequence chromosome or scaffold" uint chromStart; "Start position in chromosome" uint chromEnd; "End position in chromosome" string name; "Variant display name" uint score; "Not used, all 0" char[1] strand; "Not used, all ." uint thickStart; "Same as chromStart" uint thickEnd; "Same as chromEnd" uint reserved; "RGB color" string acmgCode; "BA1 / BS1 / PM2_Supporting" string points; "Tavtigian points" string af; "Global AF in gnomAD v4.1 exomes" - string faf; "Filter allele frequency (faf95)" + string faf; "Max faf95 across non-founder ancestry groups (BA1/BS1 metric)" string grpmax_af; "AF in the grpmax population" string grpmax_pop; "Population with grpmax AF" - string chipNote; "gnomAD CHIP annotation (if any)" + string fafGroup; "Non-founder ancestry group with the max faf95" lstring _mouseOver; "HTML mouseover" ) """ -def classify(af_global, faf, af_grpmax, grpmax_pop, an_grpmax): - # BA1/BS1 per CSpec GN009 use the filtering AF of a single non-founder - # continental ancestry group. Approximated here by AF_grpmax when the top - # group is non-founder and >=2000 alleles were tested. - grpmax_usable = (af_grpmax is not None - and grpmax_pop is not None - and grpmax_pop not in FOUNDER_GROUPS - and an_grpmax is not None and an_grpmax >= GRPMAX_MIN_AN) - if grpmax_usable and af_grpmax >= BA1_FAF: +def classify(af_global, af_grpmax, faf_ba1bs1): + # BA1/BS1 per CSpec GN009 use the filtering AF (faf95) of a single + # non-founder continental ancestry group; faf_ba1bs1 is the max such value. + if faf_ba1bs1 is not None and faf_ba1bs1 >= BA1_FAF: return 'BA1' - if grpmax_usable and BS1_FAF_LOW <= af_grpmax < BS1_FAF_HIGH: + if faf_ba1bs1 is not None and BS1_FAF_LOW <= faf_ba1bs1 < BS1_FAF_HIGH: return 'BS1' # PM2_Supporting: rare globally AND grpmax under limit if (af_global is not None and af_global < PM2_AF_GLOBAL_MAX and af_grpmax is not None and af_grpmax < PM2_AF_GRPMAX_MAX and (af_global > 0 or af_grpmax > 0)): return 'PM2_Supporting' return None def safe_float(s): - if s is None or s == '' or s == 'N/A': + if s is None or s in ('', 'N/A', '.'): return None try: return float(s) except ValueError: return None -def mouseover(display, code, ref, alt, af_global, faf, af_grpmax, grpmax_pop, chip): - chip_line = "" - if chip and chip not in ('N/A', ''): - chip_line = "
CHIP note: {}".format(chip) +def parse_info(info): + """Parse a VCF INFO column into a {key: value} dict (flag keys omitted).""" + d = {} + for kv in info.split(';'): + if '=' in kv: + k, v = kv.split('=', 1) + d[k] = v + return d + + +def nonfounder_max_faf(info): + """Return (max faf95 across non-founder continental groups, group display + name) from a parsed INFO dict, per the CSpec's continental-subpopulation + FAF rule for BA1/BS1. Returns (None, None) if no group has a faf95.""" + best_v = None + best_g = None + for abbr, name in NON_FOUNDER_FAF.items(): + v = safe_float(info.get('faf95_' + abbr)) + if v is not None and (best_v is None or v > best_v): + best_v = v + best_g = name + return best_v, best_g + + +def mouseover(display, code, ref, alt, af_global, faf, faf_group, af_grpmax, grpmax_pop): + faf_txt = "{:.2e}".format(faf) if faf is not None else "N/A" + if faf is not None and faf_group: + faf_txt = "{} ({})".format(faf_txt, faf_group) return ( "Variant: {disp} ({ref}>{alt})" "
ACMG code: {code} ({pts})" "
Rule: {rule}" "
Global AF: {af}" - "
FAF (faf95): {faf}" + "
Filtering AF (max non-founder group faf95): {faf}" "
grpmax AF: {gmax} ({gpop})" - "{chip}" "
Source: gnomAD v4.1 exomes" ).format( disp=display, ref=ref, alt=alt, code=code, pts=POINTS[code], rule=RULES[code], af="{:.2e}".format(af_global) if af_global is not None else "N/A", - faf="{:.2e}".format(faf) if faf is not None else "N/A", + faf=faf_txt, gmax="{:.2e}".format(af_grpmax) if af_grpmax is not None else "N/A", gpop=grpmax_pop or "N/A", - chip=chip_line, ) def classify_and_build_rows(tx, chrom): - """Query gnomAD v4.1 exomes on hg38 and emit a list of classified rows - keyed by an immutable hg38 identifier. The hg38 identifier is used as - the 'name' field so the hg19 build can look up the same row after liftOver - and rewrite the display text to reflect hg19 coords.""" - raw_bed = "/tmp/tp53_gnomad_{}.bed".format(os.getpid()) - subprocess.run( - ["bigBedToBed", GNOMAD_BB, - "-chrom=" + chrom, - "-start=" + str(tx['txStart']), - "-end=" + str(tx['txEnd']), - raw_bed], - check=True) - with open(raw_bed) as f: - rows = [line.rstrip("\n").split("\t") for line in f] - os.remove(raw_bed) - print(" {} variants in TP53 region (hg38)".format(len(rows))) - - # Build rows keyed on the hg38 display id + """Read the source gnomAD v4.1 exomes VCF on hg38 via tabix and emit a list + of classified rows keyed by an immutable hg38 identifier. The hg38 id is used + as the 'name' field so the hg19 build can look up the same row after liftOver + and rewrite the display text to reflect hg19 coords. Reading the source VCF + (not the /gbdb bigBed) gives the per-ancestry faf95 that BA1/BS1 require.""" + region = "{}:{}-{}".format(chrom, tx['txStart'] + 1, tx['txEnd']) + out = subprocess.run([TABIX, GNOMAD_VCF, region], + capture_output=True, text=True, check=True).stdout + vcf_lines = [ln for ln in out.splitlines() if ln and not ln.startswith('#')] + print(" {} variants in TP53 region (hg38)".format(len(vcf_lines))) + classified = [] # list of dicts with all fields; hg38 coords fixed all_records = [] # every gnomAD variant, for the PM2 "present in gnomAD" set - stats = dict(total=len(rows), BA1=0, BS1=0, PM2=0, skipped=0) - for r in rows: - c_start = int(r[1]) - c_end = int(r[2]) - ref = r[9] - alt = r[10] + stats = dict(total=len(vcf_lines), BA1=0, BS1=0, PM2=0, skipped=0, multi=0) + for ln in vcf_lines: + f = ln.split('\t') + pos = int(f[1]) # 1-based VCF POS + ref = f[3] + alt = f[4] + if ',' in alt: # multi-allelic (none expected in v4.1 sites VCF) + stats['multi'] += 1 + continue + info = parse_info(f[7]) + c_start = pos - 1 + c_end = c_start + len(ref) all_records.append({'chrom': chrom, 'hg38_start': c_start, 'hg38_end': c_end, 'ref': ref, 'alt': alt}) - af_global = safe_float(r[14]) - faf = safe_float(r[15]) - af_grpmax = safe_float(r[26]) - an_grpmax = safe_float(r[25]) - grpmax_pop = r[23] if r[23] not in ('N/A', '') else None - chip = r[28] if len(r) > 28 else '' - hg38_name = "{}-{}-{}-{}".format(chrom, c_start + 1, ref, alt) - - code = classify(af_global, faf, af_grpmax, grpmax_pop, an_grpmax) + af_global = safe_float(info.get('AF')) + af_grpmax = safe_float(info.get('AF_grpmax')) + grpmax_pop = GRPMAX_POP_NAMES.get(info.get('grpmax'), info.get('grpmax')) + faf_ba1bs1, faf_group = nonfounder_max_faf(info) + hg38_name = "{}-{}-{}-{}".format(chrom, pos, ref, alt) + + code = classify(af_global, af_grpmax, faf_ba1bs1) if code is None: stats['skipped'] += 1 continue stats[code if code in ('BA1', 'BS1') else 'PM2'] += 1 classified.append({ 'hg38_name': hg38_name, 'hg38_start': c_start, 'hg38_end': c_end, 'chrom': chrom, 'ref': ref, 'alt': alt, - 'af_global': af_global, 'faf': faf, 'af_grpmax': af_grpmax, - 'grpmax_pop': grpmax_pop, 'chip': chip, + 'af_global': af_global, 'faf': faf_ba1bs1, 'faf_group': faf_group, + 'af_grpmax': af_grpmax, 'grpmax_pop': grpmax_pop, 'code': code, }) - print(" classified: BA1={BA1} BS1={BS1} PM2={PM2} skipped={skipped}".format(**stats)) + print(" classified: BA1={BA1} BS1={BS1} PM2={PM2} skipped={skipped} " + "multiallelic_skipped={multi}".format(**stats)) return classified, all_records def write_present_set(all_records, db, outdir): """Write the set of every gnomAD variant key ('chrom-pos1-ref-alt') present at the TP53 locus, in the coordinates of the requested assembly. The Provisional track uses this to tell 'absent from gnomAD' (PM2_Supporting applies) apart from 'present but not rare enough to be coded' (no PM2). For hg19 the hg38 coords are lifted so the keys match the hg19 Provisional build.""" path = os.path.join(outdir, "TP53AF_present_{}.txt".format(db)) keys = [] if db == 'hg38': for r in all_records: keys.append("{}-{}-{}-{}".format( r['chrom'], r['hg38_start'] + 1, r['ref'], r['alt'])) @@ -245,42 +270,42 @@ None — rows use their native hg38 coords.""" lines = [] for rec in classified: if assembly == 'hg38': chrom = rec['chrom'] start = rec['hg38_start'] end = rec['hg38_end'] else: if rec['hg38_name'] not in coord_lookup: continue chrom, start, end = coord_lookup[rec['hg38_name']] # Use assembly-appropriate display name so hg19 viewers see hg19 pos display = "{}-{}-{}-{}".format(chrom, start + 1, rec['ref'], rec['alt']) color = COLORS[rec['code']] mo = mouseover(display, rec['code'], rec['ref'], rec['alt'], - rec['af_global'], rec['faf'], rec['af_grpmax'], - rec['grpmax_pop'], rec['chip']) + rec['af_global'], rec['faf'], rec['faf_group'], + rec['af_grpmax'], rec['grpmax_pop']) lines.append("\t".join([ chrom, str(start), str(end), display, "0", ".", str(start), str(end), color, rec['code'], POINTS[rec['code']], "{:.2e}".format(rec['af_global']) if rec['af_global'] is not None else "N/A", "{:.2e}".format(rec['faf']) if rec['faf'] is not None else "N/A", "{:.2e}".format(rec['af_grpmax']) if rec['af_grpmax'] is not None else "N/A", rec['grpmax_pop'] or "N/A", - rec['chip'] or "", + rec['faf_group'] or "N/A", mo, ])) return lines def liftover_hg38_to_hg19(classified, outdir): """Lift each hg38 coord to hg19, returning dict hg38_name → (chrom,start,end).""" chain = "/cluster/data/hg38/bed/liftOver/hg38ToHg19.over.chain.gz" input_bed = os.path.join(outdir, ".tp53af_lift_in.bed") output_bed = os.path.join(outdir, ".tp53af_lift_out.bed") unmapped = os.path.join(outdir, ".tp53af_unmapped.bed") with open(input_bed, 'w') as f: for rec in classified: f.write("{}\t{}\t{}\t{}\n".format( rec['chrom'], rec['hg38_start'], rec['hg38_end'], rec['hg38_name']))