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/cmpVCEPClinVar506161.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinVar506161.py index 116125b3ef9..d41eaa2bed9 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinVar506161.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinVar506161.py @@ -1,21 +1,21 @@ #!/usr/bin/env python3 """ -B.8 — Cardiomyopathy VCEP ClinVar submitter 506161 track builder. +B.8 - Cardiomyopathy VCEP ClinVar submitter 506161 track builder. Renders all 199 ClinVar submissions from "ClinGen Cardiomyopathy Variant Curation -Expert Panel" — verified at A.3 to be 100% "reviewed by expert panel" (no relabel +Expert Panel" - verified at A.3 to be 100% "reviewed by expert panel" (no relabel branch needed). Input: /hive/users/lrnassar/claude/RM37446/cmp_downloads/clinvar/cardiomyopathyVCEP_submissions.tsv (199 submission rows, schema = ClinVar submission_summary.txt.gz) /hive/data/outside/otto/clinvar/downloads/2026-05-30/variant_summary.txt.gz (per-VariationID hg19 + hg38 coords; joined by col 31 = VariationID) Outputs: cmpVCEPClinVar506161/cmpVCEPClinVar506161.as cmpVCEPClinVar506161/cmpVCEPClinVar506161Hg38.bed + .bb cmpVCEPClinVar506161/cmpVCEPClinVar506161Hg19.bed + .bb Usage: python3 cmpVCEPClinVar506161.py --db hg38 --db hg19 \ @@ -26,31 +26,31 @@ SUBMISSIONS_TSV = '/hive/users/lrnassar/claude/RM37446/cmp_downloads/clinvar/cardiomyopathyVCEP_submissions.tsv' VARIANT_SUMMARY = '/hive/data/outside/otto/clinvar/downloads/2026-05-30/variant_summary.txt.gz' # Standard 5-tier ACMG ramp (same as EvRepo / Walsh 2019); the "with codes vs no codes" # distinction is carried by the track name and mouseover, not by color. COLORS = { 'Pathogenic': '210,0,0', 'Likely pathogenic': '245,152,152', 'Uncertain significance': '0,0,136', 'Likely benign': '213,247,213', 'Benign': '0,210,0', } DEFAULT_COLOR = '136,136,136' -# Per-gene EvRepo VariationIDs — built once at build-time from EvRepo JSON +# Per-gene EvRepo VariationIDs - built once at build-time from EvRepo JSON EVREPO_VAR_IDS = set() EVREPO_JSON = '/hive/users/lrnassar/claude/RM37446/cmp_downloads/erepo/cardiomyopathyVCEP_classifications.json' try: _erepo = json.load(open(EVREPO_JSON)) EVREPO_VAR_IDS = set(v.get('variationId', '') for v in _erepo['variantInterpretations'] if v.get('variationId')) print(f' loaded {len(EVREPO_VAR_IDS)} EvRepo VariationIDs (will tag overlapping rows)') except Exception as e: print(f' WARNING: could not load EvRepo JSON: {e}', file=sys.stderr) CHROM_SIZES = { 'hg38': '/cluster/data/hg38/chrom.sizes', 'hg19': '/cluster/data/hg19/chrom.sizes', } AUTOSQL = """table cmpVCEPClinVar506161 @@ -69,75 +69,75 @@ string hgvsName; "Canonical HGVS Name (NM_xx(GENE):c.x>y (p.xxx))" string classification; "VCEP classification" string reviewStatus; "ClinVar review status (verified all-expert-panel at A.3)" string lastEvaluated; "Date last evaluated" string scvAccession; "SCV submission accession" string variationId; "ClinVar VariationID" string rcvAccession; "ClinVar RCV accession" string inEvRepo; "yes if also in EvRepo (variant is in both sources)" string clinvarUrl; "ClinVar variant URL" lstring _mouseOver; "Tooltip HTML" ) """ def load_submissions(): - """Parse 199-row TSV. Returns dict by VariationID (string) → submission record.""" + """Parse 199-row TSV. Returns dict by VariationID (string) -> submission record.""" submissions = {} header = None for line in open(SUBMISSIONS_TSV): line = line.rstrip('\n') if line.startswith('##'): continue if line.startswith('#'): header = line.lstrip('#').split('\t') continue f = line.split('\t') if len(f) < 12: continue var_id = f[0] submissions[var_id] = { 'classification': f[1], 'last_evaluated': f[2], 'description': f[3], 'review_status': f[6], 'submitter': f[9], 'scv': f[10], 'gene': f[11] if len(f) > 11 else '', } print(f' loaded {len(submissions)} ClinVar submitter-506161 records') return submissions def load_variant_summary_coords(submissions): """Stream variant_summary.txt.gz, picking rows whose VariationID (col 31) is in our set. - Returns dict: var_id → {assembly: {chrom, start_bed, end_bed, ref, alt, hgvs_name, rcv}}.""" + Returns dict: var_id -> {assembly: {chrom, start_bed, end_bed, ref, alt, hgvs_name, rcv}}.""" var_ids = set(submissions.keys()) coords = {vid: {} for vid in var_ids} # Schema (1-based positions in awk; 0-based here): # 3=Name, 5=GeneSymbol, 12=RCVaccession, 17=Assembly, 18=ChromosomeAccession, # 19=Chromosome, 20=Start, 21=Stop, 22=ReferenceAllele, 23=AlternateAllele, 31=VariationID n_rows_seen = 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 - var_id = f[30] # VariationID (1-based col 31 → 0-based index 30) + var_id = f[30] # VariationID (1-based col 31 -> 0-based index 30) if var_id not in var_ids: continue n_rows_seen += 1 assembly = f[16] # GRCh38 / GRCh37 chrom_num = f[18] # 1, 2, ..., X try: start1 = int(f[19]) # 1-based start stop1 = int(f[20]) # 1-based stop (inclusive) except ValueError: continue ref = f[21] alt = f[22] chrom = f'chr{chrom_num}' db = 'hg38' if assembly == 'GRCh38' else ('hg19' if assembly == 'GRCh37' else None) if db is None: @@ -204,31 +204,31 @@ cls, sub['review_status'], sub['last_evaluated'], sub['scv'], var_id, c['rcv'], in_evrepo, clinvar_url, mouseover, ])) rows.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) with open(out_path, 'w') as f: for line in rows: f.write(line + '\n') - print(f' wrote {len(rows)} BED features → {out_path} (skipped {skipped})') + print(f' wrote {len(rows)} BED features -> {out_path} (skipped {skipped})') return len(rows) def make_bigbed(bed_path, db, as_path, bb_path): cmd = ['bedToBigBed', '-tab', '-type=bed9+11', '-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() @@ -255,21 +255,21 @@ f.write(AUTOSQL) counts = {} for db in args.db: bed_path = os.path.join(out_dir, f'cmpVCEPClinVar506161Hg{"38" if db=="hg38" else "19"}.bed') bb_path = os.path.join(out_dir, f'cmpVCEPClinVar506161Hg{"38" if db=="hg38" else "19"}.bb') n = emit_bed(submissions, coords, db, bed_path) make_bigbed(bed_path, db, as_path, bb_path) counts[db] = n 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"]}', + print(f' WARNING: parity FAILED - hg38={counts["hg38"]} hg19={counts["hg19"]}', file=sys.stderr) if __name__ == '__main__': main()