615721361f4baf75c0715bb931c5fc1015101622 lrnassar Tue Aug 11 18:21:56 2026 -0700 Address code-review findings on the Cardiomyopathy VCEP scripts. refs #37446 - Gate PM1 to missense variants per the CSpec ("applicable to missense variants"); a positional-only test wrongly gave synonymous/truncating/splice variants PM1 and let it collide with BA1/BP7. PM1 firing 1,293 -> 700. - Transcript-gate the Walsh-2019 ClinVar coordinate lookup so a classic-vs-MANE c.notation collision no longer mis-places TNNT2 R92Q (was drawn ~331 nt off with a different variant's VariationID); the gate applies only to the WALSH_TX genes. - Show the amino-acid change in the REVEL mouseover (computed from the MANE CDS) so the per-alt genomic-forward-strand score is not misread on minus-strand genes. - Resolve every build input relative to --output-dir (sibling track outputs and cmp_downloads sources) for otto portability; canonical build byte-identical. - Also key the diseaseTag off the counted PM1 code, and makedoc corrections (worked example REVEL/gnomAD values, PM1 count, universe and EvRepo notes). diff --git src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py index 5805fe469e1..d1449537f69 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPProvisionalClass.py @@ -330,30 +330,41 @@ 'splice_donor_variant', 'missense_variant', 'inframe_deletion', 'inframe_insertion', 'initiator_codon_variant', 'splice_region_variant', 'synonymous_variant', '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 (debug only)') args = ap.parse_args() + # Resolve every input relative to --output-dir (the sibling track outputs and the + # cmp_downloads sources are all produced under it), so a different --output-dir does + # not silently mix fresh outputs with stale inputs from the default WORKDIR. + global B3_BED, B4_BED, B1_BED, EVREPO_JSON, ANNOT_TSV + wd = args.output_dir + B3_BED = f'{wd}/cmpVCEPAFfrequencies/cmpVCEPAFfrequenciesHg38.bed' + B4_BED = f'{wd}/cmpVCEPRevel/cmpVCEPRevelHg38.bed' + B1_BED = f'{wd}/cmpVCEPClinDomains/cmpVCEPClinDomainsHg38.bed' + EVREPO_JSON = f'{wd}/cmp_downloads/erepo/cardiomyopathyVCEP_classifications.json' + ANNOT_TSV = f'{wd}/cmpVCEPAnnotate/cmpVCEPAnnotations.hg38.tsv' + out_dir = os.path.join(args.output_dir, 'cmpVCEPProvisionalClass') os.makedirs(out_dir, exist_ok=True) print(' [B.11 Variant Evidence 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'])) @@ -392,32 +403,35 @@ thr = '0.0002' if gene == 'MYBPC3' else '0.0001' af_ev = ('BS1_Strong', f'FAF95 {faf:.2e}, ≥ {thr}') elif v['af_code'] == 'PM2_supporting': codes.add('PM2_Supporting') af_ev = ('PM2_Supporting', f'FAF95 {faf:.2e}, rare (≤ 4e-05)') # REVEL PP3/BP4 - missense only if is_missense: rc = revel.get((chrom, v['start'], alt)) if rc: code, score = rc codes.add(code) thr = '≥ 0.70' if code.startswith('PP3') else '≤ 0.40' revel_ev = (code, f'{score} ({thr})') - # PM1 hotspot (HCM-calibrated) - pm1_hit = in_pm1_region(chrom, v['start'], pm1) + # PM1 hotspot (HCM-calibrated). The CSpec applies PM1 to MISSENSE variants only + # ("Applicable to missense variants ... in the specific regions listed"), so gate on + # consequence, not position alone - otherwise synonymous/truncating/splice variants in + # the hotspot codons wrongly earn PM1 (and collide with BA1/BS1/BP7). + pm1_hit = is_missense and in_pm1_region(chrom, v['start'], pm1) if pm1_hit: codes.add('PM1_Moderate') # 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') ps1_codon = codon if pm5: @@ -460,31 +474,32 @@ 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. splice_safety = 'yes' if sa_score >= SPLICE_SAFETY_THRESHOLD else 'no' n_features += 1 for c in codes: code_counts[c] += 1 color = TRACK_COLOR disease_tag = '' if gene in ('MYH7', 'TNNT2'): - disease_tag = 'PM1 HCM-calibrated' if pm1_hit else 'HCM/DCM' + # key off the counted code, not raw pm1_hit, so a suppressed PM1 is not labelled + disease_tag = 'PM1 HCM-calibrated' if 'PM1_Moderate' in codes else 'HCM/DCM' applied_str = ';'.join(sorted(codes)) or 'no codes' kind = variant_kind(so) # ---- evidence-first mouseover: the data leads; the ACMG code is a subordinate tag ---- ev = [] if af_ev: ev.append(f'<b>Population frequency (gnomAD v4.1):</b> {af_ev[1]} → supports {af_ev[0]}') if revel_ev: ev.append(f'<b>Missense predictor (REVEL):</b> {revel_ev[1]} → supports {revel_ev[0]}') 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 '