aa5669fe641fb39d7711eb81ec05092d14f267fb
lrnassar
Tue Jun 30 15:20:56 2026 -0700
Adding Cardiomyopathy VCEP track hub build scripts and makedoc. refs #37446
Adds the 12 per-track build scripts under
src/hg/makeDb/scripts/cardiomyopathyVCEP/ (gnomAD v4.1 allele frequencies,
REVEL, CardioBoost, the hgVai consequence/HGVSp annotation layer, ClinGen
EvRepo, ClinVar submitter 506161, Walsh 2019 curations, PM1 clinical-domain
hotspots, MYBPC3 PVS1 caveats, Walsh 2017 PS4 odds-ratio track, Atlas PS4
per-variant OR, and the NON-FINAL provisional classifier) plus the build
documentation at src/hg/makeDb/doc/Cardiomyopathy.txt.
All ACMG thresholds are taken directly from the ClinGen Cardiomyopathy CSpecs
(8 genes, affiliation 50002); no thresholds are invented in the build.
diff --git src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinDomains.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinDomains.py
new file mode 100644
index 00000000000..aab352f9a8a
--- /dev/null
+++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPClinDomains.py
@@ -0,0 +1,339 @@
+#!/usr/bin/env python3
+"""
+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)
+
+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
+
+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)"
+ uint score; "Score (always 0)"
+ char[1] strand; "Strand"
+ uint thickStart; "Same as chromStart"
+ uint thickEnd; "Same as chromEnd"
+ uint itemRgb; "Display color"
+ string gene; "Gene symbol"
+ string transcript; "MANE Select transcript"
+ string aaRange; "Amino acid range (e.g. 167-931)"
+ string exonInfo; "Exon containing this segment"
+ lstring _mouseOver; "Tooltip HTML"
+ )
+"""
+
+
+# ------------------------------------------------------------------
+# Codon-to-genomic conversion (clean rewrite of InSiGHT helper).
+# ------------------------------------------------------------------
+
+def parse_mane_record(gene_symbol):
+ """Extract bigGenePred fields for one gene from /gbdb/hg38/mane/mane.bb."""
+ out = subprocess.check_output(['bigBedToBed', MANE_BB, 'stdout'], text=True)
+ for line in out.splitlines():
+ f = line.split('\t')
+ # bigGenePred col 19 (0-indexed 18) is geneName2 (HGNC symbol)
+ if len(f) >= 20 and f[18] == gene_symbol:
+ return {
+ 'chrom': f[0],
+ '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."""
+ 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.
+
+ Algorithm:
+ 1. Compute mRNA-CDS nucleotide range (1-based inclusive): [nt_lo, nt_hi]
+ nt_lo = (aa_start - 1) * 3 + 1
+ nt_hi = aa_end * 3
+ 2. Walk CDS exons in TRANSCRIPT order (high-to-low genomic for minus strand,
+ low-to-high for plus strand), tracking cumulative mRNA position.
+ 3. For each exon, find overlap with [nt_lo, nt_hi]; project back to genomic
+ coords using the appropriate strand mapping.
+
+ Returns list of (genomic_start, genomic_end, exon_idx_in_transcript_order).
+ """
+ nt_lo = (aa_start - 1) * 3 + 1
+ nt_hi = aa_end * 3
+
+ exons_genomic = cds_exons(mane) # genomic-order
+ if mane['strand'] == '+':
+ exons_tx = list(enumerate(exons_genomic, start=1)) # tx-order = genomic-order
+ else:
+ # tx-order is reverse genomic (highest coord = first in transcript)
+ exons_tx = list(enumerate(reversed(exons_genomic), start=1))
+
+ segments = []
+ mrna_consumed = 0 # nt count already walked through earlier exons in transcript-order
+ for exon_idx, (gs, ge) in exons_tx:
+ exon_len = ge - gs
+ exon_mrna_lo = mrna_consumed + 1 # 1-based inclusive
+ exon_mrna_hi = mrna_consumed + exon_len # 1-based inclusive
+ mrna_consumed = exon_mrna_hi
+
+ if exon_mrna_hi < nt_lo or exon_mrna_lo > nt_hi:
+ continue # this exon entirely outside our range
+
+ # overlap region in mRNA-coords
+ ov_lo = max(nt_lo, exon_mrna_lo)
+ ov_hi = min(nt_hi, exon_mrna_hi)
+ # offset within the exon (0-based, inclusive lo, inclusive hi)
+ off_lo = ov_lo - exon_mrna_lo
+ off_hi = ov_hi - exon_mrna_lo
+
+ if mane['strand'] == '+':
+ # tx-order = genomic-order: low offset -> low genomic
+ seg_start = gs + off_lo
+ seg_end = gs + off_hi + 1
+ else:
+ # tx-order = reverse genomic: high genomic = low transcript-offset
+ seg_end = ge - off_lo
+ seg_start = ge - off_hi - 1
+
+ segments.append((seg_start, seg_end, exon_idx))
+
+ return segments
+
+
+# ------------------------------------------------------------------
+# Unit tests on every minus-strand gene + the plus-strand gene
+# ------------------------------------------------------------------
+
+def unit_test_codon_conversion():
+ """Verify codon-to-genomic conversion against known landmarks for each gene.
+
+ Tests:
+ - First codon (codon 1) maps to first 3 nt of CDS
+ - Last codon maps to the stop codon (last 3 nt of CDS)
+ - Length sanity: aa_to_genomic_segments() total bp = (aa_end - aa_start + 1) * 3
+
+ All 7 minus-strand genes + TPM1 (only plus-strand) tested.
+ """
+ print(' [B.1 unit tests]')
+ test_genes = ['MYH7', 'MYBPC3', 'TNNT2', 'TNNI3', 'TPM1', 'ACTC1', 'MYL2', 'MYL3']
+ failures = 0
+ for gene in test_genes:
+ mane = parse_mane_record(gene)
+ cds_total_nt = sum(ee - es for es, ee in cds_exons(mane))
+ last_aa = cds_total_nt // 3 # includes stop
+
+ # Test 1: codon 1 spans 3 bp
+ segs1 = aa_to_genomic_segments(1, 1, mane)
+ bp1 = sum(e - s for s, e, _ in segs1)
+ ok1 = bp1 == 3
+
+ # Test 2: last codon spans 3 bp
+ segs_last = aa_to_genomic_segments(last_aa, last_aa, mane)
+ bp_last = sum(e - s for s, e, _ in segs_last)
+ ok_last = bp_last == 3
+
+ # Test 3: full CDS
+ segs_full = aa_to_genomic_segments(1, last_aa, mane)
+ bp_full = sum(e - s for s, e, _ in segs_full)
+ ok_full = bp_full == cds_total_nt
+
+ # Test 4 (minus-strand only): first codon should be at the HIGH end of the CDS
+ if mane['strand'] == '-':
+ ok_strand = max(s for s, e, _ in segs1) > min(s for s, e, _ in segs_last)
+ else:
+ ok_strand = min(s for s, e, _ in segs1) < min(s for s, e, _ in segs_last)
+
+ status = '✓' if (ok1 and ok_last and ok_full and ok_strand) else '✗'
+ print(f' {status} {gene} ({mane["strand"]}): cds_nt={cds_total_nt} last_aa={last_aa} '
+ f'codon1={bp1}bp last={bp_last}bp full={bp_full}bp strand_order={"OK" if ok_strand else "WRONG"}')
+ if not (ok1 and ok_last and ok_full and ok_strand):
+ failures += 1
+ if failures:
+ print(f' [B.1 unit tests FAILED for {failures} gene(s)]', file=sys.stderr)
+ sys.exit(1)
+ print(' [B.1 unit tests PASSED for all 8 genes]')
+
+
+# ------------------------------------------------------------------
+# Build BED + bigBed
+# ------------------------------------------------------------------
+
+def emit_bed(output_path):
+ """Build hg38 BED from PM1 region definitions."""
+ bed_lines = []
+ for gene, (transcript, ranges) in PM1_REGIONS.items():
+ mane = parse_mane_record(gene)
+ if mane['refSeqAcc'] != transcript:
+ print(f' WARNING: {gene} MANE refSeqAcc {mane["refSeqAcc"]} != expected {transcript}',
+ file=sys.stderr)
+ for aa_start, aa_end, label in ranges:
+ segments = aa_to_genomic_segments(aa_start, aa_end, mane)
+ for seg_start, seg_end, exon_idx in segments:
+ name = f'{gene}_{label.replace(" ", "_")}_aa{aa_start}-{aa_end}_ex{exon_idx}'
+ aa_range = f'{aa_start}-{aa_end}'
+ exon_info = f'exon-{exon_idx}-of-tx'
+ mouseover = (
+ f'PM1 hotspot - Cardiomyopathy VCEP
'
+ f'Gene: {gene} ({transcript})
'
+ f'Region: codons {aa_range} ({label})
'
+ f'Source: ClinGen Cardiomyopathy CSpec; calibrated by Walsh 2019 (PMID 30696458)'
+ )
+ bed_lines.append('\t'.join([
+ mane['chrom'],
+ str(seg_start),
+ str(seg_end),
+ name,
+ '0',
+ mane['strand'],
+ 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}')
+ 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
+ cmd = ['liftOver', '-bedPlus=9', '-tab', hg38_bed, LIFTOVER_HG38_TO_HG19, hg19_bed, unmapped]
+ print(f' $ {" ".join(cmd)}')
+ subprocess.run(cmd, check=True)
+ 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
+ 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}',
+ file=sys.stderr)
+ else:
+ print(f' cross-assembly parity OK: {n_hg38} features each')
+
+
+if __name__ == '__main__':
+ main()