c9446936b026a803d97e7f5128404a327829fe2e max Sun Jul 19 00:01:22 2026 -0700 ClinVar Mapped: add Pfam-domain mapping method Maps ClinVar coding variants through Pfam protein domains as a second method under the clinvarMapped container. Each Pfam family's profile HMM gives a shared coordinate, so a variant inside a domain is projected to the residue in the same HMM match-state column in every other MANE Select gene carrying that domain, then mapped back to that gene's genomic codon. This links genes that share a domain without being paralogs (e.g. the SCN sodium channels and the CACNA1 calcium channels via Ion_trans). Each variant is shown at its 25 most conserved equivalents. Uses Pfam-A 38.2. Adds two alpha subtracks (Pfam Domain Variants, Pfam Domains), the hmmsearch/hmmalign build scripts, autoSql, makeDoc and description page, refs #37883 diff --git src/hg/makeDb/scripts/clinvarMapped/clinvarMappedPfamProject.py src/hg/makeDb/scripts/clinvarMapped/clinvarMappedPfamProject.py new file mode 100644 index 00000000000..79688449993 --- /dev/null +++ src/hg/makeDb/scripts/clinvarMapped/clinvarMappedPfamProject.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python3 +"""clinvarMappedPfamProject.py - project ClinVar coding variants onto every other +gene that shares the same Pfam domain, at the equivalent HMM match-state column. + +For each ClinVar variant already assigned to a MANE codon (clinvarCodons.tsv) we +look up which Pfam domain column the source residue occupies (alnCols.tsv). Every +other instance of that domain that has a residue in the same column is an +equivalent position; we map that residue back to its gene's genomic codon and emit +one bed12+ feature there. The Pfam HMM is the shared coordinate hub, so this links +all genes carrying a domain, not only annotated paralog pairs. + +Inputs (in dataDir): + clinvarCodons.tsv ensg enst sym aaPos refRes chrom gPos vcv variantId + clinSignCode clinSign stars molConseq protChange + pfam/alnCols.tsv ensg pfamAcc pfamName aaPos matchCol + maneMeta.tsv enst ensg sym refTx refProt chrom strand geneType + maneSelect.gp genePred + geneProt.faa protein per gene (">ENSG...") +Output (-o): unsorted bed12+ (see clinvarMappedPfam.as) + +Usage: clinvarMappedPfamProject.py <dataDir> -o out.bed [--maxFanout N] [-j N] +""" +import argparse, sys, os, array, multiprocessing as mp +from Bio.Align import substitution_matrices + +BLOSUM = substitution_matrices.load("BLOSUM62") + +RGB = { + "PG": "210,0,0", "LP": "210,0,0", + "BN": "0,210,0", "LB": "0,210,0", + "VUS": "0,0,128", + "CF": "137,121,212", + "OT": "128,128,128", "RF": "128,128,128", +} + +prot = {} # ensg -> protein sequence +cdsG = {} # ensg -> array('i') genomic coords in transcription order +strandOf = {} +symOf = {} +chromOf = {} +byRes = {} # (ensg, aaPos) -> list of (pfamAcc, pfamName, col) +col2inst = {} # (pfamAcc, col) -> list of (ensg, aaPos) +maxFanout = 0 + +def loadFasta(path): + d, name, buf = {}, None, [] + for line in open(path): + if line[0] == ">": + if name: d[name] = "".join(buf) + name = line[1:].split()[0]; buf = [] + else: buf.append(line.strip()) + if name: d[name] = "".join(buf) + return d + +def buildCdsG(gpPath): + out = {} + for line in open(gpPath): + f = line.rstrip("\n").split("\t") + enst, chrom, strand = f[0], f[1], f[2] + cdsStart, cdsEnd = int(f[5]), int(f[6]) + if cdsStart >= cdsEnd: continue + exStarts = [int(x) for x in f[8].rstrip(",").split(",")] + exEnds = [int(x) for x in f[9].rstrip(",").split(",")] + segs = [] + for s, e in zip(exStarts, exEnds): + cs, ce = max(s, cdsStart), min(e, cdsEnd) + if cs < ce: segs.append((cs, ce)) + a = array.array('i') + if strand == "+": + for (s, e) in segs: a.extend(range(s, e)) + else: + for (s, e) in reversed(segs): a.extend(range(e - 1, s - 1, -1)) + out[enst] = a + return out + +def shorten(s, n): + return s if len(s) <= n else s[:n - 3] + "..." + +def resClass(a, b): + if a == b: return "identical" + try: + return "similar" if BLOSUM[a, b] > 0 else "different" + except (KeyError, IndexError): + return "different" + +def codonBlocks(coords): + cs = sorted(coords) + start, end = cs[0], cs[-1] + 1 + sizes, starts = [], [] + i = 0 + while i < len(cs): + j = i + while j + 1 < len(cs) and cs[j + 1] == cs[j] + 1: + j += 1 + starts.append(cs[i] - start); sizes.append(cs[j] - cs[i] + 1) + i = j + 1 + return start, end, sizes, starts + +RANK = {"identical": 0, "similar": 1, "different": 2} + +def projectChunk(lines): + outRows = [] + for line in lines: + f = line.rstrip("\n").split("\t") + srcEnsg, srcSym, aaPos, refRes = f[0], f[2], int(f[3]), f[4] + vcv, clinCode, clinSign, stars, molConseq, protChange = \ + f[7], f[9], f[10], f[11], f[12], f[13] + srcChrom, srcGpos = f[5], f[6] + cols = byRes.get((srcEnsg, aaPos)) + if not cols: continue # source residue not in any domain + # gather unique candidate targets across all domains/columns of this residue + cand = {} # (tgtEnsg, tgtAaPos) -> (pfamAcc, pfamName, col) + for (pfamAcc, pfamName, col) in cols: + for (tgtEnsg, tgtAaPos) in col2inst.get((pfamAcc, col), ()): + if tgtEnsg == srcEnsg: continue # skip self + cand.setdefault((tgtEnsg, tgtAaPos), (pfamAcc, pfamName, col)) + if not cand: continue + # rank each candidate by residue conservation vs the source residue + ranked = [] + for (tgtEnsg, tgtAaPos), (pfamAcc, pfamName, col) in cand.items(): + pprot = prot.get(tgtEnsg, "") + if tgtAaPos > len(pprot): continue + tgtRes = pprot[tgtAaPos - 1] + rc = resClass(refRes, tgtRes) + ranked.append((RANK[rc], tgtEnsg, tgtAaPos, tgtRes, rc, pfamAcc, pfamName, col)) + # keep the most conserved first (identical > similar > different), deterministic + ranked.sort(key=lambda x: (x[0], x[1], x[2])) + if maxFanout: + ranked = ranked[:maxFanout] + + rgb = RGB.get(clinCode, "0,0,0") + try: starsN = int(stars) + except ValueError: starsN = 0 + score = min(1000, starsN * 250) + sourceLocus = "%s:%s" % (srcChrom, int(srcGpos) + 1) + pChange = protChange or ("p.%s%d" % (refRes, aaPos)) + name = shorten("%s %s" % (srcSym, pChange), 120) + for (_rank, tgtEnsg, tgtAaPos, tgtRes, rc, pfamAcc, pfamName, col) in ranked: + g = cdsG.get(tgtEnsg) + base = 3 * (tgtAaPos - 1) + if g is None or base + 3 > len(g): continue + coords = [g[base], g[base + 1], g[base + 2]] + start, end, sizes, starts = codonBlocks(coords) + tgtSym = symOf.get(tgtEnsg, tgtEnsg) + row = [chromOf[tgtEnsg], start, end, name, score, strandOf[tgtEnsg], + start, end, rgb, len(sizes), + ",".join(map(str, sizes)) + ",", ",".join(map(str, starts)) + ",", + srcSym, tgtSym, shorten(pChange, 200), aaPos, refRes, tgtRes, rc, + clinSign, clinCode, starsN, molConseq, + pfamName, pfamAcc, col, vcv, sourceLocus] + outRows.append("\t".join(map(str, row))) + return outRows + +def main(): + ap = argparse.ArgumentParser() + ap.add_argument("dataDir") + ap.add_argument("-o", required=True) + ap.add_argument("--maxFanout", type=int, default=25, + help="cap targets per source variant, keeping the most conserved (0 = unlimited)") + ap.add_argument("-j", type=int, default=32) + args = ap.parse_args() + D = args.dataDir + + global prot, cdsG, strandOf, symOf, chromOf, byRes, col2inst, maxFanout + maxFanout = args.maxFanout + prot = loadFasta(os.path.join(D, "geneProt.faa")) + cdsG = buildCdsG(os.path.join(D, "maneSelect.gp")) + enst2ensg = {} + for line in open(os.path.join(D, "maneMeta.tsv")): + f = line.rstrip("\n").split("\t") + enst2ensg[f[0]] = f[1] + strandOf[f[1]] = f[6]; symOf[f[1]] = f[2]; chromOf[f[1]] = f[5] + cdsG = {enst2ensg[k]: v for k, v in cdsG.items() if k in enst2ensg} + + nrow = 0 + for line in open(os.path.join(D, "pfam", "alnCols.tsv")): + f = line.rstrip("\n").split("\t") + ensg, pfamAcc, pfamName, aaPos, col = f[0], f[1], f[2], int(f[3]), int(f[4]) + byRes.setdefault((ensg, aaPos), []).append((pfamAcc, pfamName, col)) + col2inst.setdefault((pfamAcc, col), []).append((ensg, aaPos)) + nrow += 1 + sys.stderr.write("indexed %d residue-column rows, %d residues, %d columns\n" + % (nrow, len(byRes), len(col2inst))) + + lines = open(os.path.join(D, "clinvarCodons.tsv")).readlines() + sys.stderr.write("projecting %d variant-codon rows on %d workers (maxFanout=%d)\n" + % (len(lines), args.j, maxFanout)) + chunk = 20000 + chunks = [lines[i:i + chunk] for i in range(0, len(lines), chunk)] + nout = 0 + with open(args.o, "w") as out, mp.Pool(args.j) as pool: + for rows in pool.imap_unordered(projectChunk, chunks): + for r in rows: out.write(r + "\n") + nout += len(rows) + sys.stderr.write("wrote %d projected features\n" % nout) + +if __name__ == "__main__": + main()