6796b337521a06767307b5c273bf89903d5a3bd6 lrnassar Tue Aug 11 17:05:08 2026 -0700 Fix Variant Evidence Summary SpliceAI/BP7 issues per code review. refs #37446 - Use the released spliceAIsnvsMasked.bb instead of a superseded symlink; the older /gbdb/hg38/bbi/spliceAi.bb pointed at a stale file the live track no longer serves. - Distinguish a SpliceAI lookup miss from a measured 0.00 in the mouseover. A miss now reads "no record (below the 0.02 reporting floor)"; BP7 still applies because the file's 0.02 floor means a missing score is below 0.1. - Correct the BP7 mouseover threshold text from "<= 0.1" to "< 0.1" to match the code, docstring, and makedoc. - Record PM1 suppression in the codeNotes field so it matches the autoSql description instead of leaving the field empty on every row. diff --git src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py index 0007e034422..4d991c786b6 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py @@ -33,31 +33,31 @@ import argparse, json, os, re, subprocess, sys from collections import defaultdict sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from cmpVCEPClinDomains import parse_mane_record OUR_GENES = ['MYH7', 'MYBPC3', 'TNNT2', 'TNNI3', 'TPM1', 'ACTC1', 'MYL2', 'MYL3'] WORKDIR = '/hive/users/lrnassar/claude/RM37446' B3_BED = f'{WORKDIR}/cmpVCEPAFfrequencies/cmpVCEPAFfrequenciesHg38.bed' B4_BED = f'{WORKDIR}/cmpVCEPRevel/cmpVCEPRevelHg38.bed' B1_BED = f'{WORKDIR}/cmpVCEPClinDomains/cmpVCEPClinDomainsHg38.bed' EVREPO_JSON = f'{WORKDIR}/cmp_downloads/erepo/cardiomyopathyVCEP_classifications.json' ANNOT_TSV = f'{WORKDIR}/cmpVCEPAnnotate/cmpVCEPAnnotations.hg38.tsv' -SPLICEAI_BB = '/gbdb/hg38/bbi/spliceAi.bb' +SPLICEAI_BB = '/gbdb/hg38/bbi/spliceAIsnvsMasked.bb' # released masked-SNV file (what the live UCSC track serves) # Per-gene thresholds (from CSpec - NOT invented here) BS1_THRESHOLDS = {'MYBPC3': 0.0002} DEFAULT_BS1 = 0.0001 BA1_THRESHOLD = 0.001 PM2_SUPPORTING_THRESHOLD = 0.00004 SPLICE_SAFETY_THRESHOLD = 0.20 # SpliceAI delta flagged as possible splice impact (informational) # BP7: synonymous with no predicted splice impact. The CM VCEP set this at SpliceAI < 0.1 and # removed the conservation requirement, per Walker 2023 (PMID 37352859). BP7_SPLICE_MAX = 0.10 # NC_ accession (hg38) -> chrom, for parsing EvRepo genomic HGVS (leave-one-out keys) NC_HG38 = { 'NC_000001.11': 'chr1', 'NC_000003.12': 'chr3', 'NC_000011.10': 'chr11', 'NC_000012.12': 'chr12', 'NC_000014.9': 'chr14', 'NC_000015.10': 'chr15', @@ -426,31 +426,36 @@ if gene != 'MYBPC3' and (so & TRUNCATING_SO): cdna, J = a.get('cdnaPos'), nmd_junction.get(gene) nmd_escape = (a.get('exonNum') is not None and a.get('exonNum') == a.get('exonTotal')) \ or (cdna is not None and J is not None and cdna > J - 50) if nmd_escape: codes.add('PM4_Supporting') last = a.get('exonNum') == a.get('exonTotal') where = 'last exon' if last else 'within 50 nt of the last exon-exon junction' pm4_ev = f'truncating, escapes NMD ({where})' elif gene != 'MYBPC3' and 'stop_lost' in so: codes.add('PM4_Supporting') pm4_ev = 'stop-loss variant' # BP7 - synonymous with no predicted splice impact (SpliceAI < 0.1, per Walker 2023 # PMID 37352859). The CM VCEP removed the conservation requirement, so no phyloP gate. - sa_score = spliceai.get((chrom, pos1, ref, alt), 0.0) + # The SpliceAI file has a 0.02 reporting floor, so a lookup miss means the true + # score is below 0.02 (hence below 0.1) - BP7 still applies. Track presence so the + # mouseover does not print a missing record as a measured "0.00". + sa_hit = spliceai.get((chrom, pos1, ref, alt)) + sa_present = sa_hit is not None + sa_score = sa_hit if sa_present else 0.0 if is_synonymous and not args.no_spliceai: if sa_score < BP7_SPLICE_MAX: codes.add('BP7_Supporting') # CSpec exclusion: PM1 must NOT be combined with PM5. GN002 PM5: "use of PM5 is most # appropriate since it is variant specific" -> keep PM5, drop PM1. pm1_suppressed = False if 'PM1_Moderate' in codes and 'PM5_Moderate' in codes: codes.discard('PM1_Moderate') pm1_suppressed = True # CSpec is silent on PM1+PS1; flag as a possible double-count (do not suppress). if 'PM1_Moderate' in codes and 'PS1_Strong' in codes: notes.append('PM1 and PS1 co-occur here, a possible double-count (open VCEP question)') # Splice signal (informational): SpliceAI >= 0.20 flags possible splice impact. @@ -477,57 +482,65 @@ if ps1_codon: ev.append(f'<b>Known variants at residue {ps1_codon} (EvRepo):</b> the same amino-acid change is ' f'VCEP P/LP → supports PS1_Strong') if pm5_codon: ev.append(f'<b>Known variants at residue {pm5_codon} (EvRepo):</b> a different VCEP P/LP change ' f'exists at this residue → supports PM5_Moderate') if pm1_hit: if pm1_suppressed: ev.append(f'<b>PM1 hotspot region:</b> within the {gene} hotspot (HCM only); not counted here, ' f'PM5 (variant-specific) is preferred per CSpec') else: ev.append(f'<b>PM1 hotspot region:</b> within the {gene} hotspot (HCM only) → supports PM1_Moderate') if pm4_ev: ev.append(f'<b>Protein-truncating:</b> {pm4_ev} → supports PM4_Supporting') if 'BP7_Supporting' in codes: - ev.append(f'<b>Synonymous, splicing (SpliceAI):</b> {sa_score:.2f}, no predicted impact ' - f'(≤ 0.1, per Walker 2023) → supports BP7_Supporting') + sa_txt = f'{sa_score:.2f}' if sa_present else 'no record (below the 0.02 reporting floor)' + ev.append(f'<b>Synonymous, splicing (SpliceAI):</b> {sa_txt}, no predicted impact ' + f'(< 0.1, per Walker 2023) → supports BP7_Supporting') if splice_safety == 'yes': ev.append(f'<b>Splicing (SpliceAI):</b> {sa_score:.2f}, possible splice impact (informational)') mo = ['<b>Variant evidence</b> (not a VCEP classification)<br>', f'<b>{gene}</b> · {chrom}:{pos1} {ref}>{alt}'] if a.get('hgvsp'): mo.append(f' · <i>{a["hgvsp"]}</i>') mo.append(f' · {kind}<br>') if ev: for line in ev: mo.append(line + '<br>') else: mo.append('No computable evidence at this position.<br>') if notes: mo.append('<span style="color:#a00">' + ' | '.join(notes) + '</span><br>') mo.append('<br><i>Shows the computable evidence and the ACMG criteria it would support. ' 'Not a VCEP review or classification; clinical/functional evidence ' '(PS2/PS3/PS4/PP1/PP4/BS3/BS4) is not included.</i>') mouseover = ''.join(mo) + # codeNotes data field: contested/double-count notes plus any suppressed code. + # The mouseover renders suppression inline as an evidence line, so it is kept out + # of the red-span `notes` above to avoid duplication; it belongs in the data field. + code_notes = list(notes) + if pm1_suppressed: + code_notes.append('PM1 present but not combined with PM5 (PM5 preferred per CSpec)') + name = f'{gene}_{pos1}_{ref}>{alt}' bed_lines.append('\t'.join([ chrom, str(v['start']), str(v['end']), name, '0', v['strand'], str(v['start']), str(v['end']), color, gene, ref, alt, kind, - applied_str, disease_tag, ' | '.join(notes), splice_safety, mouseover, + applied_str, disease_tag, ' | '.join(code_notes), splice_safety, mouseover, ])) print(f' features: {n_features}') print(f' code firing counts: {dict(sorted(code_counts.items()))}') bed_lines.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) as_path = os.path.join(out_dir, 'cmpVCEPProvisionalClass.as') with open(as_path, 'w') as f: f.write(AUTOSQL) hg38_bed = os.path.join(out_dir, 'cmpVCEPProvisionalClassHg38.bed') with open(hg38_bed, 'w') as f: f.write('\n'.join(bed_lines) + '\n') print(f' wrote {len(bed_lines)} BED features -> {hg38_bed}')