7fb0ec3d9354f0bb4b059b0106ecaf5c40889a05 lrnassar Thu Aug 13 15:13:03 2026 -0700 Tier TP53 FLOSSIES BS2 by carrier count; ClinGen wording, per VCEP meeting. refs #37399 Replace the flat BS2 (-4) with the VCEP's carrier-count tiers: 2-3 carriers BS2_Supporting (-1), 4-7 BS2_Moderate (-2), 8+ BS2 (-4); a single carrier does not meet BS2. Carriers = allele_count - hom_count. tp53Flossies.py assigns the tier and color per observation; tp53ProvisionalClass.py reads the tier from the FLOSSIES bed (load_flossies_lookup returns label+points, bs2_evidence picks the strongest matching tier) and feeds the matching point value into the Provisional point sum. Also change the Provisional disclaimer from 'NOT a VCEP classification' to 'NOT a ClinGen classification' in the mouseover, docstring, and AUTOSQL, and drop a now-unused POINTS['BS2'] entry. diff --git src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py index 50982ce9e63..82d3ed6c5b2 100644 --- src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py +++ src/hg/makeDb/scripts/tp53/tp53ProvisionalClass.py @@ -10,31 +10,31 @@ * PP3 / BP4 (CSpec GN009 v2.4.0 Table S2) * BA1 / BS1 / PM2 (gnomAD v4.1 exomes, per CSpec FAF/grpmax thresholds) * BS2 (FLOSSIES healthy-women-over-70 cohort) * 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 VCEP classification — the warning is in every mouseover +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. """ 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 @@ -56,59 +56,58 @@ '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 POINTS = { 'PM1_Moderate': 2, 'PM1_Supporting': 1, 'PS3_Strong': 4, 'PS3': 4, 'PS3_Moderate': 2, 'PS3_Supporting': 1, 'BS3_Strong': -4, 'BS3': -4, 'BS3_Supporting': -1, 'PP3_Moderate': 2, 'PP3_moderate': 2, 'PP3': 1, 'BP4_Moderate': -2, 'BP4_moderate': -2, 'BP4': -1, 'BA1': 0, # stand-alone B -> not summed; forces class 'BS1': -4, 'PM2_Supporting': 1, - 'BS2': -4, 'No evidence': 0, 'Indeterminate': 0, } SPLICE_PP3_THRESHOLD = 0.2 # CSpec GN009: SpliceAI >= 0.2 -> PP3 splicing def bucket(pts): if pts >= 10: return 'Pathogenic' if 6 <= pts <= 9: return 'Likely pathogenic' if -6 <= pts <= -2: return 'Likely benign' if pts <= -7: return 'Benign' return 'Uncertain significance' CLASS_COLOR = { 'Pathogenic': '210,0,0', 'Likely pathogenic': '245,152,152', 'Uncertain significance': '0,0,136', 'Likely benign': '213,247,213', 'Benign': '0,210,0', } AUTOSQL = """table TP53ProvisionalClass -"TP53 NON-FINAL Provisional Classification by Tavtigian point sum (NOT a VCEP classification)" +"TP53 NON-FINAL Provisional Classification by Tavtigian point sum (NOT a ClinGen classification)" ( string chrom; "Reference sequence chromosome or scaffold" uint chromStart; "Start position in chromosome" uint chromEnd; "End position in chromosome" 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)" @@ -243,51 +242,56 @@ 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 +BS2_TIER_POINTS = {'BS2': -4, 'BS2_Moderate': -2, 'BS2_Supporting': -1} + + def load_flossies_lookup(db): - """Read TP53Flossies_<db>.bed and return set of (chrom, pos1based, ref, alt) - keys for rows with BS2 applicability = 'BS2'. Match is variant-specific - (a FLOSSIES observation supports BS2 only for the exact nt change observed, - not for any change at the same codon).""" - out = set() + """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: flds = line.rstrip("\n").split("\t") if len(flds) < 19: continue - if flds[9] != 'BS2': + pts = BS2_TIER_POINTS.get(flds[9]) + if pts is None: continue chrom = flds[0] start = int(flds[1]) ref = flds[17] alt = flds[18] - out.add((chrom, start + 1, ref, alt)) + out[(chrom, start + 1, ref, alt)] = (flds[9], pts) return out def pm1_code_and_points(wt, codon, alt, hotspot_occ): if codon in PM1_HARDCODED_CODONS: return ('PM1_Moderate', 2) count = hotspot_occ.get((codon, alt), 0) if count >= 10: return ('PM1_Moderate', 2) if count >= 2: return ('PM1_Supporting', 1) return ('', 0) def ps3_bs3_label_and_points(code): @@ -339,63 +343,65 @@ 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 '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') return ('', 0, '') -def bs2_observed(wt, codon, paths, flossies_lookup, tx): - """Return True if any (c.X>Y) path resolves to a FLOSSIES variant that - matches by exact genomic position AND ref/alt. 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.""" +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 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 - if (chrom, g + 1, g_ref, g_alt) in flossies_lookup: - return True - return False + hit = flossies_lookup.get((chrom, g + 1, g_ref, g_alt)) + if hit and hit[1] < best[1]: + best = hit + return best 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 VCEP classification — preliminary point-sum only." + "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): splice_section = "" if 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>" @@ -468,81 +474,81 @@ 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) af_lbl = "{} ({})".format(af_code, af_qty) if af_code else '' ba1 = (af_code == 'BA1') - # BS2 from FLOSSIES - bs2_applies = bs2_observed(wt, codon, paths, flossies_lookup, tx) - bs2_pts = POINTS['BS2'] if bs2_applies else 0 - bs2_lbl = "BS2 (-4 pts)" if bs2_applies else "Not observed" + # 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] 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("BS2 (-4)") + 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) 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" if bs2_applies else "-", + 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) flossies_lookup = load_flossies_lookup(db) print(" S3 entries: {} S2 entries: {} " @@ -555,42 +561,42 @@ 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() - bs2_cnt = 0 + 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] == 'BS2': - bs2_cnt += 1 + if flds[16] in BS2_TIER_POINTS: + bs2_cnt[flds[16]] += 1 print(" Class distribution:") for k, n in cnt.most_common(): print(" {}: {}".format(k, n)) print(" AF code distribution: {}".format(dict(af_cnt))) - print(" BS2 applied: {} rows".format(bs2_cnt)) + print(" BS2 applied: {} rows {}".format(sum(bs2_cnt.values()), dict(bs2_cnt))) 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) if __name__ == "__main__": main()