5b4e8221e9a3ce8c7cd18d4554ee95dc4b4abfbe lrnassar Thu Aug 13 17:28:29 2026 -0700 Restore BP4_Moderate strength in TP53 VCEP Curated track from Table S2. refs #37399 EvRepo's structured export emits a bare 'BP4' tag with no BP4_Moderate variant (unlike PP3_Moderate / PM2_Supporting, which it keeps), so the VCEP Curated track dropped the BP4 strength the panel actually applied. Restore it from the VCEP's Table S2 bioinformatic worksheet, matched on the exact nt change: tp53VCEPClinVar.py loads Table S2's BP4_Moderate set and upgrades a bare Met 'BP4' to 'BP4_Moderate' for the matching variant. Validated against the ClinVar interpretation prose: where that text is parseable its BP4 strength matches Table S2 for every missense variant. The enrichment changes 48 variants, all exactly BP4 -> BP4_Moderate with no off-target code changes, and removes all 34 BP4 strength disagreements in the EvRepo cross-check (full code agreement 37% -> 43%). diff --git src/hg/makeDb/scripts/tp53/tp53VCEPClinVar.py src/hg/makeDb/scripts/tp53/tp53VCEPClinVar.py index 75fc4276adf..6dd7c81a4ee 100644 --- src/hg/makeDb/scripts/tp53/tp53VCEPClinVar.py +++ src/hg/makeDb/scripts/tp53/tp53VCEPClinVar.py @@ -1,477 +1,518 @@ #!/usr/bin/env python3 """ TP53 VCEP Curated Variants track generator. Fetches TP53 VCEP classifications from the ClinGen Evidence Repository (EvRepo). EvRepo is the authoritative source for VCEP classifications (176 TP53 variants as of 2026-04-17) and is used here instead of ClinVar esearch because ClinVar's NCBI esearch index silently drops legacy VCV records like VCV 12374 / R175H. Emits bigBed 9+9 with the per-variant list of applied evidence codes plus the overall classification (P / LP / VUS / LB / B). hg38 and hg19 coordinates parsed directly from the HGVS list; liftOver fallback where absent. """ import argparse import html import json import os import re import subprocess import sys import tempfile import time import urllib.parse import urllib.request EVREPO_URL = ("https://erepo.genome.network/evrepo/api/classifications" "?gene=TP53&matchLimit=2000&format=json") CLINVAR_EFETCH = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi" +# The VCEP bioinformatic worksheet (Table S2). Used to restore the strength that +# EvRepo's structured export drops from BP4: EvRepo emits a bare "BP4" tag with +# no BP4_Moderate variant (unlike PP3_Moderate / PM2_Supporting, which it keeps). +# Validated against the ClinVar interpretation prose: where that text is +# parseable, its BP4 strength matches Table S2 for every missense variant. +SRC_S2 = "/hive/users/lrnassar/claude/RM37399/tp53_downloads/bioinformatic_worksheet.xlsx" +_S2_BP4MOD = None + + +def _norm_c(s): + """Strip a transcript prefix so 'NM_000546.6:c.886C>T' -> 'c.886C>T'.""" + if not s: + return '' + return re.sub(r'^NM_000546\.[0-9]*:', '', s.strip()) + + +def load_s2_bp4_moderate(): + """Return the set of c. notations that Table S2 assigns BP4_Moderate.""" + global _S2_BP4MOD + if _S2_BP4MOD is not None: + return _S2_BP4MOD + out = set() + try: + import openpyxl + wb = openpyxl.load_workbook(SRC_S2, data_only=True) + ws = wb["Supplementary Table S2"] + for row in ws.iter_rows(min_row=4, values_only=True): + hgvsc = row[0] + code = str(row[4]).strip() if row[4] is not None else '' + if isinstance(hgvsc, str) and code.lower() == 'bp4_moderate': + out.add(_norm_c(hgvsc)) + except Exception as ex: + log(" WARNING: Table S2 BP4-strength recovery unavailable: {}".format(ex)) + _S2_BP4MOD = out + return out + # HGVS genomic regex for hg38 (NC_000017.11) and hg19 (NC_000017.10) HGVS_G_RE = re.compile(r'^NC_0000(17)\.(\d+):g\.(\d+)([ACGT])>([ACGT])$') HGVS_G_INDEL_RE = re.compile(r'^NC_0000(17)\.(\d+):g\.(\d+)_?(\d+)?(del|ins|dup).*$') LIFTOVER_CHAINS = { 'hg19_to_hg38': '/cluster/data/hg19/bed/liftOver/hg19ToHg38.over.chain.gz', 'hg38_to_hg19': '/cluster/data/hg38/bed/liftOver/hg38ToHg19.over.chain.gz', } CHROM_SIZES = { 'hg19': '/cluster/data/hg19/chrom.sizes', 'hg38': '/cluster/data/hg38/chrom.sizes', } COLORS = { 'pathogenic': '210,0,0', 'likely pathogenic': '245,152,152', 'uncertain significance': '0,0,136', 'likely benign': '213,247,213', 'benign': '0,210,0', } DEFAULT_COLOR = '136,136,136' CLASSIFICATION_CANONICAL = { 'pathogenic': 'Pathogenic', 'likely pathogenic': 'Likely pathogenic', 'uncertain significance': 'Uncertain significance', 'likely benign': 'Likely benign', 'benign': 'Benign', } AUTOSQL = """table TP53VCEPCuratedVars "TP53 VCEP expert-panel-reviewed classifications from ClinGen EvRepo" ( string chrom; "Reference sequence chromosome or scaffold" uint chromStart; "Start position in chromosome" uint chromEnd; "End position in chromosome" string name; "HGVS (NM_000546.6 preferred)" uint score; "Not used, set to 0" char[1] strand; "Not used, set to ." uint thickStart; "Same as chromStart" uint thickEnd; "Same as chromEnd" uint reserved; "RGB color" string classification; "Final VCEP classification (P/LP/VUS/LB/B)" string hgvsc; "c. notation (NM_000546.6)" string hgvsp; "p. notation if available" string varId; "ClinVar variation ID" string caid; "ClinGen Canonical Allele ID" string publishedDate; "Date published to EvRepo" string metCodes; "Semicolon-separated list of Met evidence codes" string notMetCodes; "Count / summary of Not Met codes" lstring _mouseOver; "HTML mouseover" ) """ def log(msg): print(msg, file=sys.stderr) def escape_html(s): if not s: return "" # HTML-escape, then encode non-ASCII as numeric entities so UCSC # hgTracks' mouseover pipeline renders cleanly (raw UTF-8 ends up as # mojibake like 'â€"'). escaped = html.escape(str(s)) return "".join( ch if ord(ch) < 128 else "&#{};".format(ord(ch)) for ch in escaped ) def fetch_evrepo(): """Download the full TP53 VCEP classification list from EvRepo.""" log("Fetching TP53 VCEP classifications from EvRepo...") req = urllib.request.Request(EVREPO_URL, headers={'User-Agent': 'UCSC-kent/TP53-hub'}) with urllib.request.urlopen(req, timeout=60) as resp: data = json.loads(resp.read()) entries = data.get('variantInterpretations', []) log(" Found {} TP53 VCEP entries".format(len(entries))) return entries def pick_hgvs(hgvs_list): """Select preferred HGVS notations from the HGVS list. Returns dict with hg38_g, hg19_g (genomic NC format), hgvsc (NM_000546.6), hgvsp (NM_000546.6 with p.), display name. """ out = {'hg38_g': None, 'hg19_g': None, 'hgvsc': None, 'hgvsp_full': None, 'display': None} for h in hgvs_list: if h.startswith('NC_000017.11:g.'): out['hg38_g'] = h elif h.startswith('NC_000017.10:g.'): out['hg19_g'] = h elif h.startswith('NM_000546.6:c.'): out['hgvsc'] = h elif h.startswith('NM_000546.5(TP53):c.') or h.startswith('NM_000546.6(TP53):c.'): # The prose form includes p. — e.g. # "NM_000546.5(TP53):c.245C>T (p.Pro82Leu)" out['hgvsp_full'] = h if out['display'] is None: out['display'] = h if out['display'] is None: out['display'] = out['hgvsc'] or (hgvs_list[0] if hgvs_list else 'unknown') return out def hgvs_to_bed_coords(hgvs_g): """Convert NC_000017:g.POS[ACGT]>[ACGT] or indel-style g. to (start, end).""" if not hgvs_g: return None m = HGVS_G_RE.match(hgvs_g) if m: pos = int(m.group(3)) return (pos - 1, pos) m = HGVS_G_INDEL_RE.match(hgvs_g) if m: start = int(m.group(3)) end_s = m.group(4) end = int(end_s) if end_s else start kind = m.group(5) # For a dup, the coord is the duplicated region. For del, the deleted range. # For ins, zero-width between positions; use start-1 to start in BED. if kind == 'ins': return (start, start) return (start - 1, end) return None def parse_entry(entry): """Turn one EvRepo interpretation into our internal dict.""" hgvs_list = entry.get('hgvs', []) or [] picks = pick_hgvs(hgvs_list) hg38 = hgvs_to_bed_coords(picks['hg38_g']) hg19 = hgvs_to_bed_coords(picks['hg19_g']) outcome = None met = [] not_met_count = 0 for g in entry.get('guidelines', []) or []: if g.get('outcome'): outcome = g['outcome'].get('label') for a in g.get('agents', []) or []: for c in a.get('evidenceCodes', []) or []: label = c.get('label') or '' status = (c.get('status') or '').strip().lower() if status == 'met': met.append(label) elif status == 'not met': not_met_count += 1 cls_raw = outcome or '' cls = CLASSIFICATION_CANONICAL.get(cls_raw.strip().lower(), cls_raw) # Pull a short p. from the prose form, if present. hgvsp = '' if picks['hgvsp_full']: m = re.search(r'\(p\.[^)]+\)', picks['hgvsp_full']) if m: hgvsp = m.group(0)[1:-1] hgvsc = '' if picks['hgvsc']: hgvsc = picks['hgvsc'] elif picks['hgvsp_full']: # Fallback to c. from the prose form m = re.search(r'c\.[^\s(]+', picks['hgvsp_full']) if m: hgvsc = 'NM_000546.6:' + m.group(0) + # Restore the BP4 strength EvRepo drops (bare "BP4" -> "BP4_Moderate") from + # Table S2, matched on the exact nt change. + if 'BP4' in met and _norm_c(hgvsc) in load_s2_bp4_moderate(): + met = ['BP4_Moderate' if c == 'BP4' else c for c in met] + return { 'var_id': entry.get('variationId', ''), 'caid': entry.get('caid', ''), 'published': entry.get('publishedDate', ''), 'classification': cls, 'hgvsc': hgvsc, 'hgvsp': hgvsp, 'display': picks['display'], 'hg38_bed': hg38, # (start,end) or None 'hg19_bed': hg19, 'met_codes': met, 'not_met_count': not_met_count, } def liftover_coords(coords, chain): if not coords: return {} with tempfile.NamedTemporaryFile(mode='w', suffix='.bed', delete=False) as f: in_bed = f.name for vid, (chrom, s, e) in coords.items(): f.write("{}\t{}\t{}\t{}\n".format(chrom, s, e, vid)) out_bed = in_bed.replace('.bed', '.lifted.bed') un_bed = in_bed.replace('.bed', '.unmapped.bed') try: subprocess.run( ["liftOver", in_bed, chain, out_bed, un_bed], check=True, stderr=subprocess.DEVNULL) except Exception: pass lifted = {} if os.path.exists(out_bed): with open(out_bed) as f: for line in f: flds = line.strip().split('\t') if len(flds) >= 4: lifted[flds[3]] = (flds[0], int(flds[1]), int(flds[2])) for p in [in_bed, out_bed, un_bed]: if os.path.exists(p): os.remove(p) return lifted def create_bed(variants, assembly): lines = [] unmapped = [] # Which coord set is native on this assembly? native_key = 'hg38_bed' if assembly == 'hg38' else 'hg19_bed' other_key = 'hg19_bed' if assembly == 'hg38' else 'hg38_bed' chain = LIFTOVER_CHAINS['hg19_to_hg38'] if assembly == 'hg38' \ else LIFTOVER_CHAINS['hg38_to_hg19'] # Build needs_liftover map for variants missing native coords. needs = {} for v in variants: if v[native_key] is None and v[other_key] is not None: other_start, other_end = v[other_key] needs[v['var_id']] = ('chr17', other_start, other_end) lifted = liftover_coords(needs, chain) if needs else {} if needs: log(" liftOver ({} -> {}): {} lifted of {} needing".format( 'hg19' if assembly == 'hg38' else 'hg38', assembly, len(lifted), len(needs))) for v in variants: start, end = None, None if v[native_key] is not None: start, end = v[native_key] chrom = 'chr17' elif v['var_id'] in lifted: chrom, start, end = lifted[v['var_id']] else: unmapped.append(v) continue color = COLORS.get((v['classification'] or '').strip().lower(), DEFAULT_COLOR) cv_url = "https://www.ncbi.nlm.nih.gov/clinvar/variation/{}/".format(v['var_id']) met_summary = ", ".join(v['met_codes']) if v['met_codes'] else "(none Met)" mo = ( "<b>Final — VCEP ClinVar submission</b><br>" "<b>Variant:</b> {disp}<br>" "<b>Classification:</b> {cls}<br>" "<b>ClinGen CAid:</b> {caid}<br>" "<b>ClinVar ID:</b> <a href=\"{url}\" target=\"_blank\">{vid}</a><br>" "<b>Published:</b> {pub}<br>" "<b>Met evidence codes:</b> {met}<br>" "<b>Not Met codes:</b> {notmet}<br>" "<b>Source:</b> ClinGen Evidence Repository (EvRepo)" ).format( disp=escape_html(v['display']), cls=escape_html(v['classification']), caid=escape_html(v['caid']), url=cv_url, vid=v['var_id'], pub=escape_html(v['published']), met=escape_html(met_summary), notmet=v['not_met_count'], ) # Use HGVSp short if available, else c. notation, else full display. # ASCII-encode every free-text field to avoid UCSC mouseover mojibake. def _safe(s): if s is None: return '' return "".join( ch if ord(ch) < 128 else "&#{};".format(ord(ch)) for ch in str(s) ) name = v['hgvsp'] or v['hgvsc'] or v['display'] if len(name) > 200: name = name[:197] + "..." lines.append("\t".join([ chrom, str(start), str(end), _safe(name), '0', '.', str(start), str(end), color, v['classification'] or 'Unknown', _safe(v['hgvsc'] or ''), _safe(v['hgvsp'] or ''), v['var_id'], _safe(v['caid']), v['published'], _safe("; ".join(v['met_codes'])), str(v['not_met_count']), mo, ])) return lines, unmapped def build_assembly(variants, assembly, outdir): log("\n=== {} ===".format(assembly)) entries, unmapped = create_bed(variants, assembly) log(" mapped: {} unmapped: {}".format(len(entries), len(unmapped))) if not entries: return None bed = os.path.join(outdir, "TP53VCEPCuratedVars_{}.bed".format(assembly)) with open(bed, 'w') as f: f.write("\n".join(entries) + "\n") subprocess.run(["sort", "-k1,1", "-k2,2n", bed, "-o", bed], check=True) as_file = os.path.join(outdir, "TP53VCEPCuratedVars.as") bb = os.path.join(outdir, "TP53VCEPCuratedVars{}.bb".format(assembly.capitalize())) subprocess.run( ["bedToBigBed", "-as=" + as_file, "-type=bed9+9", "-tab", bed, CHROM_SIZES[assembly], bb], check=True) log(" wrote {}".format(bb)) return bb def write_tsv(variants, path): cols = ['var_id', 'caid', 'display', 'hgvsc', 'hgvsp', 'classification', 'published', 'met_codes', 'not_met_count', 'hg38_start', 'hg38_end', 'hg19_start', 'hg19_end'] with open(path, 'w') as f: f.write("\t".join(cols) + "\n") for v in variants: h38 = v['hg38_bed'] or (None, None) h19 = v['hg19_bed'] or (None, None) row = [ v['var_id'], v['caid'], v['display'], v['hgvsc'] or '', v['hgvsp'] or '', v['classification'] or '', v['published'], ';'.join(v['met_codes']), str(v['not_met_count']), str(h38[0] if h38[0] is not None else ''), str(h38[1] if h38[1] is not None else ''), str(h19[0] if h19[0] is not None else ''), str(h19[1] if h19[1] is not None else ''), ] f.write("\t".join(row) + "\n") def backfill_clinvar_for_unlinked(variants): """For records that came in without a ClinVar var_id, look them up by HGVSc and backfill var_id + hgvsp from ClinVar esummary. Some EvRepo records (e.g. complex indels like c.322_332delinsATTCA) come through the variantInterpretations endpoint with hgvsp empty and no variationId populated; we recover both via ClinVar's esearch by VarName + esummary. """ eutils_search = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi" eutils_summary = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi" fixed = 0 # Pace requests to stay well under NCBI's 3 req/sec limit and avoid 429s # on first-hit cold cache. time.sleep(1.0) for v in variants: if v['var_id'] or not v['hgvsc']: continue m = re.match(r'NM_000546(?:\.\d+)?:(c\.[^\s(]+)', v['hgvsc']) if not m: continue c_form = m.group(1) success = False for attempt in range(3): try: url = (eutils_search + "?db=clinvar&retmode=json&term=" + urllib.parse.quote('"{}"[VarName] AND tp53[gene]'.format(c_form))) with urllib.request.urlopen(url, timeout=20) as resp: hits = json.loads(resp.read())['esearchresult']['idlist'] if not hits: break vid = hits[0] time.sleep(0.5) url = "{}?db=clinvar&id={}&retmode=json".format(eutils_summary, vid) with urllib.request.urlopen(url, timeout=20) as resp: summ = json.loads(resp.read())['result'][vid] v['var_id'] = vid title = summ.get('title', '') mp = re.search(r'\(p\.[^)]+\)', title) if mp and not v['hgvsp']: v['hgvsp'] = mp.group(0)[1:-1] fixed += 1 success = True break except urllib.error.HTTPError as e: if e.code == 429 and attempt < 2: time.sleep(2.0 * (attempt + 1)) continue log(" backfill skipped for {}: {}".format(c_form, e)) break except Exception as e: log(" backfill skipped for {}: {}".format(c_form, e)) break if success: time.sleep(0.5) if fixed: log("Backfilled ClinVar var_id for {} previously-unlinked records".format(fixed)) def fetch_all_variants(): entries = fetch_evrepo() variants = [parse_entry(e) for e in entries] variants = [v for v in variants if v['classification']] log("Parsed {} variants with classifications".format(len(variants))) backfill_clinvar_for_unlinked(variants) return variants def main(): p = argparse.ArgumentParser(description=__doc__) p.add_argument('-o', '--output-dir', default=os.path.dirname(os.path.abspath(__file__)), help='Output directory (default: script directory)') p.add_argument('--db', action='append', help='hg38 or hg19 (repeat); default both') args = p.parse_args() outdir = args.output_dir os.makedirs(outdir, exist_ok=True) dbs = args.db if args.db else ['hg38', 'hg19'] variants = fetch_all_variants() if not variants: log("ERROR: no variants parsed from EvRepo") sys.exit(1) tsv = os.path.join(outdir, "tp53_vcep_clinvar_variants.tsv") write_tsv(variants, tsv) log("TSV: {}".format(tsv)) as_file = os.path.join(outdir, "TP53VCEPCuratedVars.as") with open(as_file, 'w') as f: f.write(AUTOSQL) for db in dbs: build_assembly(variants, db, outdir) from collections import Counter cls_counts = Counter(v['classification'] for v in variants) log("\nClassification counts:") for c, n in cls_counts.most_common(): log(" {}: {}".format(c, n)) if __name__ == "__main__": main()