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/cmpVCEPCardioBoost.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPCardioBoost.py new file mode 100644 index 00000000000..545bec84d68 --- /dev/null +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPCardioBoost.py @@ -0,0 +1,136 @@ +#!/usr/bin/env python3 +""" +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, +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 +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. + +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 + +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. +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" + uint itemRgb; "Color by CardioBoost class" + string gene; "Gene symbol" + char[1] refAllele; "Reference nucleotide" + char[1] altAllele; "Alternate nucleotide" + string hgvsc; "CardioBoost HGVS c. (Ensembl transcript)" + string hgvsp; "CardioBoost HGVS p. (Ensembl protein)" + 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(): + rows = [] + 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') + os.makedirs(out_dir, exist_ok=True) + 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 = 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') + + 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) + + # 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') + + +if __name__ == '__main__': + main()