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/cmpVCEPAFfrequencies.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAFfrequencies.py index fcd524866fa..7be77b8538b 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAFfrequencies.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAFfrequencies.py @@ -1,60 +1,60 @@ #!/usr/bin/env python3 """ -B.3 — gnomAD v4.1 Allele Frequencies track builder. +B.3 - gnomAD v4.1 Allele Frequencies track builder. -For each variant in the 8 cardiomyopathy gene coding regions (±20 nt splice padding), +For each variant in the 8 cardiomyopathy gene coding regions (+/-20 nt splice padding), parse gnomAD v4.1 exomes for the per-variant FAF95 (filtering allele frequency, 95% CI lower bound, popmax) and apply per-gene CSpec thresholds: BA1 if FAF95 >= 0.001 (all 8 genes) - BS1 if FAF95 >= 0.0001 (≥0.0002 for MYBPC3 only) + BS1 if FAF95 >= 0.0001 (>=0.0002 for MYBPC3 only) PM2_supporting if FAF95 <= 0.00004 no-code otherwise (still emitted, useful baseline) Source: /hive/data/outside/gnomAD.4/v4.1/exomes/gnomad.exomes.v4.1.sites.chr{N}.vcf.bgz Field: fafmax_faf95_max (max FAF95 across genetic ancestry groups) Outputs: cmpVCEPAFfrequencies/cmpVCEPAFfrequencies.as cmpVCEPAFfrequencies/cmpVCEPAFfrequenciesHg{38,19}.bed + .bb """ import argparse, os, re, subprocess, sys OUR_GENES = ['MYH7', 'MYBPC3', 'TNNT2', 'TNNI3', 'TPM1', 'ACTC1', 'MYL2', 'MYL3'] -# Gene-specific BS1 threshold (per CSpec — MYBPC3 outlier) +# Gene-specific BS1 threshold (per CSpec - MYBPC3 outlier) BS1_THRESHOLDS = { 'MYBPC3': 0.0002, } DEFAULT_BS1 = 0.0001 BA1_THRESHOLD = 0.001 PM2_SUPPORTING_THRESHOLD = 0.00004 SPLICE_PADDING = 20 # nt up/downstream of CDS exons for splice-region inclusion GNOMAD_VCF_PATTERN = '/hive/data/outside/gnomAD.4/v4.1/exomes/gnomad.exomes.v4.1.sites.{chrom}.vcf.bgz' TABIX = '/cluster/bin/x86_64/tabix' # Colors (matching plan) COLORS = { - 'BA1': '0,160,0', # dark green — strong benign frequency - 'BS1': '120,200,120', # light green — benign frequency - 'PM2_supporting': '250,160,160', # salmon — rarity (weak pathogenic). Salmon (not the old + 'BA1': '0,160,0', # dark green - strong benign frequency + 'BS1': '120,200,120', # light green - benign frequency + 'PM2_supporting': '250,160,160', # salmon - rarity (weak pathogenic). Salmon (not the old # fuchsia) keeps AF distinct; PM1 regions use magenta-rose 230,3,131. - 'no-code': '180,180,180', # light gray — no AF-based code + 'no-code': '180,180,180', # light gray - no AF-based code } CHROM_SIZES = { 'hg38': '/cluster/data/hg38/chrom.sizes', 'hg19': '/cluster/data/hg19/chrom.sizes', } LIFTOVER_HG38_TO_HG19 = '/cluster/data/hg38/bed/liftOver/hg38ToHg19.over.chain.gz' AUTOSQL = """table cmpVCEPAFfrequencies "gnomAD v4.1 allele frequencies in cardiomyopathy gene coding regions, with applied ACMG codes" ( string chrom; "Chromosome" uint chromStart; "Position (BED 0-based)" uint chromEnd; "End" @@ -80,34 +80,34 @@ sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from cmpVCEPClinDomains import parse_mane_record, cds_exons def parse_vcf_info(info_str): """Parse VCF INFO field into dict. Returns dict[key] = first-value-as-string.""" result = {} for kv in info_str.split(';'): if '=' in kv: k, v = kv.split('=', 1) result[k] = v return result def fetch_gene_variants(gene, mane): - """tabix-query gnomAD VCF for gene CDS region (±SPLICE_PADDING). Returns list of variant dicts.""" + """tabix-query gnomAD VCF for gene CDS region (+/-SPLICE_PADDING). Returns list of variant dicts.""" chrom = mane['chrom'] exons = cds_exons(mane) - # Build region list: each CDS exon ± SPLICE_PADDING + # Build region list: each CDS exon +/- SPLICE_PADDING regions = [] for ex_start, ex_end in exons: regions.append((max(0, ex_start - SPLICE_PADDING), ex_end + SPLICE_PADDING)) # Merge overlapping regions regions.sort() merged = [] for s, e in regions: if merged and s <= merged[-1][1]: merged[-1] = (merged[-1][0], max(merged[-1][1], e)) else: merged.append((s, e)) vcf = GNOMAD_VCF_PATTERN.format(chrom=chrom) variants = [] for s, e in merged: @@ -158,40 +158,40 @@ if faf95 <= PM2_SUPPORTING_THRESHOLD: return 'PM2_supporting', COLORS['PM2_supporting'], bs1_thresh return 'no-code', COLORS['no-code'], bs1_thresh 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() out_dir = os.path.join(args.output_dir, 'cmpVCEPAFfrequencies') os.makedirs(out_dir, exist_ok=True) print(' [B.3 gnomAD v4.1 Allele Frequencies]') - print(f' thresholds: BA1≥{BA1_THRESHOLD}, BS1≥{DEFAULT_BS1} (or 0.0002 MYBPC3), PM2_sup≤{PM2_SUPPORTING_THRESHOLD}') + print(f' thresholds: BA1>={BA1_THRESHOLD}, BS1>={DEFAULT_BS1} (or 0.0002 MYBPC3), PM2_sup<={PM2_SUPPORTING_THRESHOLD}') bed_lines = [] counts = {'BA1': 0, 'BS1': 0, 'PM2_supporting': 0, 'no-code': 0} for gene in OUR_GENES: mane = parse_mane_record(gene) print(f' fetching {gene} ({mane["chrom"]} {mane["strand"]})...') variants = fetch_gene_variants(gene, mane) - print(f' {len(variants)} PASS variants in CDS ±{SPLICE_PADDING} nt') + print(f' {len(variants)} PASS variants in CDS +/-{SPLICE_PADDING} nt') for v in variants: code, color, bs1_thresh = apply_code(v['faf95'], gene) counts[code] += 1 start_bed = v['pos'] - 1 end_bed = v['pos'] - 1 + len(v['ref']) mouseover = ( f'gnomAD v4.1 - {code}
' f'{gene} {v["chrom"]}:{v["pos"]} {v["ref"]}>{v["alt"]}
' f'FAF95 (popmax): {v["faf95"]:.2e}
' f'AF (overall): {v["af"]:.2e}
' f'Popmax group: {v["grpmax"]}
' f'CSpec thresholds for {gene}: BA1≥{BA1_THRESHOLD}, BS1≥{bs1_thresh}, PM2_sup≤{PM2_SUPPORTING_THRESHOLD}' ) name = f'{gene}_{v["ref"]}>{v["alt"]}_{code[:6]}' @@ -210,31 +210,31 @@ mouseover, ])) print(f' total counts: {counts}') bed_lines.sort(key=lambda l: (l.split('\t')[0], int(l.split('\t')[1]))) as_path = os.path.join(out_dir, 'cmpVCEPAFfrequencies.as') with open(as_path, 'w') as f: f.write(AUTOSQL) hg38_bed = os.path.join(out_dir, 'cmpVCEPAFfrequenciesHg38.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, 'cmpVCEPAFfrequenciesHg38.bb') cmd = ['bedToBigBed', '-tab', '-type=bed9+9', '-as=' + as_path, hg38_bed, CHROM_SIZES['hg38'], hg38_bb] print(f' $ {" ".join(cmd)}') subprocess.run(cmd, check=True) print(f' hg38 bigBed: {hg38_bb}') if 'hg19' in args.db: hg19_bed = os.path.join(out_dir, 'cmpVCEPAFfrequenciesHg19.bed') unmapped = hg19_bed + '.unmapped' cmd = ['liftOver', '-bedPlus=9', '-tab', hg38_bed, LIFTOVER_HG38_TO_HG19, hg19_bed, unmapped] print(f' $ {" ".join(cmd)}') subprocess.run(cmd, check=True)