615721361f4baf75c0715bb931c5fc1015101622 lrnassar Tue Aug 11 18:21:56 2026 -0700 Address code-review findings on the Cardiomyopathy VCEP scripts. refs #37446 - Gate PM1 to missense variants per the CSpec ("applicable to missense variants"); a positional-only test wrongly gave synonymous/truncating/splice variants PM1 and let it collide with BA1/BP7. PM1 firing 1,293 -> 700. - Transcript-gate the Walsh-2019 ClinVar coordinate lookup so a classic-vs-MANE c.notation collision no longer mis-places TNNT2 R92Q (was drawn ~331 nt off with a different variant's VariationID); the gate applies only to the WALSH_TX genes. - Show the amino-acid change in the REVEL mouseover (computed from the MANE CDS) so the per-alt genomic-forward-strand score is not misread on minus-strand genes. - Resolve every build input relative to --output-dir (sibling track outputs and cmp_downloads sources) for otto portability; canonical build byte-identical. - Also key the diseaseTag off the counted PM1 code, and makedoc corrections (worked example REVEL/gnomAD values, PM1 count, universe and EvRepo notes). diff --git src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPRevel.py src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPRevel.py index d88aa4f1e0f..f08a8b11a79 100644 --- src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPRevel.py +++ src/hg/makeDb/scripts/cardiomyopathyVCEP/cmpVCEPRevel.py @@ -48,30 +48,88 @@ uint thickStart; "Same as chromStart" uint thickEnd; "Same as chromEnd" uint itemRgb; "PP3 light-purple or BP4 light-orange" string gene; "Gene symbol" char[1] altAllele; "Alternate nucleotide" double revelScore; "REVEL score" string acmgCode; "PP3_Supporting or BP4_Supporting" lstring _mouseOver; "Tooltip HTML" ) """ # Re-use B.1's MANE parsing sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from cmpVCEPClinDomains import parse_mane_record, cds_exons +# Standard genetic code + amino-acid 3-letter names, for the mouseover aa change. +# REVEL is per-alt-nucleotide on the genomic forward strand; without the resulting aa change +# a curator on a minus-strand gene reads the wrong allele's score (e.g. reads R870L for R870H). +CODON1 = { + 'TTT':'F','TTC':'F','TTA':'L','TTG':'L','CTT':'L','CTC':'L','CTA':'L','CTG':'L', + 'ATT':'I','ATC':'I','ATA':'I','ATG':'M','GTT':'V','GTC':'V','GTA':'V','GTG':'V', + 'TCT':'S','TCC':'S','TCA':'S','TCG':'S','CCT':'P','CCC':'P','CCA':'P','CCG':'P', + 'ACT':'T','ACC':'T','ACA':'T','ACG':'T','GCT':'A','GCC':'A','GCA':'A','GCG':'A', + 'TAT':'Y','TAC':'Y','TAA':'*','TAG':'*','CAT':'H','CAC':'H','CAA':'Q','CAG':'Q', + 'AAT':'N','AAC':'N','AAA':'K','AAG':'K','GAT':'D','GAC':'D','GAA':'E','GAG':'E', + 'TGT':'C','TGC':'C','TGA':'*','TGG':'W','CGT':'R','CGC':'R','CGA':'R','CGG':'R', + 'AGT':'S','AGC':'S','AGA':'R','AGG':'R','GGT':'G','GGC':'G','GGA':'G','GGG':'G', +} +AA3 = {'A':'Ala','R':'Arg','N':'Asn','D':'Asp','C':'Cys','Q':'Gln','E':'Glu','G':'Gly', + 'H':'His','I':'Ile','L':'Leu','K':'Lys','M':'Met','F':'Phe','P':'Pro','S':'Ser', + 'T':'Thr','W':'Trp','Y':'Tyr','V':'Val','*':'Ter'} +COMP = str.maketrans('ACGTacgt', 'TGCAtgca') +TWOBIT = '/gbdb/hg38/hg38.2bit' + + +def _fetch_seq(chrom, start, end): + out = subprocess.check_output(['twoBitToFa', f'{TWOBIT}:{chrom}:{start}-{end}', 'stdout'], text=True) + return ''.join(out.splitlines()[1:]).upper() + + +def build_codon_index(mane): + """Return (cds coding sequence, {genomic 0-based pos -> CDS index}) for this MANE transcript. + Positions are in transcription order; minus-strand exons are reverse-complemented.""" + chrom, strand = mane['chrom'], mane['strand'] + seq_parts, order = [], [] + exons = cds_exons(mane) # ascending genomic (start, end), 0-based half-open + for s, e in (exons if strand == '+' else reversed(exons)): + seg = _fetch_seq(chrom, s, e) + if strand == '+': + seq_parts.append(seg); order.extend(range(s, e)) + else: + seq_parts.append(seg.translate(COMP)[::-1]); order.extend(range(e - 1, s - 1, -1)) + cds = ''.join(seq_parts) + return cds, {p: i for i, p in enumerate(order)} + + +def revel_hgvsp(cds, pos2idx, strand, gpos0, alt_fwd): + """(hgvsp, short) for a forward-strand SNV at genomic 0-based gpos0, or None if not missense.""" + idx = pos2idx.get(gpos0) + if idx is None: + return None + cnum, cp = idx // 3, idx % 3 + ref_codon = cds[cnum * 3: cnum * 3 + 3] + if len(ref_codon) < 3: + return None + coding_alt = alt_fwd if strand == '+' else alt_fwd.translate(COMP) + alt_codon = ref_codon[:cp] + coding_alt + ref_codon[cp + 1:] + ra, aa = CODON1.get(ref_codon), CODON1.get(alt_codon) + if ra is None or aa is None or ra == aa: + return None + n = cnum + 1 + return (f'p.{AA3[ra]}{n}{AA3[aa]}', f'{ra}{n}{aa}') + def fetch_revel_bedgraph(chrom, start, end, alt_nt): """Return list of (genomic_start, genomic_end, score) for REVEL alt=alt_nt in this region. bigWigToBedGraph collapses runs of identical scores at adjacent positions into a single multi-bp BED row. REVEL is per-position-per-alt: each position has its own REF allele, so the bedGraph compaction is wrong for our purposes. Split any multi-bp run into N consecutive 1-bp records before returning. (FULL audit P0 #1 fix, 2026-04-28.) """ cmd = ['bigWigToBedGraph', f'-chrom={chrom}', f'-start={start}', f'-end={end}', REVEL_BW[alt_nt], 'stdout'] out = subprocess.check_output(cmd, text=True) rows = [] @@ -96,54 +154,59 @@ os.makedirs(out_dir, exist_ok=True) print(' [B.4 REVEL PP3/BP4]') print(f' thresholds: PP3 if REVEL >= {PP3_THRESHOLD}; BP4 if REVEL <= {BP4_THRESHOLD}') bed_lines = [] n_pp3 = 0 n_bp4 = 0 n_dropped = 0 for gene in OUR_GENES: mane = parse_mane_record(gene) chrom = mane['chrom'] strand = mane['strand'] exons = cds_exons(mane) + cds, pos2idx = build_codon_index(mane) for ex_start, ex_end in exons: for alt_nt in 'acgt': rows = fetch_revel_bedgraph(chrom, ex_start, ex_end, alt_nt) for s, e, score in rows: if score == 0: continue # 0 = not missense / no REVEL score if score >= PP3_THRESHOLD: code = 'PP3_Supporting' color = PP3_COLOR n_pp3 += 1 elif score <= BP4_THRESHOLD: code = 'BP4_Supporting' color = BP4_COLOR n_bp4 += 1 else: n_dropped += 1 continue # indeterminate band - drop per InSiGHT precedent name = f'{gene}_{alt_nt.upper()}_{score:.3f}_{code[:3]}' + hp = revel_hgvsp(cds, pos2idx, strand, s, alt_nt.upper()) + aa_html = f'{hp[0]} ({hp[1]})
' if hp else '' mouseover = ( f'REVEL - {code}
' - f'{gene} {chrom}:{s+1} alt={alt_nt.upper()}
' - f'REVEL score: {score:.3f}
' - f'CSpec threshold: PP3 ≥ {PP3_THRESHOLD}; BP4 ≤ {BP4_THRESHOLD}' + f'{gene} {chrom}:{s+1}
' + + aa_html + + f'REVEL score: {score:.3f} ' + f'(genomic forward-strand alt {alt_nt.upper()})
' + + f'CSpec threshold: PP3 ≥ {PP3_THRESHOLD}; BP4 ≤ {BP4_THRESHOLD}' ) bed_lines.append('\t'.join([ chrom, str(s), str(e), name, '0', strand, str(s), str(e), color, gene, alt_nt.upper(), f'{score:.3f}', code, mouseover, ])) print(f' {gene}: scanned {len(exons)} CDS exons') print(f' total: PP3={n_pp3}, BP4={n_bp4}, dropped indeterminate={n_dropped}')