5fe93cbef0d7692552e37e62bbf394cc7d222a28
max
  Thu Jul 16 00:32:55 2026 -0700
ClinVar Mapped track: map ClinVar coding variants through paralog alignments

#Preview2 week - bugs introduced now will need a build patch to fix
Adds the clinvarMapped composite (hg38) with two subtracks:
- clinvarMappedParalog: every protein-changing ClinVar variant projected onto
the equivalent (aligned) residue of each of its gene's paralogs
- clinvarMappedParalogAln: the pairwise protein alignments used to do the
mapping, as bigPsl, so the evidence for each projection can be inspected

Pipeline: paralog pairs from Ensembl BioMart, one representative transcript per
gene from MANE Select, pairwise global protein alignment (BLOSUM62) for pairs at
>=20% identity, variants projected residue-to-residue and mapped back to the
paralog's genomic codon. Colors and filters mirror the ClinVar track; the map
track is filterable by source gene, classification, review stars, residue
conservation and percent identity.

refs #37883

diff --git src/hg/makeDb/scripts/clinvarMapped/clinvarMappedParalogProject.py src/hg/makeDb/scripts/clinvarMapped/clinvarMappedParalogProject.py
new file mode 100644
index 00000000000..69bcb1a008c
--- /dev/null
+++ src/hg/makeDb/scripts/clinvarMapped/clinvarMappedParalogProject.py
@@ -0,0 +1,207 @@
+#!/usr/bin/env python3
+"""clinvarMappedParalogProject.py - project ClinVar coding variants onto their
+paralogs' equivalent residues.
+
+For each ClinVar variant assigned to a MANE codon (clinvarCodons.tsv) and each
+annotated paralog of the source gene (alnBlocks.tsv), we look up the aligned
+residue in the paralog, map that residue back to the paralog's genomic codon,
+and emit one bed12+ feature on the paralog. Everything needed to reproduce the
+mapping - source gene, source variant, ClinVar link, residue conservation,
+identities - travels on the feature.
+
+Inputs:
+  clinvarCodons.tsv  ensg enst sym aaPos refRes chrom gPos vcv variantId
+                     clinSignCode clinSign stars molConseq protChange
+  alnBlocks.tsv      geneA geneB lenA lenB nIdent matchCols blockCount
+                     aStarts bStarts blockSizes   (protein 0-based)
+  pairs.ge20.tsv     geneA geneB symA symB pidAB pidBA maxPid
+  maneMeta.tsv       enst ensg sym refTx refProt chrom strand geneType
+  maneSelect.gp      genePred
+  geneProt.faa       protein per gene (">ENSG...")
+Output (-o): unsorted bed12+14 (see clinvarMappedParalog.as)
+
+Usage: clinvarMappedParalogProject.py <dataDir> -o out.bed [-j N]
+"""
+import argparse, sys, os, array, multiprocessing as mp
+from Bio.Align import substitution_matrices
+
+BLOSUM = substitution_matrices.load("BLOSUM62")
+
+# colors by clinical significance, matching the UCSC ClinVar track (clinVarToBed)
+RGB = {
+    "PG": "210,0,0",     "LP": "210,0,0",     # pathogenic / likely pathogenic
+    "BN": "0,210,0",     "LB": "0,210,0",     # benign / likely benign
+    "VUS": "0,0,128",                          # uncertain
+    "CF": "137,121,212",                       # conflicting
+    "OT": "128,128,128", "RF": "128,128,128",  # other / risk factor
+}
+
+# module globals shared with workers via fork
+prot = {}          # ensg -> protein sequence
+cdsG = {}          # ensg -> array('i') genomic coords in transcription order
+strandOf = {}      # ensg -> '+'/'-'
+symOf = {}         # ensg -> symbol
+partners = {}      # ensg -> list of (partnerEnsg, srcStarts, srcSizes, tgtStarts, alnPid)
+bioPid = {}        # (a,b) sorted -> max BioMart percent id
+
+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):
+    """ensg -> array of genomic coords for every CDS base in transcription order."""
+    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):
+    """Keep bigBed string fields under the 255-char limit (long polyQ ins/dels)."""
+    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):
+    """3 genomic coords (any order) -> (start, end, blockSizes, blockStarts)."""
+    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
+
+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]
+        plist = partners.get(srcEnsg)
+        if not plist: continue
+        aa0 = aaPos - 1
+        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)
+        for (pEnsg, srcStarts, srcSizes, tgtStarts, alnPid) in plist:
+            # find block containing aa0 in the source coordinate
+            tgt0 = None
+            for bi in range(len(srcStarts)):
+                s = srcStarts[bi]
+                if s <= aa0 < s + srcSizes[bi]:
+                    tgt0 = tgtStarts[bi] + (aa0 - s); break
+            if tgt0 is None: continue                  # aligns to a gap
+            pprot = prot.get(pEnsg, "")
+            if tgt0 >= len(pprot): continue
+            pRes = pprot[tgt0]
+            g = cdsG.get(pEnsg)
+            base = 3 * tgt0
+            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)
+            pChrom = chromOf[pEnsg]
+            pChange = protChange or ("p.%s%d" % (refRes, aaPos))
+            name = shorten("%s %s" % (srcSym, pChange), 120)
+            rc = resClass(refRes, pRes)
+            bpid = bioPid.get((srcEnsg, pEnsg) if srcEnsg < pEnsg else (pEnsg, srcEnsg), 0.0)
+            row = [pChrom, start, end, name, score, strandOf[pEnsg],
+                   start, end, rgb, len(sizes),
+                   ",".join(map(str, sizes)) + ",", ",".join(map(str, starts)) + ",",
+                   srcSym, shorten(pChange, 200), aaPos, refRes, pRes, rc,
+                   clinSign, clinCode, starsN, molConseq,
+                   "%.1f" % bpid, "%.1f" % alnPid, vcv, sourceLocus]
+            outRows.append("\t".join(map(str, row)))
+    return outRows
+
+chromOf = {}
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("dataDir")
+    ap.add_argument("-o", required=True); ap.add_argument("-j", type=int, default=32)
+    args = ap.parse_args()
+    D = args.dataDir
+
+    global prot, cdsG, strandOf, symOf, partners, bioPid, chromOf
+    prot = loadFasta(os.path.join(D, "geneProt.faa"))
+    cdsG = buildCdsG(os.path.join(D, "maneSelect.gp"))
+    for line in open(os.path.join(D, "maneMeta.tsv")):
+        f = line.rstrip("\n").split("\t")
+        # map by ENSG using this gene's transcript; enst is f[0]
+        strandOf[f[1]] = f[6]; symOf[f[1]] = f[2]; chromOf[f[1]] = f[5]
+    # cdsG is keyed by ENST; re-key to ENSG (one transcript per gene)
+    enst2ensg = {}
+    for line in open(os.path.join(D, "maneMeta.tsv")):
+        f = line.rstrip("\n").split("\t"); enst2ensg[f[0]] = f[1]
+    cdsG = {enst2ensg[k]: v for k, v in cdsG.items() if k in enst2ensg}
+    prot = prot   # already ENSG-keyed
+
+    for line in open(os.path.join(D, "pairs.ge20.tsv")):
+        f = line.rstrip("\n").split("\t")
+        a, b = f[0], f[1]
+        bioPid[(a, b) if a < b else (b, a)] = float(f[6])
+
+    # build the partner index from alnBlocks (both orientations)
+    n = 0
+    for line in open(os.path.join(D, "alnBlocks.tsv")):
+        f = line.rstrip("\n").split("\t")
+        a, b = f[0], f[1]
+        nIdent, matchCols = int(f[4]), int(f[5])
+        pid = 100.0 * nIdent / matchCols if matchCols else 0.0
+        aStarts = [int(x) for x in f[7].split(",")]
+        bStarts = [int(x) for x in f[8].split(",")]
+        sizes   = [int(x) for x in f[9].split(",")]
+        partners.setdefault(a, []).append((b, aStarts, sizes, bStarts, pid))
+        partners.setdefault(b, []).append((a, bStarts, sizes, aStarts, pid))
+        n += 1
+    sys.stderr.write("indexed %d pairs, %d genes with partners, %d proteins\n"
+                     % (n, len(partners), len(prot)))
+
+    lines = open(os.path.join(D, "clinvarCodons.tsv")).readlines()
+    sys.stderr.write("projecting %d variant-codon rows on %d workers\n"
+                     % (len(lines), args.j))
+    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()