0d12a0f167c83fb30ff8e001de6aff57c0982447 lrnassar Wed Jul 8 15:55:13 2026 -0700 Reframe the Cardiomyopathy VCEP provisional track to a Computable ACMG Criteria Summary. refs #37446 Per the CM VCEP chair, the cmpVCEPProvisionalClass track no longer computes an overall ACMG classification. It now lists, per variant, only the computable ACMG codes that are triggered, each with its underlying value: gnomAD FAF95 for BA1/BS1/PM2_supporting, REVEL score for PP3/BP4, the PM1 hotspot region, the EvRepo codon reference for PS1/PM5, the NMD-escape basis for PM4, and SpliceAI/phyloP for BP7. Display is a single neutral color; a SpliceAI score >= 0.20 is recorded as an informational splice flag rather than overriding a call. The GN002 combining logic is retained in the script as classify() but is retired and no longer called. Makedoc B.11 updated to match. diff --git src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py index 435f5fcaa64..d712286d124 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py @@ -1,22 +1,22 @@ #!/usr/bin/env python3 """ -B.11 — NON-FINAL Provisional Classification track (NOT a VCEP classification). +B.11: Computable ACMG Criteria Summary track (NOT a VCEP classification). For every gnomAD-observed variant in the 8 cardiomyopathy gene CDS regions ±20 nt -splice padding, computes a NON-FINAL provisional ACMG/AMP classification from the -subset of evidence codes a hub can compute automatically: +splice padding, lists the subset of ACMG/AMP evidence codes a hub can compute +automatically. No overall classification is calculated: - BA1 / BS1 / PM2_Supporting (gnomAD v4.1 FAF95; B.3) - PP3 / BP4 (REVEL; B.4; missense only, per hgVai consequence) - PM1 (B.1 hotspot regions; HCM-scoped; NOT combined with PM5) - PS1 / PM5 (EvRepo P/LP reference, LEAVE-ONE-OUT; see caveat below) - PM4 (NMD-escaping truncating variants, non-MYBPC3; CSpec disease-specific) - BP7 (synonymous + SpliceAI no-impact + not conserved) - Splice safety net (SpliceAI ≥ 0.20 overrides a benign-leaning call to VUS) - HCM/DCM tag (MYH7, TNNT2: PM1 is HCM-calibrated) CONSEQUENCE/CODON/AA come from the Phase-1 hgVai annotation TSV (cmpVCEPAnnotate). HONESTY / KNOWN LIMITS (surfaced to the VCEP, not hidden): * This mockup CANNOT compute the clinical/functional codes (PS2, PS3, PS4, PP1, PP4, BS3, BS4). Many true P/LP calls rest on those, so this track structurally under-calls pathogenicity. It is best read as a benign/VUS-axis + "flag for @@ -65,72 +65,67 @@ 'NC_000012.12': 'chr12', 'NC_000014.9': 'chr14', 'NC_000015.10': 'chr15', 'NC_000019.10': 'chr19', } # Variants with established splicing impact — excluded from the PS1/PM5 reference per GN002 PS1. PS1_SPLICE_EXCLUDE = {'NM_000256.3:c.2308G>A'} NC_G_RE = re.compile(r'^(NC_\d+\.\d+):g\.(\d+)([ACGT]+)>([ACGT]+)$') PROT_MISSENSE_RE = re.compile(r'p\.([A-Z][a-z]{2})(\d+)([A-Z][a-z]{2})') AA3TO1 = { 'Ala': 'A', 'Arg': 'R', 'Asn': 'N', 'Asp': 'D', 'Cys': 'C', 'Gln': 'Q', 'Glu': 'E', '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', 'Ter': '*', } -COLORS = { - 'Pathogenic': '210,0,0', - 'Likely Pathogenic': '245,152,152', - 'Uncertain Significance': '0,0,136', - 'Likely Benign': '213,247,213', - 'Benign': '0,210,0', -} +TRACK_COLOR = '91,107,122' # neutral slate; no classification encoded (evidence-only track) CHROM_SIZES = {'hg38': '/cluster/data/hg38/chrom.sizes', 'hg19': '/cluster/data/hg19/chrom.sizes'} LIFTOVER_HG38_TO_HG19 = '/cluster/data/hg38/bed/liftOver/hg38ToHg19.over.chain.gz' AUTOSQL = """table cmpVCEPProvisionalClass -"NON-FINAL provisional ACMG classification (computable-codes subset; NOT a VCEP classification)" +"Computable ACMG criteria per variant (evidence summary; NOT a VCEP classification)" ( string chrom; "Chromosome" uint chromStart; "Position" uint chromEnd; "End" string name; "Display name" uint score; "0" char[1] strand; "Strand" uint thickStart; "Same" uint thickEnd; "Same" - uint itemRgb; "ACMG color" + uint itemRgb; "Display color (neutral; no classification encoded)" string gene; "Gene" string refAllele; "Ref" string altAllele; "Alt" string variantKind; "Predicted consequence (hgVai)" - string classification; "NON-FINAL provisional call" - string appliedCodes; "ACMG codes met (semicolon-separated, with strengths)" - string ruleMatch; "Combination rule that fired" + string appliedCodes; "Computable ACMG codes triggered (semicolon-separated, with strengths)" string diseaseTag; "HCM/DCM phenotype scoping note (MYH7, TNNT2)" string codeNotes; "Suppressed/contested codes (e.g. PM5 not combined with PM1)" - string splice_safety; "yes if SpliceAI >= 0.20 triggered the safety net" + string splice_safety; "yes if SpliceAI >= 0.20 (possible splice impact; informational)" lstring _mouseOver; "Tooltip" ) """ # ============================================================ # Combination rules — transcribed verbatim from CSpec GN002 # ============================================================ +# RETIRED 2026-07-08 (per CM VCEP chair L. Bronicki): this track no longer computes an overall +# ACMG classification, only the computable codes that fire. The GN002 combining logic below is +# kept for reference and is intentionally NOT called. def classify(codes): """Apply the Cardiomyopathy CSpec (GN002) combining rules. `codes` is an iterable of code strings carrying explicit strengths where relevant, e.g. {'PM1_Moderate', 'PP3_Supporting', 'PS1_Strong', 'BA1', 'BS1_Strong'}. Returns (classification, rule_match). Strong: PS1, PS2, PS3, PS4, PP1_Strong Moderate: PS3_Moderate, PS4_Moderate, PM1, PM4(_Moderate), PM5, PM6, PP1_Moderate Supporting: PS3_Supporting, PS4_Supporting, PM2_Supporting, PM5_Supporting, PP1, PP3, PM4_Supporting Benign: BA1 (stand-alone); BS1/BS3/BS4 (strong); BP4/BP7 (supporting) """ p_s = p_m = p_sup = 0 b_sa = b_s = b_sup = 0 for code in codes: @@ -204,31 +199,31 @@ if len(f) < 18: continue rows.append({'chrom': f[0], 'start': int(f[1]), 'end': int(f[2]), 'strand': f[5], 'gene': f[9], 'ref': f[10], 'alt': f[11], 'faf95': float(f[12]), 'af_code': f[15]}) print(f' B.3 universe: {len(rows)} variants', file=sys.stderr) return rows def load_b4_revel_lookup(): lookup = {} for line in open(B4_BED): f = line.rstrip('\n').split('\t') if len(f) < 14: continue - lookup[(f[0], int(f[1]), f[10])] = f[12] # (chrom, start, alt) -> PP3_/BP4_Supporting + lookup[(f[0], int(f[1]), f[10])] = (f[12], f[11]) # (chrom,start,alt) -> (code, REVEL score) print(f' REVEL lookup: {len(lookup)} keys', file=sys.stderr) return lookup def load_pm1_intervals(): intervals = [] for line in open(B1_BED): f = line.rstrip('\n').split('\t') intervals.append((f[0], int(f[1]), int(f[2]))) print(f' PM1 intervals: {len(intervals)}', file=sys.stderr) return intervals def in_pm1_region(chrom, pos0, intervals): for c, s, e in intervals: @@ -354,201 +349,213 @@ 'intron_variant', '5_prime_UTR_variant', '3_prime_UTR_variant'): if k in so: return k return ','.join(sorted(so)) if so else 'unknown' def main(): ap = argparse.ArgumentParser() ap.add_argument('--db', action='append', required=True, choices=['hg38', 'hg19']) ap.add_argument('--output-dir', required=True) ap.add_argument('--no-spliceai', action='store_true', help='Skip SpliceAI/BP7-conservation (debug only)') args = ap.parse_args() out_dir = os.path.join(args.output_dir, 'cmpVCEPProvisionalClass') os.makedirs(out_dir, exist_ok=True) - print(' [B.11 NON-FINAL Provisional Classification]') + print(' [B.11 Computable ACMG Criteria Summary]') b3 = load_b3_variants() revel = load_b4_revel_lookup() pm1 = load_pm1_intervals() ann = load_annotation() evref = load_evrepo_reference() regions = [] nmd_junction = {} # gene -> cDNA coord of the last exon-exon junction (transcript len - last exon len) for gene in OUR_GENES: m = parse_mane_record(gene) regions.append((m['chrom'], m['chromStart'], m['chromEnd'])) bs = m['blockSizes'] last_exon = bs[0] if m['strand'] == '-' else bs[-1] # 3'-most transcript exon nmd_junction[gene] = sum(bs) - last_exon if args.no_spliceai: spliceai, phylop = {}, {} else: spliceai = batch_spliceai(regions) phylop = batch_phylop(regions) - counts = defaultdict(int) + n_features = 0 code_counts = defaultdict(int) bed_lines = [] for v in b3: gene = v['gene'] if gene not in OUR_GENES: continue chrom, pos1, ref, alt = v['chrom'], v['start'] + 1, v['ref'], v['alt'] a = ann.get((chrom, pos1, ref, alt), {}) so = a.get('so', set()) is_missense = 'missense_variant' in so is_synonymous = 'synonymous_variant' in so codes = set() + code_why = {} notes = [] + faf = v['faf95'] # gnomAD AF (B.3) if v['af_code'] == 'BA1': codes.add('BA1') + code_why['BA1'] = f'gnomAD FAF95 (popmax) {faf:.2e} ≥ 0.001' elif v['af_code'] == 'BS1': codes.add('BS1_Strong') + thr = '0.0002' if gene == 'MYBPC3' else '0.0001' + code_why['BS1_Strong'] = f'gnomAD FAF95 (popmax) {faf:.2e} ≥ {thr}' elif v['af_code'] == 'PM2_supporting': codes.add('PM2_Supporting') + code_why['PM2_Supporting'] = f'gnomAD FAF95 (popmax) {faf:.2e} ≤ 4e-05 (rare)' # REVEL PP3/BP4 — missense only if is_missense: rc = revel.get((chrom, v['start'], alt)) if rc: - codes.add(rc) + code, score = rc + codes.add(code) + thr = '≥ 0.70' if code.startswith('PP3') else '≤ 0.40' + code_why[code] = f'REVEL {score} ({thr})' # PM1 hotspot (HCM-calibrated) pm1_hit = in_pm1_region(chrom, v['start'], pm1) if pm1_hit: codes.add('PM1_Moderate') + code_why['PM1_Moderate'] = f'in the {gene} PM1 hotspot region (HCM-calibrated)' # PS1 / PM5 — EvRepo P/LP reference, LEAVE-ONE-OUT (exclude self by genomic key) if is_missense and a.get('codon') and a.get('aaAlt'): codon, aaalt = a['codon'], a['aaAlt'] gkey = (chrom, pos1, ref, alt) ps1 = any(e['gene'] == gene and e['codon'] == codon and e['alt_aa1'] == aaalt and e['gkey'] != gkey for e in evref) pm5 = any(e['gene'] == gene and e['codon'] == codon and e['alt_aa1'] != aaalt and e['gkey'] != gkey for e in evref) if ps1: codes.add('PS1_Strong') + code_why['PS1_Strong'] = f'same amino-acid change as a VCEP EvRepo P/LP variant at codon {codon}' if pm5: codes.add('PM5_Moderate') + code_why['PM5_Moderate'] = f'a different missense at codon {codon} is classified P/LP in the VCEP EvRepo set' # PM4 — NMD-escaping truncating, non-MYBPC3 (CSpec disease-specific; PVS1 N/A for these genes). # NMD escapes if the PTC is in the last exon OR within 50 nt of the last exon-exon junction # (J = transcript length - last exon length); cDNA position is in transcript orientation. 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 last exon-exon junction' - notes.append(f'<b>PM4:</b> NMD-escaping truncating ({where}); strength per VCEP') + where = 'last exon' if last else 'within 50 nt of the last exon-exon junction' + code_why['PM4_Supporting'] = f'NMD-escaping truncating variant ({where})' elif gene != 'MYBPC3' and 'stop_lost' in so: codes.add('PM4_Supporting') - notes.append('<b>PM4:</b> stop-loss; strength per VCEP') + code_why['PM4_Supporting'] = 'stop-loss variant' # BP7 — synonymous + SpliceAI no-impact + not conserved (conservation cutoff PROVISIONAL) sa_score = spliceai.get((chrom, pos1, ref, alt), 0.0) if is_synonymous and not args.no_spliceai: phy = phylop.get((chrom, pos1)) if sa_score < BP7_SPLICE_MAX and phy is not None and phy <= BP7_PHYLOP_MAX: codes.add('BP7_Supporting') - notes.append(f'<b>BP7:</b> SpliceAI {sa_score:.2f}<{BP7_SPLICE_MAX}, phyloP {phy:.2f}<=0 (cutoff provisional — VCEP to confirm)') + code_why['BP7_Supporting'] = (f'synonymous; SpliceAI {sa_score:.2f} < {BP7_SPLICE_MAX}; ' + f'phyloP {phy:.2f} ≤ 0 (conservation cutoff provisional)') # 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. if 'PM1_Moderate' in codes and 'PM5_Moderate' in codes: codes.discard('PM1_Moderate') notes.append('<b>PM1 suppressed:</b> CSpec says PM5 (variant-specific) is preferred over PM1') # CSpec is silent on PM1+PS1; flag as possible double-count for VCEP (do not suppress) if 'PM1_Moderate' in codes and 'PS1_Strong' in codes: notes.append('<b>Note:</b> PM1+PS1 co-occur — possible double-counting (VCEP question)') - classification, rule_match = classify(codes) - - # Splice safety net: SpliceAI >= 0.20 on a benign-leaning call -> VUS + # Splice signal (informational): SpliceAI >= 0.20 flags possible splice impact. + # No longer overrides a call — this track computes no overall classification. splice_safety = 'no' - if sa_score >= SPLICE_SAFETY_THRESHOLD and classification in ('Likely Benign', 'Benign'): - notes.append(f'<b>Splice safety net:</b> SpliceAI {sa_score:.2f} >= {SPLICE_SAFETY_THRESHOLD} overrode {classification} → VUS') - classification = 'Uncertain Significance' - rule_match = f'splice safety net (SpliceAI={sa_score:.2f})' + if sa_score >= SPLICE_SAFETY_THRESHOLD: + notes.append(f'<b>Splice flag:</b> SpliceAI {sa_score:.2f} >= {SPLICE_SAFETY_THRESHOLD} (possible splice impact; informational)') splice_safety = 'yes' - counts[classification] += 1 + n_features += 1 for c in codes: code_counts[c] += 1 - color = COLORS.get(classification, '136,136,136') + color = TRACK_COLOR disease_tag = '' if gene in ('MYH7', 'TNNT2'): disease_tag = 'PM1 HCM-calibrated' if pm1_hit else 'HCM/DCM' applied_str = ';'.join(sorted(codes)) or 'no codes' notes_str = ' | '.join(notes) kind = variant_kind(so) - mo = [f'<b>NON-FINAL Provisional Classification</b> (NOT a VCEP classification)<br>', + mo = [f'<b>Computable ACMG Criteria Summary</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'<br><b>Consequence:</b> {kind}<br>') - mo.append(f'<b>Codes met:</b> {applied_str}<br>') - mo.append(f'<b>Rule:</b> {rule_match}<br>') - mo.append(f'<b>Provisional call:</b> {classification}<br>') + if codes: + mo.append('<b>Computable codes triggered:</b><br>') + for c in sorted(codes): + why = code_why.get(c, '') + mo.append(f' <b>{c}:</b> {why}<br>' if why else f' <b>{c}</b><br>') + else: + mo.append('<b>Computable codes triggered:</b> none<br>') if notes_str: mo.append(f'<span style="color:#a00">{notes_str}</span><br>') - if 'PS1_Strong' in codes or 'PM5_Moderate' in codes: - mo.append('<i>PS1/PM5 from VCEP EvRepo P/LP set (leave-one-out); reference DB is an open VCEP question.</i><br>') - mo.append('<i>Mockup of computable codes only — excludes clinical/functional PS2/PS3/PS4/PP1/PP4/BS3/BS4; ' - 'best used to flag variants for expert review, not as a classification.</i>') + mo.append('<br><i>Shows only the computable ACMG codes that are triggered. Excludes clinical/functional ' + 'PS2/PS3/PS4/PP1/PP4/BS3/BS4.</i>') mouseover = ''.join(mo) - name = f'{gene}_{pos1}_{ref}>{alt}_{classification[:1]}' + 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, - classification, applied_str, rule_match, disease_tag, notes_str, splice_safety, mouseover, + applied_str, disease_tag, notes_str, splice_safety, mouseover, ])) - print(f' classifications: {dict(counts)}') + 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}') if 'hg38' in args.db: hg38_bb = os.path.join(out_dir, 'cmpVCEPProvisionalClassHg38.bb') - subprocess.run(['bedToBigBed', '-tab', '-type=bed9+11', '-as=' + as_path, + subprocess.run(['bedToBigBed', '-tab', '-type=bed9+9', '-as=' + as_path, hg38_bed, CHROM_SIZES['hg38'], hg38_bb], check=True) print(f' hg38 bigBed: {hg38_bb}') if 'hg19' in args.db: hg19_bed = os.path.join(out_dir, 'cmpVCEPProvisionalClassHg19.bed') unmapped = hg19_bed + '.unmapped' subprocess.run(['liftOver', '-bedPlus=9', '-tab', hg38_bed, LIFTOVER_HG38_TO_HG19, hg19_bed, unmapped], check=True) if os.path.getsize(unmapped) > 0: n = sum(1 for line in open(unmapped) if not line.startswith('#')) print(f' WARNING: {n} unmapped in hg19 liftOver: {unmapped}', file=sys.stderr) hg19_bb = os.path.join(out_dir, 'cmpVCEPProvisionalClassHg19.bb') - subprocess.run(['bedToBigBed', '-tab', '-type=bed9+11', '-as=' + as_path, + subprocess.run(['bedToBigBed', '-tab', '-type=bed9+9', '-as=' + as_path, hg19_bed, CHROM_SIZES['hg19'], hg19_bb], check=True) print(f' hg19 bigBed: {hg19_bb}') if __name__ == '__main__': main()