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/cmpVCEPClinDomains.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinDomains.py index aab352f9a8a..1da7a393ed8 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinDomains.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinDomains.py @@ -1,57 +1,57 @@ #!/usr/bin/env python3 """ -B.1 — Cardiomyopathy VCEP PM1 Hotspot Regions track builder. +B.1 - Cardiomyopathy VCEP PM1 Hotspot Regions track builder. Renders PM1 hotspot codon regions per the per-gene CSpec for the 4 of 8 genes where PM1 is applicable (MYH7, MYBPC3, TNNT2, TNNI3). PM1 is NOT specified for ACTC1, MYBPC3*, MYL2, MYL3 (where * = MYBPC3 has PM1 applicable but only on missense). PM1 ranges (per CSpec, see ../doc/Cardiomyopathy.txt §A.10): - MYH7 167–931 (NM_000257.4 / Walsh 2019 calibration) - MYBPC3 485–502 + 1248–1266 (NM_000256.3) - TNNT2 89–189 (NM_001276345.2) - TNNI3 141–209 (NM_000363.5) + MYH7 167-931 (NM_000257.4 / Walsh 2019 calibration) + MYBPC3 485-502 + 1248-1266 (NM_000256.3) + TNNT2 89-189 (NM_001276345.2) + TNNI3 141-209 (NM_000363.5) Reads MANE Select bigGenePred from /gbdb/hg38/mane/mane.bb to extract CDS exon structure; converts amino-acid ranges to genomic coordinates with explicit unit tests against the first/last codon of each gene (catches the InSiGHT PMS2 off-by-one bug class). Outputs (per --output-dir): cmpVCEPClinDomains/cmpVCEPClinDomains.as cmpVCEPClinDomains/cmpVCEPClinDomainsHg38.bed cmpVCEPClinDomains/cmpVCEPClinDomainsHg38.bb cmpVCEPClinDomains/cmpVCEPClinDomainsHg19.bed (via liftOver) cmpVCEPClinDomains/cmpVCEPClinDomainsHg19.bb Usage: python3 cmpVCEPClinDomains.py --db hg38 --db hg19 \ --output-dir /hive/users/lrnassar/claude/RM37446 """ import argparse, os, subprocess, sys, tempfile # PM1 hotspot regions (gene, transcript, [(aa_start, aa_end, label), ...]) PM1_REGIONS = { 'MYH7': ('NM_000257.4', [(167, 931, 'Head/neck/converter')]), 'MYBPC3': ('NM_000256.3', [(485, 502, 'Hotspot 1'), (1248, 1266, 'Hotspot 2')]), 'TNNT2': ('NM_001276345.2', [(89, 189, 'Hotspot')]), 'TNNI3': ('NM_000363.5', [(141, 209, 'Hotspot')]), } -PM1_COLOR = '230,3,131' # magenta-rose — matches InSiGHT clinDomains + TP53 clinical-domains convention +PM1_COLOR = '230,3,131' # magenta-rose - matches InSiGHT clinDomains + TP53 clinical-domains convention MANE_BB = '/gbdb/hg38/mane/mane.bb' LIFTOVER_HG38_TO_HG19 = '/cluster/data/hg38/bed/liftOver/hg38ToHg19.over.chain.gz' CHROM_SIZES = { 'hg38': '/cluster/data/hg38/chrom.sizes', 'hg19': '/cluster/data/hg19/chrom.sizes', } AUTOSQL = """table cmpVCEPClinDomains "Cardiomyopathy VCEP PM1 hotspot regions per ClinGen CSpec" ( string chrom; "Chromosome" uint chromStart; "Start position" uint chromEnd; "End position" string name; "Display name (gene + region label)" @@ -85,31 +85,31 @@ 'chromStart': int(f[1]), 'chromEnd': int(f[2]), 'name': f[3], 'strand': f[5], 'thickStart': int(f[6]), 'thickEnd': int(f[7]), 'blockSizes': [int(s) for s in f[10].rstrip(',').split(',')], 'chromStarts': [int(s) for s in f[11].rstrip(',').split(',')], 'refSeqAcc': f[21], } raise KeyError(f"Gene {gene_symbol} not found in {MANE_BB}") def cds_exons(mane): """Return list of (genomic_start, genomic_end) for CDS portions of each exon, - in genomic order (low → high coord). BED half-open semantics.""" + in genomic order (low -> high coord). BED half-open semantics.""" exons = [] for size, rstart in zip(mane['blockSizes'], mane['chromStarts']): es = mane['chromStart'] + rstart ee = es + size # clip to CDS if ee <= mane['thickStart'] or es >= mane['thickEnd']: continue exons.append((max(es, mane['thickStart']), min(ee, mane['thickEnd']))) return exons def aa_to_genomic_segments(aa_start, aa_end, mane): """Convert 1-based inclusive amino-acid range [aa_start, aa_end] to genomic BED half-open intervals. @@ -254,31 +254,31 @@ str(seg_start), str(seg_end), PM1_COLOR, gene, transcript, aa_range, exon_info, mouseover, ])) # sort by chrom + start bed_lines.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) with open(output_path, 'w') as f: for line in bed_lines: f.write(line + '\n') - print(f' wrote {len(bed_lines)} BED features → {output_path}') + print(f' wrote {len(bed_lines)} BED features -> {output_path}') return len(bed_lines) def make_bigbed(bed_path, db, as_path, bb_path): """Convert BED to bigBed using the autoSql schema.""" cmd = ['bedToBigBed', '-tab', '-type=bed9+5', '-as=' + as_path, bed_path, CHROM_SIZES[db], bb_path] print(f' $ {" ".join(cmd)}') subprocess.run(cmd, check=True) def liftover_to_hg19(hg38_bed, hg19_bed): """liftOver hg38 BED to hg19; warn on dropped features.""" unmapped = hg19_bed + '.unmapped' # Use bedExtraFieldsToBigBed style: liftOver supports -bedPlus to keep extra fields @@ -288,52 +288,52 @@ if os.path.getsize(unmapped) > 0: n = sum(1 for l in open(unmapped) if not l.startswith('#')) print(f' WARNING: {n} features failed liftOver to hg19; see {unmapped}', file=sys.stderr) def main(): ap = argparse.ArgumentParser() ap.add_argument('--db', action='append', required=True, choices=['hg38', 'hg19'], help='Assembly to build (specify --db hg38 --db hg19 for both)') ap.add_argument('--output-dir', required=True) args = ap.parse_args() out_dir = os.path.join(args.output_dir, 'cmpVCEPClinDomains') os.makedirs(out_dir, exist_ok=True) - # Run unit tests first — must pass before emitting any features + # Run unit tests first - must pass before emitting any features unit_test_codon_conversion() # Write autoSql schema as_path = os.path.join(out_dir, 'cmpVCEPClinDomains.as') with open(as_path, 'w') as f: f.write(AUTOSQL) # Build hg38 first (native), then hg19 via liftOver hg38_bed = os.path.join(out_dir, 'cmpVCEPClinDomainsHg38.bed') hg38_bb = os.path.join(out_dir, 'cmpVCEPClinDomainsHg38.bb') n_features = emit_bed(hg38_bed) if 'hg38' in args.db: make_bigbed(hg38_bed, 'hg38', as_path, hg38_bb) print(f' hg38 bigBed: {hg38_bb}') if 'hg19' in args.db: hg19_bed = os.path.join(out_dir, 'cmpVCEPClinDomainsHg19.bed') hg19_bb = os.path.join(out_dir, 'cmpVCEPClinDomainsHg19.bb') liftover_to_hg19(hg38_bed, hg19_bed) make_bigbed(hg19_bed, 'hg19', as_path, hg19_bb) print(f' hg19 bigBed: {hg19_bb}') # Cross-assembly parity if 'hg38' in args.db and 'hg19' in args.db: n_hg38 = sum(1 for _ in open(hg38_bed)) n_hg19 = sum(1 for _ in open(os.path.join(out_dir, 'cmpVCEPClinDomainsHg19.bed'))) if n_hg38 != n_hg19: - print(f' WARNING: cross-assembly parity FAILED — hg38={n_hg38} hg19={n_hg19}', + print(f' WARNING: cross-assembly parity FAILED - hg38={n_hg38} hg19={n_hg19}', file=sys.stderr) else: print(f' cross-assembly parity OK: {n_hg38} features each') if __name__ == '__main__': main()