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/cmpVCEPAnnotate.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAnnotate.py new file mode 100644 index 00000000000..ff57aa546b2 --- /dev/null +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPAnnotate.py @@ -0,0 +1,189 @@ +#!/usr/bin/env python3 +""" +Phase 1 — hgVai annotation layer for the Cardiomyopathy VCEP hub (RM37446). + +Annotates the gnomAD-observed variant universe (the SAME universe used by the +AF + Provisional tracks) with protein consequence via the internal UCSC tool +hgVai (vai.pl), on the single MANE Select transcript per gene (ncbiRefSeqSelect). +Reuses cmpVCEPAFfrequencies.fetch_gene_variants so the universe is identical. + +Output (consumed by Phase 2 cmpVCEPProvisionalClass): + cmpVCEPAnnotate/cmpVCEPAnnotations.hg38.tsv one row per (chrom,pos,ref,alt) + cmpVCEPAnnotate/phase1_unmapped.tsv universe variants with no VEP line + cmpVCEPAnnotate/<gene>.vcf / <gene>.vep per-gene intermediates + +Annotation is computed on hg38 only; Phase 2 codes carry to hg19 via the existing +liftOver of the universe (annotation is transcript-based, assembly-independent). + +Usage: + python3 cmpVCEPAnnotate.py --output-dir /hive/users/lrnassar/claude/RM37446 +""" +import argparse +import os +import re +import subprocess +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import cmpVCEPAFfrequencies as af + +VAI = os.path.expanduser('~/bin/scripts/vai.pl') +HGVAI = '/usr/local/apache/cgi-bin/hgVai' +GENE_TRACK = 'ncbiRefSeqSelect' + +# We set the VCF ID column to "chrom:pos:ref:alt"; vai.pl echoes it into the +# "Uploaded Variation" column verbatim, giving an exact join even for indels +# (whose default chrom_pos_ref/alt name VEP reformats). + + +def write_gene_vcf(variants, path): + """Write a sorted, deduped minimal VCF for one gene's universe variants.""" + seen = set() + rows = [] + for v in variants: + key = (v['chrom'], v['pos'], v['ref'], v['alt']) + if key in seen: + continue + seen.add(key) + rows.append(key) + rows.sort(key=lambda k: k[1]) + with open(path, 'w') as fh: + fh.write('##fileformat=VCFv4.2\n') + fh.write('#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n') + for chrom, pos, ref, alt in rows: + vid = f'{chrom}:{pos}:{ref}:{alt}' + fh.write(f'{chrom}\t{pos}\t{vid}\t{ref}\t{alt}\t.\t.\t.\n') + return rows + + +def run_vai(vcf_path, vep_path, chrom, tx_start, tx_end): + """Run hgVai over the gene tx region; write VEP output to vep_path.""" + region = f'{chrom}:{tx_start + 1}-{tx_end}' + cmd = [VAI, '--variantLimit=200000000', f'--hgVai={HGVAI}', + f'--position={region}', f'--geneTrack={GENE_TRACK}', + '--hgvsG=off', '--hgvsCN=off', '--hgvsP=on', 'hg38', vcf_path] + with open(vep_path, 'w') as out: + subprocess.run(cmd, stdout=out, stderr=subprocess.DEVNULL, check=True) + + +def parse_extra(extra): + """Pull HGVSP and EXON out of the VEP Extra column (key=value;key=value).""" + hgvsp, exon_n, exon_total = '', '', '' + for kv in extra.split(';'): + if kv.startswith('HGVSP='): + hgvsp = kv[len('HGVSP='):] + elif kv.startswith('EXON='): + m = re.match(r'(\d+)/(\d+)', kv[len('EXON='):]) + if m: + exon_n, exon_total = m.group(1), m.group(2) + return hgvsp, exon_n, exon_total + + +def parse_vep(vep_path, gene): + """Parse one gene's VEP output → dict key (chrom,pos,ref,alt) → annotation. + Aggregates SO terms across lines; keeps the coding fields from the line(s) + that carry them.""" + ann = {} + for line in open(vep_path): + if line.startswith('#') or line.startswith('Uploaded'): + continue + f = line.rstrip('\n').split('\t') + if len(f) < 14: + continue + parts = f[0].split(':') + if len(parts) != 4: + continue + chrom, pos, ref, alt = parts[0], int(parts[1]), parts[2], parts[3] + key = (chrom, pos, ref, alt) + so = f[6] + cdna_pos = f[7] # "Position in cDNA" (transcript orientation; for NMD-escape) + protein_pos = f[9] + aa_change = f[10] # "R/H" (missense), "R" (syn), or "-" + codon = f[11] + hgvsp, exon_n, exon_total = parse_extra(f[13]) + + rec = ann.setdefault(key, { + 'gene': gene, 'so': set(), 'proteinPos': '', 'aaRef': '', 'aaAlt': '', + 'codon': '', 'exonNum': '', 'exonTotal': '', 'hgvsp': '', 'cdnaPos': '', + }) + rec['so'].add(so) + # Capture coding fields from whichever line carries them (non-'-'). + if cdna_pos not in ('', '-') and not rec['cdnaPos']: + rec['cdnaPos'] = cdna_pos.split('-')[0] # range start for indels + if protein_pos not in ('', '-') and not rec['proteinPos']: + rec['proteinPos'] = protein_pos + if codon not in ('', '-') and not rec['codon']: + rec['codon'] = codon + if exon_total and not rec['exonTotal']: + rec['exonNum'], rec['exonTotal'] = exon_n, exon_total + if hgvsp and not rec['hgvsp']: + rec['hgvsp'] = hgvsp + if aa_change not in ('', '-') and not rec['aaRef']: + if '/' in aa_change: + rec['aaRef'], rec['aaAlt'] = aa_change.split('/', 1) + else: + rec['aaRef'] = rec['aaAlt'] = aa_change # synonymous: ref==alt + return ann + + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument('--output-dir', required=True) + args = ap.parse_args() + + work = os.path.join(args.output_dir, 'cmpVCEPAnnotate') + os.makedirs(work, exist_ok=True) + + print(' [Phase 1 hgVai annotation]') + all_ann = {} + unmapped = [] + for gene in af.OUR_GENES: + mane = af.parse_mane_record(gene) + variants = af.fetch_gene_variants(gene, mane) + vcf_path = os.path.join(work, f'{gene}.vcf') + vep_path = os.path.join(work, f'{gene}.vep') + universe_keys = write_gene_vcf(variants, vcf_path) + run_vai(vcf_path, vep_path, mane['chrom'], mane['chromStart'], mane['chromEnd']) + ann = parse_vep(vep_path, gene) + # Detect universe variants that got no VEP annotation + ann_keys = set(ann.keys()) + for k in universe_keys: + if k not in ann_keys: + unmapped.append((gene,) + k) + all_ann.update(ann) + print(f' {gene} ({mane["chrom"]} {mane["strand"]}): ' + f'{len(universe_keys)} universe, {len(ann)} annotated, ' + f'{sum(1 for k in universe_keys if k not in ann_keys)} unmapped') + + tsv = os.path.join(work, 'cmpVCEPAnnotations.hg38.tsv') + with open(tsv, 'w') as fh: + fh.write('#chrom\tpos\tref\talt\tgene\tsoTerms\tproteinPos\taaRef\taaAlt\t' + 'codonChange\texonNum\texonTotal\thgvsp\tcdnaPos\n') + for (chrom, pos, ref, alt), r in sorted(all_ann.items(), key=lambda x: (x[0][0], x[0][1])): + fh.write('\t'.join([ + chrom, str(pos), ref, alt, r['gene'], + ','.join(sorted(r['so'])), r['proteinPos'], r['aaRef'], r['aaAlt'], + r['codon'], r['exonNum'], r['exonTotal'], r['hgvsp'], r['cdnaPos'], + ]) + '\n') + print(f' wrote {len(all_ann)} annotations → {tsv}') + + unmapped_path = os.path.join(work, 'phase1_unmapped.tsv') + with open(unmapped_path, 'w') as fh: + fh.write('#gene\tchrom\tpos\tref\talt\n') + for row in unmapped: + fh.write('\t'.join(str(x) for x in row) + '\n') + print(f' {len(unmapped)} unmapped universe variants logged → {unmapped_path}') + + # SO-term summary for sanity + from collections import Counter + so_counter = Counter() + for r in all_ann.values(): + for s in r['so']: + so_counter[s] += 1 + print(' SO-term distribution:') + for s, c in so_counter.most_common(): + print(f' {c:6d} {s}') + + +if __name__ == '__main__': + main()