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'Population frequency (gnomAD v4.1): {af_ev[1]} → supports {af_ev[0]}')
if revel_ev:
ev.append(f'Missense predictor (REVEL): {revel_ev[1]} → supports {revel_ev[0]}')
if ps1_codon:
ev.append(f'Known variants at residue {ps1_codon} (EvRepo): the same amino-acid change is '
f'VCEP P/LP → supports PS1_Strong')
if pm5_codon:
ev.append(f'Known variants at residue {pm5_codon} (EvRepo): a different VCEP P/LP change '