6577d5ee1319bbea85988c1c89179436b4a94edf lrnassar Tue Jul 14 11:27:59 2026 -0700 Address code-review feedback on the Cardiomyopathy VCEP build scripts. refs #37446 - cmpVCEPCardioBoost.py: add the standard --db/--output-dir CLI. It previously hardcoded the working directory for both its input TSV and its output (unlike the 11 sibling scripts, and contrary to the makedoc's documented interface); the build loop's flags were silently ignored. Output is unchanged (31,236 variants per assembly). - Decode leftover HTML entities (arrows, >=, <=, +/-, x) in print/stderr diagnostics, comments, and docstrings across all scripts so build logs read cleanly. The mouseOver / bigBed display strings intentionally keep their entities. - cmpVCEPWalsh2019.py: fix the stale docstring that described the ClinVar-unmatched entries as "deferred" (they are mapped via the hgvsToVcf fallback, item L) and drop the unverified "163 rows" count. Per code-review feedback on commit aa5669fe64. No track data changed. diff --git src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPWalsh2019.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPWalsh2019.py index 07a0487a4df..81881ac7f2a 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPWalsh2019.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPWalsh2019.py @@ -1,44 +1,44 @@ #!/usr/bin/env python3 """ -B.7c — Walsh 2019 Pre-EvRepo curated variants subtrack (folds into Curated Variants composite). +B.7c - Walsh 2019 Pre-EvRepo curated variants subtrack (folds into Curated Variants composite). -Renders 155 per-variant ACMG/AMP rule applications from Walsh 2019 Table S6 — pre-EvRepo +Renders 155 per-variant ACMG/AMP rule applications from Walsh 2019 Table S6 - pre-EvRepo VCEP curations that document the original calibration cohort. Folds into the Curated Variants composite track 4c, off by default. -Source: cmp_downloads/walsh/walsh2019_supplement.xlsx Table S6 (163 rows; filter to our 8 genes) -Coords: lookup against ClinVar variant_summary by (gene, c.notation); skip variants not in ClinVar - (those are novel-to-Walsh and would need MANE CDS conversion — deferred). +Source: cmp_downloads/walsh/walsh2019_supplement.xlsx Table S6, filtered to our 8 genes. +Coords: lookup against ClinVar variant_summary by (gene, c.notation); entries not in ClinVar + are mapped via hgvsToVcf on each gene's Walsh transcript (item L), so all 155 render. Outputs: cmpVCEPWalsh2019/cmpVCEPWalsh2019.as cmpVCEPWalsh2019/cmpVCEPWalsh2019Hg{38,19}.bed + .bb Usage: python3 cmpVCEPWalsh2019.py --db hg38 --db hg19 \ --output-dir /hive/users/lrnassar/claude/RM37446 """ import argparse, gzip, os, re, subprocess, sys, warnings warnings.filterwarnings('ignore') WALSH_XLSX = '/hive/users/lrnassar/claude/RM37446/cmp_downloads/walsh/walsh2019_supplement.xlsx' VARIANT_SUMMARY = '/hive/data/outside/otto/clinvar/downloads/2026-05-30/variant_summary.txt.gz' # Transcript Walsh 2019 used for c. numbering, per gene (item L: map ClinVar-unmatched entries -# via hgvsToVcf). NOTE TNNT2 is NOT MANE — Walsh used the classic cardiac transcript; MANE +# via hgvsToVcf). NOTE TNNT2 is NOT MANE - Walsh used the classic cardiac transcript; MANE # (NM_001276345.2) yields HgvsRefAssertedMismatch. Verified each gives FILTER=PASS. HGVSTOVCF = '/cluster/bin/x86_64/hgvsToVcf' WALSH_TX = {'MYBPC3': 'NM_000256.3', 'MYH7': 'NM_000257.4', 'TNNI3': 'NM_000363.5', 'TNNT2': 'NM_001001430.2'} LIFTOVER_HG38_HG19 = '/cluster/data/hg38/bed/liftOver/hg38ToHg19.over.chain.gz' def _tool_coords(gene, cdna, db): """hgvsToVcf on the gene's Walsh transcript; coords only if FILTER==PASS.""" tx = WALSH_TX.get(gene) if not tx: return None try: @@ -134,31 +134,31 @@ string variantProtein; "p. notation from Walsh 2019 Table S6" string variantType; "missense, nonsense, frameshift, splice site, etc." string classification; "VCEP classification (Pathogenic, Likely Pathogenic, VUS, etc.)" string walshUpgraded; "yes if Walsh upgraded via new EF-based PM1 rule" string acmgRules; "ACMG/AMP rules activated (e.g., PM2,PP3,PM1(s))" string cases; "Number of cases observed in Walsh 2019 cohort" string exacFaf; "ExAC filtering allele frequency" string source; "Source for curated evidence (PMID or ClinVar SCV, where applicable)" string variationId; "Matched ClinVar VariationID (if found)" lstring _mouseOver; "Tooltip HTML" ) """ def load_walsh_table_s6(): - """Parse Walsh 2019 Table S6 → list of dict records for our 8 genes.""" + """Parse Walsh 2019 Table S6 -> list of dict records for our 8 genes.""" import openpyxl wb = openpyxl.load_workbook(WALSH_XLSX, read_only=True, data_only=True) ws = wb['Table S6'] rows = list(ws.iter_rows(values_only=True)) header_idx = None for i, r in enumerate(rows): if r and r[0] == 'Gene': header_idx = i break if header_idx is None: sys.exit('Table S6: header not found') records = [] for r in rows[header_idx + 1:]: @@ -179,31 +179,31 @@ 'classification': cls, 'upgraded': upgraded, 'rules': r[6] or '', 'cases': r[7] if r[7] is not None else '', 'source': r[8] if r[8] is not None else '', }) print(f' parsed {len(records)} Walsh 2019 Table S6 entries (filtered to our 8 genes)') return records # Matches ClinVar Name field: e.g. NM_000256.3(MYBPC3):c.1504C>T (p.Arg502Trp) CV_NAME_RE = re.compile(r'^([A-Z]M_[\d\.]+)\(([A-Z0-9]+)\):c\.(\S+?)(?:\s|\(|$)') def build_clinvar_lookup(): - """Stream variant_summary; build dict (gene, c.notation) → {assembly: coords + variation_id}.""" + """Stream variant_summary; build dict (gene, c.notation) -> {assembly: coords + variation_id}.""" lookup = {} 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: @@ -215,31 +215,31 @@ 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, 'variation_id': f[30], 'rcv': f[11], } n_rows += 1 - print(f' ClinVar lookup: {len(lookup)} (gene, c.notation) keys × up to 2 assemblies = {n_rows} entries') + print(f' ClinVar lookup: {len(lookup)} (gene, c.notation) keys x up to 2 assemblies = {n_rows} entries') return lookup def emit_bed(records, lookup, db, out_path): rows_emitted = [] skipped = [] for rec in records: key = (rec['gene'], rec['cdna']) c = lookup.get(key, {}).get(db) via_tool = False if c is None: c = walsh_coords_via_tool(rec['gene'], rec['cdna'], db) via_tool = c is not None if c is None: skipped.append(rec) @@ -273,31 +273,31 @@ rec['vartype'], cls, upgraded_label, rec['rules'], str(rec['cases']), str(rec['exac_faf']), str(rec['source']), c['variation_id'], mouseover, ])) rows_emitted.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) with open(out_path, 'w') as f: for line in rows_emitted: f.write(line + '\n') - print(f' wrote {len(rows_emitted)} BED features → {out_path} (skipped {len(skipped)})') + print(f' wrote {len(rows_emitted)} BED features -> {out_path} (skipped {len(skipped)})') return len(rows_emitted), skipped 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() @@ -325,30 +325,30 @@ counts = {} skipped_summary = None for db in args.db: bed_path = os.path.join(out_dir, f'cmpVCEPWalsh2019Hg{"38" if db=="hg38" else "19"}.bed') bb_path = os.path.join(out_dir, f'cmpVCEPWalsh2019Hg{"38" if db=="hg38" else "19"}.bb') n, skipped = emit_bed(records, lookup, db, bed_path) make_bigbed(bed_path, db, as_path, bb_path) counts[db] = n skipped_summary = skipped print(f' {db} bigBed: {bb_path}') if 'hg38' in counts and 'hg19' in counts: if counts['hg38'] == counts['hg19']: print(f' cross-assembly parity OK: {counts["hg38"]} features each') else: - print(f' WARNING: parity FAILED — hg38={counts["hg38"]} hg19={counts["hg19"]}', file=sys.stderr) + print(f' WARNING: parity FAILED - hg38={counts["hg38"]} hg19={counts["hg19"]}', file=sys.stderr) # Save skipped list for follow-up if skipped_summary: skip_path = os.path.join(out_dir, 'walsh2019_unmatched.txt') with open(skip_path, 'w') as f: - f.write('# Walsh 2019 Table S6 entries NOT found in ClinVar variant_summary\n') - f.write('# These need MANE CDS coordinate conversion to render — deferred to v2 of B.7c\n') + f.write('# Walsh 2019 Table S6 entries not placed via ClinVar or the hgvsToVcf fallback\n') + f.write('# (expected to be empty; any listed here could not be mapped to genomic coords)\n') for r in skipped_summary: f.write(f'{r["gene"]}\t{r["cdna"]}\t{r["protein"]}\t{r["classification"]}\n') print(f' skipped variants logged to: {skip_path}') if __name__ == '__main__': main()