c05b0aa800df56a3fd1a5c5cbb903e39a54b6396 lrnassar Fri Aug 14 10:42:09 2026 -0700 TP53 Provisional: nt-aware splicing, caveat 4, and per-nt divergence flags. refs #37399 From the ClinGen TP53 VCEP review. The Provisional summary is one row per missense protein change, but splicing and allele frequency are nt-specific: - Splicing PP3 (SpliceAI >= 0.2, supersedes a missense BP4) now applies only when UNAMBIGUOUS: every nt path yielding the protein change is >= 0.2. When paths straddle the threshold the row keeps its missense code (fixes e.g. L35F, where c.105G>T=0.35 but c.105G>C=0.13 was wrongly flagged via the codon max). - Caveat (GN009): when splicing PP3 is applied, drop the protein-level functional assay code (PS3/BS3) from the point sum, since the assay tests the missense protein not the splicing effect. - New ntVerify column (bed9+10 -> 9+11, filterable): flags rows whose nt variants disagree on SpliceAI or on gnomAD AF, with a mouseover pointer to verify the exact nt change in the Bioinformatic / Allele Frequencies track. load_s2 now tracks spliceai_min; af_code_and_points returns af_divergent (AF aggregation unchanged). Reviewed clean by a fresh-eyes pass. diff --git src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py index e50911ca9c7..fe3c1f6c28d 100644 --- src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py +++ src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py @@ -13,31 +13,31 @@ * Splicing PP3 (SpliceAI >= 0.2) The point sum is bucketed into P / LP / VUS / LB / B per the CSpec classification ranges. BA1 is stand-alone Benign and forces class = Benign regardless of other evidence. DELIBERATELY EXCLUDED from the sum (documented in every mouseover): - PVS1 (null variants only; handled in separate track) - PS1 / PS2 / PS4 / PP1 / PP4 / BS4 (require clinical observations) - BP7 (computational, but synonymous/intronic only; out of scope for this missense-only track) This is NOT a ClinGen classification — the warning is in every mouseover since clinicians live in the mouseover, not the description page. -bigBed 9+10. +bigBed 9+11. """ import argparse 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" @@ -104,30 +104,31 @@ string name; "Missense change (e.g., R175H)" 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 provisionalClass; "Provisional class (Pathogenic/LP/VUS/LB/Benign)" int totalPoints; "Sum of applied Tavtigian points" string appliedCodes; "Codes contributing (semicolon-separated)" string pm1; "PM1 contribution" string ps3bs3; "PS3/BS3 contribution (from Table S3)" string pp3bp4; "PP3/BP4 contribution (Table S2 + splicing)" string af; "AF code (BA1/BS1/PM2_Supporting/none)" string bs2; "BS2 evidence (FLOSSIES cohort observation)" string spliceAI; "Max SpliceAI delta (from Table S2)" + string ntVerify; "nt-divergence flag: SpliceAI / gnomAD AF path-dependent, verify per-nt" lstring _mouseOver; "HTML mouseover" ) """ def parse_missense_short(p): if not isinstance(p, str): return None m = PROT_RE.match(p.strip()) if not m: return None return (m.group(1), int(m.group(2)), m.group(3)) def parse_hgvsp3(p): @@ -148,64 +149,68 @@ ws = wb["Supplementary Table S3"] out = {} for row in ws.iter_rows(min_row=4, values_only=True): k = parse_missense_short(row[0]) if not k: continue out[k] = { 'code': str(row[7]).strip() if row[7] is not None else 'No evidence', } return out def load_s2(path): """Table S2: per c.X>Y missense with preliminary PP3/BP4 + max SpliceAI. - Returns {(wt, codon, alt) -> {code, spliceai_max, paths: [(c_pos, ref_nt, alt_nt), ...]}}. + Returns {(wt, codon, alt) -> {code, spliceai (max), spliceai_min, + paths: [(c_pos, ref_nt, alt_nt), ...]}}. 'paths' enumerates every c.X>Y combination that yields the protein change — used downstream to look up AF at each contributing genomic coord. """ wb = openpyxl.load_workbook(path, data_only=True) ws = wb["Supplementary Table S2"] out = {} for row in ws.iter_rows(min_row=4, values_only=True): k = parse_hgvsp3(row[1]) if not k: continue hgvsc = row[0] path_tuple = None if isinstance(hgvsc, str): mc = HGVSC_RE.match(hgvsc.strip()) if mc: path_tuple = (int(mc.group(1)), mc.group(2), mc.group(3)) new_code = str(row[4]).strip() if row[4] is not None else 'No evidence' try: new_spliceai = float(row[5]) if row[5] is not None else 0.0 except (ValueError, TypeError): new_spliceai = 0.0 existing = out.get(k) if existing is None: out[k] = { 'code': new_code, - 'spliceai': new_spliceai, + 'spliceai': new_spliceai, # max across nt paths + 'spliceai_min': new_spliceai, # min across nt paths (straddle test) 'paths': [path_tuple] if path_tuple else [], } else: if _s2_priority(new_code) > _s2_priority(existing['code']): existing['code'] = new_code if new_spliceai > existing['spliceai']: existing['spliceai'] = new_spliceai + if new_spliceai < existing['spliceai_min']: + existing['spliceai_min'] = new_spliceai if path_tuple and path_tuple not in existing['paths']: existing['paths'].append(path_tuple) return out def _s2_priority(code): order = { 'PP3_moderate': 5, 'BP4_moderate': 5, 'PP3': 4, 'BP4': 4, 'No evidence': 1, } return order.get(code, 0) def load_hotspot_occurrences(path): @@ -338,66 +343,79 @@ return ('BP4_Moderate', -2) return ('', 0) 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. + + Also returns af_divergent: True when the nt paths yield different AF codes + (e.g. one path absent from gnomAD -> PM2 while another is present at moderate + frequency -> no code). The protein-level row can only show one, so these are + flagged for per-nt verification in the AF track. """ chrom = tx['chrom'] strand = tx['strand'] rcomp = {'A':'T','T':'A','C':'G','G':'C'} - seen = [] - any_present = False - any_absent = False + per_path = [] # per nt path: 'BA1' / 'BS1' / 'PM2_Supporting' / 'absent' / 'other' 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 + per_path.append(code) + elif key not in present_set: + per_path.append('absent') # absent from gnomAD -> PM2-eligible 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') + per_path.append('other') # present but not rare enough -> no code + coded = [c for c in per_path if c not in ('absent', 'other')] + any_absent = 'absent' in per_path + any_present = ('other' in per_path) or bool(coded) + + # Effective code per path (absent -> PM2, present-Other -> none); if the set + # disagrees the protein change is path-dependent for allele frequency. + eff = set('PM2_Supporting' if s == 'absent' else ('' if s == 'other' else s) + for s in per_path) + af_divergent = len(eff) > 1 + + if 'BA1' in coded: + return ('BA1', 0, 'stand-alone Benign', af_divergent) + if 'BS1' in coded: + return ('BS1', POINTS['BS1'], '-4 pts', af_divergent) + if 'PM2_Supporting' in coded: + return ('PM2_Supporting', POINTS['PM2_Supporting'], '+1 pt', af_divergent) if any_absent and not any_present: return ('PM2_Supporting', POINTS['PM2_Supporting'], - '+1 pt (absent from gnomAD)') - return ('', 0, '') + '+1 pt (absent from gnomAD)', af_divergent) + return ('', 0, '', af_divergent) 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: continue @@ -416,204 +434,256 @@ CAVEATS_STR = ( "PVS1 (null variants only — see PVS1 tracks); " "PS1/PS2/PS4/PP1/PP4/BS4 (require clinical observations); " "BP7 (synonymous/intronic; N/A to missense)." ) HEADER_WARNING = ( "<b style='color:#c00000'>" "NON-FINAL: NOT a ClinGen classification — preliminary point-sum only." "</b>" ) def mouseover(name, cls, pts, pm1, ps3, pp3, af_lbl, bs2_lbl, applied, - spliceai, splice_pp3_active, splice_overrode_bp4, codon, ba1): + spliceai, splice_pp3_active, splice_overrode_bp4, codon, ba1, + spliceai_min=0.0, splice_straddle=False, af_divergent=False, + ps3_suppressed=''): splice_section = "" - if spliceai and spliceai > 0: + if splice_straddle: + splice_section = ( + "<br><b>SpliceAI:</b> {mn:.2f}–{mx:.2f} — " + "<span style='color:#c00000'>path-dependent: the nt variants at this " + "codon straddle the 0.2 threshold, so splicing PP3 is not applied at " + "the protein level. Verify the exact nt change in the Bioinformatic " + "track.</span>" + ).format(mn=spliceai_min, mx=spliceai) + elif spliceai and spliceai > 0: if splice_pp3_active and splice_overrode_bp4: splice_section = ( "<br><b>SpliceAI:</b> {sa:.2f} — " "<span style='color:#c00000;font-weight:bold'>" "splicing PP3 applies (supersedes missense BP4); " "potential splice disruption." "</span>" ).format(sa=spliceai) elif splice_pp3_active: splice_section = ( "<br><b>SpliceAI:</b> {sa:.2f} — splicing PP3 applies" ).format(sa=spliceai) else: splice_section = "<br><b>SpliceAI:</b> {sa:.2f}".format(sa=spliceai) + ps3_suppress_section = "" + if ps3_suppressed: + ps3_suppress_section = ( + "<br><b>Note:</b> {c} functional evidence not applied — " + "splicing PP3 takes precedence per CSpec." + ).format(c=ps3_suppressed) + af_divergent_section = "" + if af_divergent: + af_divergent_section = ( + "<br><b>gnomAD AF (path-dependent):</b> " + "<span style='color:#c00000'>the nt variants at this codon differ in " + "allele frequency; verify the exact change in the Allele Frequencies " + "track.</span>" + ) codon72 = "" if codon == 72: codon72 = ( "<br><b>Note:</b> Codon 72 PS3/BS3 data is measured on the R72 " "haplotype (rs1042522); reference is P72. P72X variants are NOT " "in Table S3 (only R72X is). See description page for details." ) ba1_section = "" if ba1: ba1_section = ( "<br><b style='color:#0a0'>BA1 stand-alone Benign:</b> " "FAF ≥ 0.001 in gnomAD v4.1 forces class = Benign." ) return ( "{warn}" "<br><b>Variant:</b> p.{name} (NP_000537.3)" "<br><b>Provisional class:</b> {cls} ({pts} pts)" "<br><b>PM1:</b> {pm1}" "<br><b>PS3/BS3:</b> {ps3}" "<br><b>PP3/BP4:</b> {pp3}" "<br><b>AF (gnomAD v4.1):</b> {af}" "<br><b>BS2 (FLOSSIES):</b> {bs2}" "<br><b>Applied codes:</b> {applied}" - "{splice}{ba1_section}{codon72}" + "{splice}{ps3sup}{afdiv}{ba1_section}{codon72}" "<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, + ps3sup=ps3_suppress_section, + afdiv=af_divergent_section, ba1_section=ba1_section, codon72=codon72, cav=CAVEATS_STR) 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 + spliceai_min = 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) + spliceai_min = s2_rec.get('spliceai_min', spliceai) 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 + # Splicing rule (SpliceAI >= 0.2 -> PP3 splicing, superseding a missense + # BP4), applied only when UNAMBIGUOUS: every nt path yielding this protein + # change is >= threshold. SpliceAI is nt-specific; when the paths straddle + # the threshold this protein-level row cannot represent it, so we keep the + # missense code and flag the row for per-nt verification (Bioinformatic). + splice_all = bool(paths) and spliceai_min >= SPLICE_PP3_THRESHOLD + splice_straddle = (spliceai >= SPLICE_PP3_THRESHOLD + and spliceai_min < SPLICE_PP3_THRESHOLD) + splice_pp3_active = splice_all 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 + # CSpec caveat: do not apply a protein-level functional assay (PS3/BS3) + # when the variant is classified as splicing via PP3. + ps3_suppressed = '' + if splice_pp3_active and ps3_lbl: + ps3_suppressed = ps3_lbl + ps3_lbl, ps3_pts = ('', 0) + # Allele-frequency code (BA1 / BS1 / PM2_Supporting) - af_code, af_pts, af_qty = af_code_and_points( + af_code, af_pts, af_qty, af_divergent = af_code_and_points( 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) color = CLASS_COLOR[cls] + # nt-divergence flags: this protein-level row summarizes >1 nt variant + # whose SpliceAI or gnomAD AF disagree; verify the exact nt change in the + # Bioinformatic / Allele Frequency tracks. + nt_flags = [] + if splice_straddle: + nt_flags.append('SpliceAI') + if af_divergent: + nt_flags.append('gnomAD AF') + nt_verify = " + ".join(nt_flags) if nt_flags else "-" + applied_codes = [] if pm1_lbl: applied_codes.append("{} (+{})".format(pm1_lbl, pm1_pts)) if ps3_lbl: sign = "+" if ps3_pts > 0 else "" applied_codes.append("{} ({}{})".format(ps3_lbl, sign, ps3_pts)) if pp3_lbl: sign = "+" if pp3_pts > 0 else "" applied_codes.append("{} ({}{})".format(pp3_lbl, sign, pp3_pts)) if af_code == 'BA1': applied_codes.append("BA1 (stand-alone B)") elif af_code: sign = "+" if af_pts > 0 else "" applied_codes.append("{} ({}{})".format(af_code, sign, af_pts)) if bs2_applies: applied_codes.append("{} ({})".format(bs2_label, bs2_pts)) applied = "; ".join(applied_codes) short = "{}{}{}".format(wt, codon, alt) mo = mouseover(short, cls, total, pm1_lbl, ps3_lbl, pp3_lbl, af_lbl, bs2_lbl, applied, spliceai, splice_pp3_active, splice_overrode_bp4, - codon, ba1) + codon, ba1, spliceai_min, splice_straddle, af_divergent, + ps3_suppressed) segs = lib.aa_codon_genomic(codon, tx) for g_start, g_end, _ex in segs: lines.append("\t".join([ chrom, str(g_start), str(g_end), short, "0", ".", str(g_start), str(g_end), color, cls, str(total), applied, pm1_lbl or "-", ps3_lbl or "-", pp3_lbl or "-", af_code or "-", bs2_label if bs2_applies else "-", "{:.2f}".format(spliceai) if spliceai else "0.00", + nt_verify, 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: {} gnomAD present: {} FLOSSIES BS2: {}".format( len(s3), len(s2), len(hotspots), len(af_lookup), len(present_set), len(flossies_lookup))) tx = lib.get_transcript_info(db) 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") + lib.run_bedToBigBed(bed, as_file, bb, lib.chrom_sizes_path(db), "bed9+11") print(" wrote {}".format(bb)) from collections import Counter cnt = Counter() af_cnt = Counter() bs2_cnt = Counter() with open(bed) as f: for line in f: flds = line.split("\t") cnt[flds[9]] += 1 af_cnt[flds[15]] += 1 if flds[16] in BS2_TIER_POINTS: bs2_cnt[flds[16]] += 1 print(" Class distribution:") for k, n in cnt.most_common():