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/clinvarMappedCodons.py src/hg/makeDb/scripts/clinvarMapped/clinvarMappedCodons.py
new file mode 100644
index 00000000000..fb0685d2a84
--- /dev/null
+++ src/hg/makeDb/scripts/clinvarMapped/clinvarMappedCodons.py
@@ -0,0 +1,224 @@
+#!/usr/bin/env python3
+"""clinvarMappedCodons.py - assign ClinVar coding variants to a MANE Select
+codon position.
+
+For every ClinVar variant with a codon-level molecular consequence we find the
+MANE Select CDS exon(s) it overlaps, compute the 1-based amino-acid position in
+that transcript, and record the reference residue (from the MANE protein). A
+variant that overlaps two MANE genes (overlapping loci) yields one row per gene.
+
+The genomic<->CDS mapping is built from the same genePred that produced the
+protein FASTA, so the reference residue is self-consistent with the protein by
+construction (verified by --selftest).
+
+Inputs:
+  clinvarBb   /gbdb/hg38/bbi/clinvar/clinvarMain.bb
+  maneGp      maneSelect.gp
+  geneFaa     geneProt.faa  (headers ">ENSG...")
+  maneMeta    maneMeta.tsv  (enst ensg sym ...)
+Output (-o):
+  ensg enst sym aaPos refRes chrom gPos vcvId variantId clinSignCode clinSign stars molConseq protChange
+
+Usage: clinvarMappedCodons.py <clinvarBb> <maneGp> <geneFaa> <maneMeta> -o out [--selftest]
+"""
+import argparse, subprocess, sys, os, tempfile, json, re
+
+# protein HGVS like "p.Glu36Gly" -> capture (3-letter ref, position)
+_PPAT = re.compile(r'p\.([A-Z][a-z]{2})(\d+)')
+
+def protChangeForPos(jsonStr, aaPos):
+    """From ClinVar's HGVS table pick the p.change whose position == aaPos (the
+    MANE codon); fall back to the first p.change. Returns '' if none."""
+    try:
+        tbl = json.loads(jsonStr)
+    except Exception:
+        return ""
+    first = ""
+    for entry in tbl:
+        for cell in entry:
+            m = _PPAT.search(cell)
+            if not m:
+                continue
+            if not first:
+                first = cell[m.start():].split()[0].rstrip(",")
+            if int(m.group(2)) == aaPos:
+                return cell[m.start():].split()[0].rstrip(",")
+    return first
+
+# molecular consequences that correspond to a single residue position
+CODING = {
+    "missense variant", "nonsense", "inframe deletion", "inframe insertion",
+    "inframe indel", "initiator codon variant", "stop lost",
+}
+# clinvarMain.bed extra-field 1-based columns (12 BED + extras)
+C_MOLCONSEQ, C_CLINSIGN, C_STARS, C_VARID, C_VCV, C_CLINCODE = 18, 14, 45, 46, 48, 41
+
+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 parseGenePred(path, enst2gene):
+    """Return dict txid -> (chrom, strand, ensg, sym, [(gStart,gEnd,cumBefore)])
+    where the exon list is CDS-clipped and cumBefore is the count of CDS bases
+    preceding the exon in transcription order."""
+    tx = {}
+    for line in open(path):
+        f = line.rstrip("\n").split("\t")
+        name, 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(",")]
+        # CDS-clipped exon segments in genomic order
+        segs = []
+        for s, e in zip(exStarts, exEnds):
+            cs, ce = max(s, cdsStart), min(e, cdsEnd)
+            if cs < ce: segs.append([cs, ce])
+        # cumulative CDS bases before each exon, in TRANSCRIPTION order
+        order = range(len(segs)) if strand == "+" else range(len(segs) - 1, -1, -1)
+        cum = 0
+        cumBefore = {}
+        for i in order:
+            cumBefore[i] = cum
+            cum += segs[i][1] - segs[i][0]
+        g = enst2gene.get(name)
+        if not g: continue
+        ensg, sym = g
+        tx[name] = (chrom, strand, ensg, sym,
+                    [(segs[i][0], segs[i][1], cumBefore[i]) for i in range(len(segs))])
+    return tx
+
+def writeCdsBed(tx, path):
+    with open(path, "w") as out:
+        for name, (chrom, strand, ensg, sym, segs) in tx.items():
+            for (s, e, cb) in segs:
+                # name field packs everything the compute step needs
+                out.write("%s\t%d\t%d\t%s|%s|%s|%s|%d\t0\t%s\n" %
+                          (chrom, s, e, name, ensg, sym, strand, cb, strand))
+
+def cdsOffset(strand, exonStart, exonEnd, cumBefore, repBase):
+    if strand == "+":
+        return cumBefore + (repBase - exonStart)
+    return cumBefore + (exonEnd - 1 - repBase)
+
+def selftest(tx, prot):
+    """Verify that translating each transcript codon-by-codon from the genome-
+    derived mapping reproduces the MANE protein for a sample of transcripts."""
+    import random
+    names = list(tx.keys()); random.seed(1); random.shuffle(names)
+    checked = bad = 0
+    for name in names[:200]:
+        chrom, strand, ensg, sym, segs = tx[name]
+        p = prot.get(ensg)
+        if not p: continue
+        # first CDS base in transcription order -> must be aa position 1
+        if strand == "+":
+            first = segs[0]; base = first[0]
+        else:
+            last = segs[-1]; base = last[1] - 1
+        # find that base's exon and compute
+        for (s, e, cb) in segs:
+            if s <= base < e:
+                off = cdsOffset(strand, s, e, cb, base)
+                if off != 0: bad += 1
+                break
+        checked += 1
+    sys.stderr.write("selftest: %d transcripts, %d with first-base!=aa1\n" % (checked, bad))
+    return bad == 0
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("clinvarBb"); ap.add_argument("maneGp")
+    ap.add_argument("geneFaa"); ap.add_argument("maneMeta")
+    ap.add_argument("-o", required=True); ap.add_argument("--selftest", action="store_true")
+    ap.add_argument("--workdir", default=None)
+    args = ap.parse_args()
+
+    enst2gene = {}
+    for line in open(args.maneMeta):
+        f = line.rstrip("\n").split("\t")
+        enst2gene[f[0]] = (f[1], f[2])   # enst -> (ensg, sym)
+
+    prot = loadFasta(args.geneFaa)
+    tx = parseGenePred(args.maneGp, enst2gene)
+    sys.stderr.write("loaded %d MANE transcripts, %d proteins\n" % (len(tx), len(prot)))
+
+    if args.selftest:
+        ok = selftest(tx, prot)
+        if not ok: sys.exit("SELFTEST FAILED")
+
+    wd = args.workdir or tempfile.mkdtemp()
+    os.makedirs(wd, exist_ok=True)
+    cdsBed = os.path.join(wd, "cdsExons.bed")
+    varBed = os.path.join(wd, "codingVars.bed")
+    writeCdsBed(tx, cdsBed)
+
+    # extract coding variants from the clinvar bigBed into a BED4+
+    sys.stderr.write("extracting coding variants ...\n")
+    nvar = 0
+    jsonByVid = {}   # variantId -> HGVS json (col28) for protChange display
+    with open(varBed, "w") as vb:
+        proc = subprocess.Popen(["bigBedToBed", args.clinvarBb, "stdout"],
+                                stdout=subprocess.PIPE, text=True)
+        for line in proc.stdout:
+            f = line.rstrip("\n").split("\t")
+            if f[C_MOLCONSEQ - 1] not in CODING: continue
+            jsonByVid[f[C_VARID - 1]] = f[27]
+            # name packs variant metadata; '|' join, sanitize stray '|'
+            meta = "%s;%s;%s;%s;%s" % (
+                f[C_VCV - 1], f[C_VARID - 1], f[C_CLINCODE - 1],
+                f[C_CLINSIGN - 1].replace(";", ","), f[C_STARS - 1])
+            vb.write("%s\t%s\t%s\tv%d\t%s\t%s\n" %
+                     (f[0], f[1], f[2], nvar, meta, f[C_MOLCONSEQ - 1].replace(" ", "_")))
+            nvar += 1
+        proc.wait()
+    sys.stderr.write("coding variants: %d\n" % nvar)
+
+    # sort both then intersect
+    subprocess.run("sort -k1,1 -k2,2n %s -o %s" % (cdsBed, cdsBed), shell=True, check=True)
+    subprocess.run("sort -k1,1 -k2,2n %s -o %s" % (varBed, varBed), shell=True, check=True)
+
+    sys.stderr.write("intersecting ...\n")
+    isect = subprocess.Popen(
+        ["bedtools", "intersect", "-a", varBed, "-b", cdsBed, "-wa", "-wb", "-sorted"],
+        stdout=subprocess.PIPE, text=True)
+
+    nrow = miss = 0
+    with open(args.o, "w") as out:
+        for line in isect.stdout:
+            f = line.rstrip("\n").split("\t")
+            vStart, vEnd = int(f[1]), int(f[2])
+            meta = f[4].split(";"); molConseq = f[5].replace("_", " ")
+            # b-record starts at column 6 (0-based): chrom s e name score strand
+            bStart, bEnd = int(f[7]), int(f[8])
+            nm = f[9].split("|")   # name enst ensg sym strand cumBefore
+            enst, ensg, sym, strand, cb = nm[0], nm[1], nm[2], nm[3], int(nm[4])
+            # representative base = transcription-first affected coding base in this exon
+            if strand == "+":
+                repBase = max(vStart, bStart)
+            else:
+                repBase = min(vEnd - 1, bEnd - 1)
+            if repBase < bStart or repBase >= bEnd:
+                continue
+            off = cdsOffset(strand, bStart, bEnd, cb, repBase)
+            aaPos = off // 3 + 1
+            p = prot.get(ensg, "")
+            refRes = p[aaPos - 1] if 0 < aaPos <= len(p) else "?"
+            if refRes == "?": miss += 1
+            vcv, varId, clinCode, clinSign, stars = meta[0], meta[1], meta[2], meta[3], meta[4]
+            protChange = protChangeForPos(jsonByVid.get(varId, ""), aaPos)
+            out.write("\t".join([ensg, enst, sym, str(aaPos), refRes, f[0],
+                                 str(repBase), vcv, varId, clinCode, clinSign,
+                                 stars, molConseq, protChange]) + "\n")
+            nrow += 1
+    isect.wait()
+    sys.stderr.write("wrote %d variant-codon rows (%d with out-of-range aaPos)\n" % (nrow, miss))
+
+if __name__ == "__main__":
+    main()