1d0e2581c362c4971a20b464f5e593f23b01fba6 lrnassar Thu Aug 13 16:01:04 2026 -0700 Fix TP53 Provisional PM2 and BA1/BS1 frequency logic to match CSpec GN009. refs #37399 From the EvRepo cross-check. PM2_Supporting now applies to missense variants absent from gnomAD, not only rare-but-present ones: tp53AFfrequencies.py writes a per-assembly gnomAD present-set, and tp53ProvisionalClass.py applies PM2 when every nt path of a missense change is absent (a present-but-not-rare path blocks it, so variants seen at moderate frequency are not given PM2). BA1/BS1 now threshold on the non-founder ancestry-group max AF (AF_grpmax), per GN009's continental-subpopulation FAF rule, instead of the overall faf95, so variants common in a single ancestry group are no longer under-called. gnomAD v4.1 exposes only the single top group with no per-group filtering AF, so this uses raw AF_grpmax with a >=2000-allele floor as a documented approximation. Validated vs 137 EvRepo final classifications: PM2 disagreements 63->12, AF strength disagreements 4->2 (P47S/G360A/T312S now agree). diff --git src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py index e2437044d85..090f856c9b8 100644 --- src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py +++ src/hg/makeDb/scripts/tp53/tp53AFfrequencies.py @@ -1,98 +1,118 @@ #!/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: - BA1 FAF >= 0.001 stand-alone B - BS1 0.0003 <= FAF < 0.001 -4 pts + BA1 non-founder grpmax AF >= 0.001 stand-alone B + BS1 0.0003 <= non-founder grpmax AF < 0.001 -4 pts PM2_Supporting AF < 0.00003 global AND grpmax AF < 0.00004 +1 pt -Uses faf95 (col 16) from the UCSC gnomAD v4.1 bigBed, plus grpmax AF -(col 27) and global AF (col 15). CHIP note (col 29) surfaced in mouseover. +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 from -the per-ancestry check per CSpec — our PM2_Supporting uses grpmax as a -conservative proxy and flags when the proxy may miss qualifying variants. +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). """ 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" # 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 + 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 FAF >= 0.001 (0.1%)', - 'BS1': 'gnomAD v4.1 FAF in [0.0003, 0.001)', + '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)', '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 grpmax_af; "AF in the grpmax population" string grpmax_pop; "Population with grpmax AF" string chipNote; "gnomAD CHIP annotation (if any)" lstring _mouseOver; "HTML mouseover" ) """ -def classify(af_global, faf, af_grpmax): - if faf is not None and faf >= BA1_FAF: +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: return 'BA1' - if faf is not None and BS1_FAF_LOW <= faf < BS1_FAF_HIGH: + if grpmax_usable and BS1_FAF_LOW <= af_grpmax < 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': return None try: return float(s) except ValueError: @@ -131,60 +151,104 @@ 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 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] + 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) + code = classify(af_global, faf, af_grpmax, grpmax_pop, an_grpmax) 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, 'code': code, }) print(" classified: BA1={BA1} BS1={BS1} PM2={PM2} skipped={skipped}".format(**stats)) - return classified + 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'])) + else: + chain = "/cluster/data/hg38/bed/liftOver/hg38ToHg19.over.chain.gz" + in_bed = os.path.join(outdir, ".present_lift_in.bed") + out_bed = os.path.join(outdir, ".present_lift_out.bed") + unmapped = os.path.join(outdir, ".present_unmapped.bed") + with open(in_bed, 'w') as f: + for i, r in enumerate(all_records): + f.write("{}\t{}\t{}\t{}\n".format( + r['chrom'], r['hg38_start'], r['hg38_end'], i)) + lib.run_liftOver(in_bed, chain, out_bed, unmapped) + idxmap = {} + with open(out_bed) as f: + for line in f: + fl = line.rstrip("\n").split("\t") + if len(fl) >= 4: + idxmap[int(fl[3])] = (fl[0], int(fl[1])) + for i, r in enumerate(all_records): + if i in idxmap: + c, s = idxmap[i] + keys.append("{}-{}-{}-{}".format(c, s + 1, r['ref'], r['alt'])) + for p in (in_bed, out_bed, unmapped): + if os.path.exists(p): + os.remove(p) + with open(path, 'w') as f: + f.write("\n".join(keys) + "\n") + print(" wrote gnomAD present-set: {} keys -> {}".format(len(keys), path)) def emit_rows(classified, assembly, coord_lookup=None): """Format rows for the given assembly. coord_lookup is a dict keyed on hg38_name → (chrom, 0-based start, end) for hg19 liftOver. For hg38 pass 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 @@ -226,57 +290,59 @@ for line in f: flds = line.rstrip("\n").split("\t") if len(flds) >= 4: lookup[flds[3]] = (flds[0], int(flds[1]), int(flds[2])) for p in [input_bed, output_bed, unmapped]: if os.path.exists(p): os.remove(p) return lookup def build(db, outdir): print("=== {} ===".format(db)) os.makedirs(outdir, exist_ok=True) # We always query gnomAD on hg38 (the source), then lift to hg19 if needed tx_hg38 = lib.get_transcript_info('hg38') - classified = classify_and_build_rows(tx_hg38, tx_hg38['chrom']) + classified, all_records = classify_and_build_rows(tx_hg38, tx_hg38['chrom']) as_file = os.path.join(outdir, "TP53AF.as") lib.write_autosql(as_file, AUTOSQL) bed = os.path.join(outdir, "TP53AF_{}.bed".format(db)) bb = os.path.join(outdir, "TP53AF{}.bb".format(db.capitalize())) if db == 'hg38': lines = emit_rows(classified, 'hg38') with open(bed, 'w') as f: f.write("\n".join(lines) + "\n") lib.run_sort_bed(bed) lib.run_bedToBigBed(bed, as_file, bb, lib.chrom_sizes_path(db), "bed9+8") print(" wrote {}".format(bb)) + write_present_set(all_records, db, outdir) return # hg19 build: liftOver each record and rewrite display name lookup = liftover_hg38_to_hg19(classified, outdir) dropped = len(classified) - len(lookup) if dropped: print(" liftOver dropped {} variants".format(dropped)) lines = emit_rows(classified, 'hg19', coord_lookup=lookup) with open(bed, 'w') as f: f.write("\n".join(lines) + "\n") lib.run_sort_bed(bed) lib.run_bedToBigBed(bed, as_file, bb, lib.chrom_sizes_path(db), "bed9+8") print(" wrote {}".format(bb)) + write_present_set(all_records, db, outdir) return def main(): p = argparse.ArgumentParser(description=__doc__) p.add_argument('-o', '--output-dir', default=DEFAULT_OUTDIR) p.add_argument('--db', action='append', help='hg38 or hg19 (repeat). Default hg38.') args = p.parse_args() dbs = args.db if args.db else ['hg38'] for db in dbs: build(db, args.output_dir)