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/cmpVCEPCardioBoost.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPCardioBoost.py index 545bec84d68..53b463733bc 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPCardioBoost.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPCardioBoost.py @@ -1,43 +1,41 @@ #!/usr/bin/env python3 """ -B.5 — CardioBoost disease-specific missense predictor track. +B.5 - CardioBoost disease-specific missense predictor track. CardioBoost (Zhang et al. 2021, Genome Medicine 13:31) is a cardiomyopathy/arrhythmia disease-specific missense pathogenicity predictor, complementary to REVEL. We render the -precomputed cardiomyopathy-model scores (cm_prediction.RData → cardioboost_8genes_hg19.tsv, +precomputed cardiomyopathy-model scores (cm_prediction.RData -> cardioboost_8genes_hg19.tsv, exported with /usr/bin/Rscript) for our 8 genes. IMPORTANT: CardioBoost is NOT a CSpec-specified predictor (the CSpec uses REVEL for PP3/BP4). -This track is informational/supplementary — it does NOT fire an ACMG code. Whether the VCEP +This track is informational/supplementary - it does NOT fire an ACMG code. Whether the VCEP wants it used is an open question (see Phase 7). The displayed CardioBoost class uses the -tool's OWN published thresholds (>=0.9 pathogenic, <=0.1 benign), cited — not invented here. +tool's OWN published thresholds (>=0.9 pathogenic, <=0.1 benign), cited - not invented here. Source is hg19 (GRCh37); we build hg19 natively and liftOver hg19->hg38. Outputs: cmpVCEPCardioBoost/cmpVCEPCardioBoost.as cmpVCEPCardioBoost/cmpVCEPCardioBoostHg{38,19}.bed + .bb """ -import os, subprocess, sys +import argparse, os, subprocess, sys -WORKDIR = '/hive/users/lrnassar/claude/RM37446' -TSV = f'{WORKDIR}/cmp_downloads/cardioboost/cardioboost_8genes_hg19.tsv' CHROM_SIZES = {'hg38': '/cluster/data/hg38/chrom.sizes', 'hg19': '/cluster/data/hg19/chrom.sizes'} LIFTOVER_HG19_TO_HG38 = '/cluster/data/hg19/bed/liftOver/hg19ToHg38.over.chain.gz' -# CardioBoost's OWN published classification thresholds (Zhang 2021) — cited, not invented. +# CardioBoost's OWN published classification thresholds (Zhang 2021) - cited, not invented. CB_PATH = 0.90 CB_BENIGN = 0.10 COLORS = {'Pathogenic': '210,0,0', 'Benign': '0,160,0', 'Indeterminate': '150,150,150'} AUTOSQL = """table cmpVCEPCardioBoost "CardioBoost disease-specific missense predictor (Zhang 2021) - informational, NOT a CSpec code" ( string chrom; "Chromosome" uint chromStart; "Position (BED 0-based)" uint chromEnd; "End (BED half-open; +1 for SNV)" string name; "Display name (gene + protein change)" uint score; "0" char[1] strand; "Strand" uint thickStart; "Same as chromStart" uint thickEnd; "Same as chromEnd" @@ -50,87 +48,95 @@ double cbScore; "CardioBoost pathogenicity probability (0-1)" string cbClass; "CardioBoost class at published thresholds (>=0.9 P, <=0.1 B)" lstring _mouseOver; "Tooltip HTML" ) """ def cb_class(score): if score >= CB_PATH: return 'Pathogenic' if score <= CB_BENIGN: return 'Benign' return 'Indeterminate' -def build_hg19_bed(): +def build_hg19_bed(tsv): rows = [] - with open(TSV) as fh: + with open(tsv) as fh: header = fh.readline() for line in fh: f = line.rstrip('\n').split('\t') if len(f) < 8: continue chrom_num, pos, ref, alt, gene, hgvsc, hgvsp, score_s = f[:8] try: pos = int(pos); score = float(score_s) except ValueError: continue if len(ref) != 1 or len(alt) != 1: # CardioBoost is missense SNVs continue chrom = 'chr' + chrom_num start, end = pos - 1, pos cls = cb_class(score) pshort = hgvsp.split(':')[-1] if ':' in hgvsp else hgvsp name = f'{gene}_{pshort}_{ref}{alt}_{score:.2f}' # ref/alt keeps same-aa variants distinct mouse = (f'<b>CardioBoost</b> - disease-specific missense predictor (Zhang 2021)<br>' f'{gene} {chrom}:{pos} {ref}>{alt} <i>{pshort}</i><br>' f'<b>CardioBoost probability:</b> {score:.3f} (class {cls}; tool thresholds ≥{CB_PATH} P, ≤{CB_BENIGN} B)<br>' f'<i>Informational — NOT a CSpec-specified predictor; the CSpec uses REVEL for ' f'PP3/BP4. Does not fire an ACMG code.</i>') rows.append('\t'.join([chrom, str(start), str(end), name, '0', '+', str(start), str(end), COLORS[cls], gene, ref, alt, hgvsc, pshort, f'{score:.6f}', cls, mouse])) rows.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) return rows def main(): - out_dir = os.path.join(WORKDIR, 'cmpVCEPCardioBoost') + 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() + + out_dir = os.path.join(args.output_dir, 'cmpVCEPCardioBoost') os.makedirs(out_dir, exist_ok=True) + tsv = os.path.join(args.output_dir, 'cmp_downloads/cardioboost/cardioboost_8genes_hg19.tsv') print(' [B.5 CardioBoost predictor]') as_path = os.path.join(out_dir, 'cmpVCEPCardioBoost.as') with open(as_path, 'w') as f: f.write(AUTOSQL) - hg19_rows = build_hg19_bed() + # hg19 bed is always built: it is the native assembly and the liftOver source for hg38. + hg19_rows = build_hg19_bed(tsv) hg19_bed = os.path.join(out_dir, 'cmpVCEPCardioBoostHg19.bed') with open(hg19_bed, 'w') as f: f.write('\n'.join(hg19_rows) + '\n') print(f' hg19 (native): {len(hg19_rows)} variants') + if 'hg19' in args.db: hg19_bb = os.path.join(out_dir, 'cmpVCEPCardioBoostHg19.bb') subprocess.run(['bedToBigBed', '-tab', '-type=bed9+8', '-as=' + as_path, hg19_bed, CHROM_SIZES['hg19'], hg19_bb], check=True) + if 'hg38' in args.db: # liftOver hg19 -> hg38 hg38_bed = os.path.join(out_dir, 'cmpVCEPCardioBoostHg38.bed') unmapped = hg38_bed + '.unmapped' subprocess.run(['liftOver', '-bedPlus=9', '-tab', hg19_bed, LIFTOVER_HG19_TO_HG38, hg38_bed, unmapped], check=True) n_un = sum(1 for line in open(unmapped) if not line.startswith('#')) if os.path.getsize(unmapped) else 0 print(f' hg38 (liftOver hg19->hg38): {sum(1 for _ in open(hg38_bed))} variants, {n_un} unmapped') - # re-sort hg38 bed after liftOver lines = [l.rstrip('\n') for l in open(hg38_bed)] lines.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) with open(hg38_bed, 'w') as f: f.write('\n'.join(lines) + '\n') hg38_bb = os.path.join(out_dir, 'cmpVCEPCardioBoostHg38.bb') subprocess.run(['bedToBigBed', '-tab', '-type=bed9+8', '-as=' + as_path, hg38_bed, CHROM_SIZES['hg38'], hg38_bb], check=True) - print(f' wrote bigBeds for both assemblies') + print(' wrote bigBeds') if __name__ == '__main__': main()