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/cmpVCEPWalsh2019.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPWalsh2019.py index 81881ac7f2a..035be445800 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPWalsh2019.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPWalsh2019.py @@ -196,30 +196,38 @@ n_rows = 0 with gzip.open(VARIANT_SUMMARY, 'rt') as fh: for line in fh: if line.startswith('#'): continue f = line.rstrip('\n').split('\t') if len(f) < 31: continue name = f[2] gene = f[4] if gene not in OUR_GENES: continue m = CV_NAME_RE.match(name) if not m: continue + # For genes with a specified Walsh transcript, match only ClinVar rows on THAT + # transcript. Otherwise a classic-vs-MANE c.notation collision mis-places the + # variant (TNNT2 classic c.275G>A [R92Q] vs MANE c.275G>A [G92E]); for TNNT2 this + # drops the MANE-named rows so the entry falls to the hgvsToVcf-on-classic path. + # Genes not in WALSH_TX (ACTC1/MYL2/MYL3/TPM1) use MANE == ClinVar's Name transcript, + # so there is no collision and the (gene, c.notation) match stands. + if gene in WALSH_TX and m.group(1) != WALSH_TX[gene]: + continue cdna = 'c.' + m.group(3) assembly = f[16] db = 'hg38' if assembly == 'GRCh38' else ('hg19' if assembly == 'GRCh37' else None) if db is None: continue try: start1 = int(f[19]); stop1 = int(f[20]) except ValueError: continue chrom_num = f[18] key = (gene, cdna) lookup.setdefault(key, {})[db] = { 'chrom': f'chr{chrom_num}', 'start': start1 - 1, 'end': stop1, @@ -290,30 +298,34 @@ def make_bigbed(bed_path, db, as_path, bb_path): cmd = ['bedToBigBed', '-tab', '-type=bed9+12', '-as=' + as_path, bed_path, CHROM_SIZES[db], bb_path] print(f' $ {" ".join(cmd)}') subprocess.run(cmd, check=True) def main(): ap = argparse.ArgumentParser() ap.add_argument('--db', action='append', required=True, choices=['hg38', 'hg19']) ap.add_argument('--output-dir', required=True) args = ap.parse_args() + # Source file lives under --output-dir/cmp_downloads (see cmpVCEPProvisionalClass). + global WALSH_XLSX + WALSH_XLSX = f'{args.output_dir}/cmp_downloads/walsh/walsh2019_supplement.xlsx' + out_dir = os.path.join(args.output_dir, 'cmpVCEPWalsh2019') os.makedirs(out_dir, exist_ok=True) print(' [B.7c Walsh 2019 Pre-EvRepo curated]') records = load_walsh_table_s6() lookup = build_clinvar_lookup() from collections import Counter classes = Counter(r['classification'] for r in records) print(f' classifications: {dict(classes)}') n_upgraded = sum(1 for r in records if r['upgraded']) print(f' Walsh-upgraded (PM1 EF rule): {n_upgraded}') genes = Counter(r['gene'] for r in records) print(f' genes: {dict(genes)}')