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/cmpVCEPAtlasEF.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAtlasEF.py index 73897ceda9a..837c066e27f 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAtlasEF.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAtlasEF.py @@ -1,50 +1,50 @@ #!/usr/bin/env python3 """ -B.10 — Atlas of Cardiac Genetic Variation per-variant track (PS4 case-control). +B.10 - Atlas of Cardiac Genetic Variation per-variant track (PS4 case-control). Parses 189 cached per-variant HTML pages from the Atlas scrape (A.9), extracts -case carriers + cohort + ExAC data, computes Fisher 2×2 OR + 95% CI with +case carriers + cohort + ExAC data, computes Fisher 2x2 OR + 95% CI with Haldane 0.5 correction, and renders as bigBed. Color by computed CI_lower PS4 strength binning: STRONG CI_lower >= 20 MODERATE CI_lower >= 10 SUPPORTING CI_lower >= 5 below CI_lower < 5 Mouseover shows: case counts, ExAC carriers, computed OR + CI, etiologic fraction (from Atlas), OMGL/LMM classification. -Validation: spot-check 5 known PS4-relevant variants — recompute OR and verify -within ±5% of cardiodb.org/cmgwap/ output (deferred to D.1 audit step). +Validation: spot-check 5 known PS4-relevant variants - recompute OR and verify +within +/-5% of cardiodb.org/cmgwap/ output (deferred to D.1 audit step). Source: cmp_downloads/atlas/variants/var_*.html (cached scrape from A.9) Outputs: cmpVCEPAtlasEF/cmpVCEPAtlasEF.as cmpVCEPAtlasEF/cmpVCEPAtlasEFHg{38,19}.bed + .bb """ import argparse, math, os, re, subprocess, sys from html.parser import HTMLParser ATLAS_VAR_DIR = '/hive/users/lrnassar/claude/RM37446/cmp_downloads/atlas/variants' OUR_GENES = {'MYH7', 'MYBPC3', 'TNNT2', 'TNNI3', 'TPM1', 'ACTC1', 'MYL2', 'MYL3'} -# ExAC reference cohort size — DEFAULT (Walsh 2017 baseline; used as fallback) +# ExAC reference cohort size - DEFAULT (Walsh 2017 baseline; used as fallback) EXAC_TOTAL_INDIVIDUALS = 60706 # Per-gene Non-Truncating ExAC denominators from Walsh 2019 Table S1. # Used for non-truncating Atlas variants (missense, synonymous, inframe) to match # the cohort sizes the VCEP CSpec PS4 calibration is grounded in. # Truncating/splice variants fall back to the default; documented as a v1 limitation. WALSH_2019_EXAC_NONTRUNC = { 'MYH7': 60469, 'MYBPC3': 45794, 'TNNT2': 57018, 'TNNI3': 52607, 'TPM1': 58642, 'MYL2': 60521, 'MYL3': 60605, 'ACTC1': 60198, @@ -232,31 +232,31 @@ se = math.sqrt(1/a + 1/b + 1/c + 1/d) ci_lo = math.exp(log_or - 1.96 * se) ci_hi = math.exp(log_or + 1.96 * se) return or_val, ci_lo, ci_hi def ps4_bin(ci_lower): if ci_lower is None: return 'unable to compute', PS4_COLOR_BELOW if ci_lower >= 20: return 'Strong', PS4_COLOR_STRONG if ci_lower >= 10: return 'Moderate', PS4_COLOR_MODERATE if ci_lower >= 5: return 'Supporting', PS4_COLOR_SUPPORTING return 'below threshold', PS4_COLOR_BELOW def liftover_grch37_to_grch38(grch37_records, out_dir): - """liftOver hg19 → hg38. Returns dict of grch37_pos → (chrom, start, end) in hg38.""" + """liftOver hg19 -> hg38. Returns dict of grch37_pos -> (chrom, start, end) in hg38.""" if not grch37_records: return {} tmp_bed = os.path.join(out_dir, '_atlas_lift.bed') out_bed = os.path.join(out_dir, '_atlas_lift_out.bed') unmapped = os.path.join(out_dir, '_atlas_lift.unmapped') with open(tmp_bed, 'w') as f: for i, r in enumerate(grch37_records): # 1-based to BED 0-based start = r['grch37_pos'] - 1 end = start + 1 f.write(f'{r["grch37_chrom"]}\t{start}\t{end}\tatlas_{r["var_id"]}\t0\t+\n') cmd = ['liftOver', tmp_bed, LIFTOVER_HG19_TO_HG38, out_bed, unmapped] subprocess.run(cmd, check=True) mapping = {} for line in open(out_bed): @@ -285,59 +285,59 @@ continue var_id = int(fname.replace('var_', '').replace('.html', '')) html = open(os.path.join(ATLAS_VAR_DIR, fname)).read() rec = parse_variant_page(html, var_id) if rec is None: parse_failures += 1 continue records.append(rec) print(f' parsed {len(records)} variant pages; {parse_failures} parse failures') from collections import Counter print(f' gene distribution: {Counter(r["gene"] for r in records)}') # liftOver Atlas GRCh37 coords to hg38 lift_map = liftover_grch37_to_grch38(records, out_dir) - print(f' liftOver hg19→hg38: {len(lift_map)} of {len(records)} mapped') + print(f' liftOver hg19->hg38: {len(lift_map)} of {len(records)} mapped') # Build BED features bed_lines = [] n_strong = n_moderate = n_supporting = n_below = 0 for r in records: # Compute combined case carriers + cohort across HCM and DCM # Use the disease that has data; if both, sum diseases_with_data = [] for disease in ('HCM', 'DCM'): cd = r['case_data'][disease] if cd['omgl'] is not None or cd['lmm'] is not None: diseases_with_data.append(disease) if not diseases_with_data: - continue # no case data → skip + continue # no case data -> skip # For each disease, compute OR for disease in diseases_with_data: cd = r['case_data'][disease] omgl = cd['omgl'] or (0, 0) lmm = cd['lmm'] or (0, 0) case_carriers = omgl[0] + lmm[0] case_cohort = omgl[1] + lmm[1] if case_cohort == 0 or case_carriers == 0: continue # skip zero-case rows - # ExAC carriers — Atlas reports allele count, not individual count. - # Use exac_carriers as "individuals carrying" (approximation: rare variants ≈ heterozygous only) + # ExAC carriers - Atlas reports allele count, not individual count. + # Use exac_carriers as "individuals carrying" (approximation: rare variants ~ heterozygous only) exac_c = r['exac_carriers'] or 0 # Per D.1 audit P0 #1: use Walsh 2019 Table S1 per-gene NonTrunc denominators # for non-truncating variants; truncating/splice fall back to default. vartype_lower = (r['vartype'] or '').lower() if vartype_lower in NONTRUNC_VARTYPES: exac_cohort = WALSH_2019_EXAC_NONTRUNC.get(r['gene'], EXAC_TOTAL_INDIVIDUALS) else: exac_cohort = EXAC_TOTAL_INDIVIDUALS a = case_carriers b = case_cohort - case_carriers c = exac_c d = exac_cohort - exac_c or_val, ci_lo, ci_hi = fisher_or_haldane(a, b, c, d) strength, color = ps4_bin(ci_lo) @@ -388,73 +388,73 @@ mouseover, ])) print(f' PS4 binning: STRONG={n_strong}, MODERATE={n_moderate}, SUPPORTING={n_supporting}, below={n_below}') bed_lines.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) as_path = os.path.join(out_dir, 'cmpVCEPAtlasEF.as') with open(as_path, 'w') as f: f.write(AUTOSQL) hg38_bed = os.path.join(out_dir, 'cmpVCEPAtlasEFHg38.bed') with open(hg38_bed, 'w') as f: for l in bed_lines: f.write(l + '\n') - print(f' wrote {len(bed_lines)} BED features → {hg38_bed}') + print(f' wrote {len(bed_lines)} BED features -> {hg38_bed}') if 'hg38' in args.db: hg38_bb = os.path.join(out_dir, 'cmpVCEPAtlasEFHg38.bb') cmd = ['bedToBigBed', '-tab', '-type=bed9+16', '-as=' + as_path, hg38_bed, CHROM_SIZES['hg38'], hg38_bb] print(f' $ {" ".join(cmd)}') subprocess.run(cmd, check=True) # hg19: derive from grch37 records directly (no liftover round-trip) if 'hg19' in args.db: hg19_bed = os.path.join(out_dir, 'cmpVCEPAtlasEFHg19.bed') # Re-emit lines using GRCh37 coords from original parse - # Build a map var_id → row, then re-emit with GRCh37 coords + # Build a map var_id -> row, then re-emit with GRCh37 coords var_recs = {r['var_id']: r for r in records} hg19_lines = [] for line in bed_lines: f = line.split('\t') # the first BED9 are hg38 coords; everything else is the same # Find the var_id from atlas URL field (last data column before mouseover) atlas_url = f[-2] # '...atlas_url' is the second-to-last (mouseover is last) m = re.search(r'id=(\d+)', atlas_url) if not m: continue var_id = int(m.group(1)) r = var_recs.get(var_id) if r is None: continue - # GRCh37 coords (1-based) → BED 0-based half-open (assume 1-bp) + # GRCh37 coords (1-based) -> BED 0-based half-open (assume 1-bp) grch37_start = r['grch37_pos'] - 1 grch37_end = grch37_start + 1 f[0] = r['grch37_chrom'] f[1] = str(grch37_start) f[2] = str(grch37_end) f[6] = str(grch37_start) f[7] = str(grch37_end) hg19_lines.append('\t'.join(f)) hg19_lines.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) with open(hg19_bed, 'w') as fh: for l in hg19_lines: fh.write(l + '\n') hg19_bb = os.path.join(out_dir, 'cmpVCEPAtlasEFHg19.bb') cmd = ['bedToBigBed', '-tab', '-type=bed9+16', '-as=' + as_path, hg19_bed, CHROM_SIZES['hg19'], hg19_bb] print(f' $ {" ".join(cmd)}') subprocess.run(cmd, check=True) if 'hg38' in args.db and 'hg19' in args.db: n38 = len(bed_lines) n19 = len(hg19_lines) if n38 == n19: print(f' cross-assembly parity OK: {n38} features each') else: - print(f' WARNING: parity FAILED — hg38={n38} hg19={n19}', file=sys.stderr) + print(f' WARNING: parity FAILED - hg38={n38} hg19={n19}', file=sys.stderr) if __name__ == '__main__': main()