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/tp53ProvisionalClass.py src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py index 82d3ed6c5b2..e50911ca9c7 100644 --- src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py +++ src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py @@ -32,30 +32,31 @@ import json import os import re import sys import openpyxl sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import tp53FuncLib as lib DEFAULT_OUTDIR = "/hive/users/lrnassar/claude/RM37399/provisionalClass" SRC_S3 = "/hive/users/lrnassar/claude/RM37399/tp53_downloads/Functional-worksheet.xlsx" SRC_S2 = "/hive/users/lrnassar/claude/RM37399/tp53_downloads/bioinformatic_worksheet.xlsx" HOTSPOTS_JSON = "/hive/users/lrnassar/claude/RM37399/cancerHotspots/cancerhotspots_single.json" AF_BED_TPL = "/hive/users/lrnassar/claude/RM37399/afFrequencies/TP53AF_{}.bed" +AF_PRESENT_TPL = "/hive/users/lrnassar/claude/RM37399/afFrequencies/TP53AF_present_{}.txt" FLOSSIES_BED_TPL = "/hive/users/lrnassar/claude/RM37399/flossies/TP53Flossies_{}.bed" PM1_HARDCODED_CODONS = {175, 245, 248, 249, 273, 282} PROT_RE = re.compile(r'^([A-Z])(\d+)([A-Z])$') HGVSP3_RE = re.compile(r'^p\.([A-Z][a-z]{2})(\d+)([A-Z][a-z]{2})$') HGVSC_RE = re.compile(r'^c\.(\d+)([ACGT])>([ACGT])$') THREE_TO_ONE = { 'Ala':'A','Arg':'R','Asn':'N','Asp':'D','Cys':'C','Glu':'E','Gln':'Q', 'Gly':'G','His':'H','Ile':'I','Leu':'L','Lys':'K','Met':'M','Phe':'F', 'Pro':'P','Ser':'S','Thr':'T','Trp':'W','Tyr':'Y','Val':'V', } # Tavtigian points per code @@ -242,30 +243,48 @@ return out with open(path) as f: for line in f: flds = line.rstrip("\n").split("\t") if len(flds) < 11: continue disp = flds[3] code = flds[9] m = re.match(r'^(chr[^-]+)-(\d+)-([ACGT-]+)-([ACGT-]+)$', disp) if not m: continue out[(m.group(1), int(m.group(2)), m.group(3), m.group(4))] = code return out +def load_gnomad_present(db): + """Read TP53AF_present_<db>.txt (one 'chrN-pos-ref-alt' key per line, 1-based) + and return a set of (chrom, pos1based, ref, alt) for every gnomAD variant at + the TP53 locus. Used to tell 'absent from gnomAD' from 'present but not + coded' when deciding PM2_Supporting.""" + out = set() + path = AF_PRESENT_TPL.format(db) + if not os.path.exists(path): + return out + with open(path) as f: + for line in f: + line = line.strip() + m = re.match(r'^(chr[^-]+)-(\d+)-([ACGT-]+)-([ACGT-]+)$', line) + if m: + out.add((m.group(1), int(m.group(2)), m.group(3), m.group(4))) + return out + + BS2_TIER_POINTS = {'BS2': -4, 'BS2_Moderate': -2, 'BS2_Supporting': -1} def load_flossies_lookup(db): """Read TP53Flossies_<db>.bed and return {(chrom, pos1based, ref, alt) -> (bs2_label, points)} for rows that meet BS2 at any tier (BS2 / BS2_Moderate / BS2_Supporting). Match is variant-specific (a FLOSSIES observation supports BS2 only for the exact nt change observed, not for any change at the same codon). The tier reflects the carrier count per CSpec GN009.""" out = {} path = FLOSSIES_BED_TPL.format(db) if not os.path.exists(path): return out with open(path) as f: for line in f: @@ -308,62 +327,76 @@ return ('', 0) def pp3_bp4_label_and_points(code): if code == 'PP3_moderate': return ('PP3_Moderate', 2) if code == 'PP3': return ('PP3', 1) if code == 'BP4': return ('BP4', -1) if code == 'BP4_moderate': return ('BP4_Moderate', -2) return ('', 0) -def af_code_and_points(wt, codon, alt, paths, af_lookup, tx): +def af_code_and_points(wt, codon, alt, paths, af_lookup, present_set, tx): """For a (wt, codon, alt) protein change, look up gnomAD v4.1 AF at every c.X>Y path and return the strongest applicable code. Priority (most-benign first): BA1 > BS1 > PM2_Supporting > none. BA1 is stand-alone Benign; the caller forces classification = Benign. + + PM2_Supporting applies either when a path is coded PM2 (present but rare) or + when every path is absent from gnomAD (present_set). A path that is present + in gnomAD but not rare enough to be coded blocks the absent-based PM2, so a + variant seen in gnomAD at moderate frequency is not given PM2. """ chrom = tx['chrom'] strand = tx['strand'] rcomp = {'A':'T','T':'A','C':'G','G':'C'} seen = [] + any_present = False + any_absent = False for c_pos, c_ref, c_alt in paths: g = lib.cdna_coding_to_genomic(c_pos, tx) if g is None: continue if strand == '-': g_ref = rcomp.get(c_ref, c_ref) g_alt = rcomp.get(c_alt, c_alt) else: g_ref = c_ref g_alt = c_alt # AF lookup uses 1-based start key = (chrom, g + 1, g_ref, g_alt) code = af_lookup.get(key) if code: seen.append(code) + if key in present_set: + any_present = True + else: + any_absent = True if 'BA1' in seen: return ('BA1', 0, 'stand-alone Benign') if 'BS1' in seen: return ('BS1', POINTS['BS1'], '-4 pts') if 'PM2_Supporting' in seen: return ('PM2_Supporting', POINTS['PM2_Supporting'], '+1 pt') + if any_absent and not any_present: + return ('PM2_Supporting', POINTS['PM2_Supporting'], + '+1 pt (absent from gnomAD)') return ('', 0, '') def bs2_evidence(wt, codon, paths, flossies_lookup, tx): """Return (bs2_label, points) for the strongest FLOSSIES BS2 tier matching any (c.X>Y) path by exact genomic position AND ref/alt, else (None, 0). BS2 requires the SAME nt change — an observation of c.1120G>A does not support BS2 for c.1120G>C even though both yield p.G374R.""" chrom = tx['chrom'] strand = tx['strand'] rcomp = {'A':'T','T':'A','C':'G','G':'C'} best = (None, 0) for c_pos, c_ref, c_alt in paths: g = lib.cdna_coding_to_genomic(c_pos, tx) if g is None: @@ -439,62 +472,62 @@ "<br><b>NOT included in this sum:</b> {cav}" ).format(warn=HEADER_WARNING, name=name, cls=cls, pts=pts, pm1=pm1 or "No contribution", ps3=ps3 or "No contribution", pp3=pp3 or "No contribution", af=af_lbl or "No contribution", bs2=bs2_lbl, applied=applied or "(none)", splice=splice_section, ba1_section=ba1_section, codon72=codon72, cav=CAVEATS_STR) -def generate_bed(s3, s2, hotspot_occ, af_lookup, flossies_lookup, tx): +def generate_bed(s3, s2, hotspot_occ, af_lookup, present_set, flossies_lookup, tx): lines = [] chrom = tx['chrom'] for (wt, codon, alt), s3_rec in sorted(s3.items(), key=lambda kv: (kv[0][1], kv[0][2])): pm1_lbl, pm1_pts = pm1_code_and_points(wt, codon, alt, hotspot_occ) ps3_lbl, ps3_pts = ps3_bs3_label_and_points(s3_rec['code']) s2_rec = s2.get((wt, codon, alt)) pp3_lbl, pp3_pts = ('', 0) spliceai = 0.0 paths = [] if s2_rec: pp3_lbl, pp3_pts = pp3_bp4_label_and_points(s2_rec['code']) spliceai = s2_rec.get('spliceai', 0.0) paths = s2_rec.get('paths', []) # Megan's splicing rule: SpliceAI >= 0.2 -> PP3 splicing applies. # When the missense call is BP4 / BP4_Moderate, splicing PP3 # supersedes — replace the BP4 contribution with PP3 (+1). splice_pp3_active = spliceai >= SPLICE_PP3_THRESHOLD splice_overrode_bp4 = False if splice_pp3_active: if pp3_lbl in ('BP4', 'BP4_Moderate'): pp3_lbl = 'PP3 (splicing)' pp3_pts = 1 splice_overrode_bp4 = True elif not pp3_lbl: pp3_lbl = 'PP3 (splicing)' pp3_pts = 1 # Allele-frequency code (BA1 / BS1 / PM2_Supporting) af_code, af_pts, af_qty = af_code_and_points( - wt, codon, alt, paths, af_lookup, tx) + wt, codon, alt, paths, af_lookup, present_set, tx) af_lbl = "{} ({})".format(af_code, af_qty) if af_code else '' ba1 = (af_code == 'BA1') # BS2 from FLOSSIES (tiered by carrier count per CSpec GN009) bs2_label, bs2_pts = bs2_evidence(wt, codon, paths, flossies_lookup, tx) bs2_applies = bs2_label is not None bs2_lbl = "{} ({} pts)".format(bs2_label, bs2_pts) if bs2_applies else "Not observed" total = pm1_pts + ps3_pts + pp3_pts + af_pts + bs2_pts if ba1: cls = 'Benign' else: cls = bucket(total) @@ -538,38 +571,39 @@ af_code or "-", bs2_label if bs2_applies else "-", "{:.2f}".format(spliceai) if spliceai else "0.00", mo, ])) return lines def build(db, outdir): print("=== {} ===".format(db)) os.makedirs(outdir, exist_ok=True) s3 = load_s3(SRC_S3) s2 = load_s2(SRC_S2) hotspots = load_hotspot_occurrences(HOTSPOTS_JSON) af_lookup = load_af_lookup(db) + present_set = load_gnomad_present(db) flossies_lookup = load_flossies_lookup(db) print(" S3 entries: {} S2 entries: {} " - "cancerhotspots: {} AF: {} FLOSSIES BS2: {}".format( + "cancerhotspots: {} AF: {} gnomAD present: {} FLOSSIES BS2: {}".format( len(s3), len(s2), len(hotspots), - len(af_lookup), len(flossies_lookup))) + len(af_lookup), len(present_set), len(flossies_lookup))) tx = lib.get_transcript_info(db) - bed_lines = generate_bed(s3, s2, hotspots, af_lookup, flossies_lookup, tx) + bed_lines = generate_bed(s3, s2, hotspots, af_lookup, present_set, flossies_lookup, tx) print(" {} BED rows".format(len(bed_lines))) as_file = os.path.join(outdir, "TP53ProvisionalClass.as") lib.write_autosql(as_file, AUTOSQL) bed = os.path.join(outdir, "TP53ProvisionalClass_{}.bed".format(db)) with open(bed, 'w') as f: f.write("\n".join(bed_lines) + "\n") lib.run_sort_bed(bed) bb = os.path.join(outdir, "TP53ProvisionalClass{}.bb".format(db.capitalize())) lib.run_bedToBigBed(bed, as_file, bb, lib.chrom_sizes_path(db), "bed9+10") print(" wrote {}".format(bb)) from collections import Counter cnt = Counter() af_cnt = Counter()