824df26b6320b692d629566c5a10b15004da82ce lrnassar Tue Sep 29 16:09:00 2026 -0700 addProteinSequence in mavemdLib translates each transcript's CDS from hg38.2bit so makeMaveMdVariants can check every projected codon against the reference residue its own HGVS term asserts; the existing comparison against MaveDB's genomic mapping only reaches the 3% of projected items that carry both terms, because 18 of the 40 protein accessions have no genomic-route variants at all. 39 of 40 accessions match at 0.000%; NP_689629.2 (FKRP) has 99 nonsense terms numbered one codon downstream of their own reference residue, which still reach mavemdVar through MaveDB's genomic mapping but are dropped from mavemdMap, which places columns from the protein term and has no fallback. The haplotype test now also reads hgvs_nt, since PTEN 00000054-a-1 states 1,236 haplotypes as c.[1207G>T;1209C>T] with no protein term and they were counted as rejected submissions, making both figures in the makeDoc wrong. assayLine runs the heatmap legend through asciiText because bedField turns the en dash in three MaveDB titles into – and the legend is drawn as raster text; clinGenId links to by_canonicalid rather than /allele, which serves JSON to a browser, matching human/civic.ra; the generated filter fragment no longer emits the blank line after each group that the makeDoc itself warns ends a stanza; and runBuild.sh tails the log on failure instead of dying silently under set -e. Also reworded the grey legend entry, which said no threshold was reached in either direction but covers 6,656 normal and 145 abnormal items, alphabetized the references, and fixed stale counts in the makeDoc. Caught by Claude review of 29af14b, fcf788d and 97c7de5. refs #38407 refs #37800 diff --git src/hg/makeDb/scripts/mavemd/mavemdLib.py src/hg/makeDb/scripts/mavemd/mavemdLib.py index 24b8102165a..40464c30f3f 100644 --- src/hg/makeDb/scripts/mavemd/mavemdLib.py +++ src/hg/makeDb/scripts/mavemd/mavemdLib.py @@ -3,30 +3,31 @@ MaveDB's VRS mapper resolves some score sets all the way to the genome and others only to a protein sequence. About a third of MaveMD variants arrive with a genomic HGVS term and can be placed directly; the rest carry only an NP_ protein term and have to be walked back to their codon here. The projection is: NP_ accession -> NM_ transcript (hg38 ncbiRefSeqLink) -> genePred (hg38 ncbiRefSeqCurated) -> the CDS bases in transcription order -> the three genomic bases of codon N. makeMaveMdVariants.py validates it against every variant that carries both a genomic term and a resolvable protein term, so a drift in either RefSeq or MaveDB's mapper shows up as a coordinate disagreement rather than as silently wrong placements. That set is smaller than the genomic-route placement count, because some variants have a genomic term and no protein term; the makeDoc records both figures per build. """ +import os import re import subprocess import sys # Standard amino acids ordered by class, matching the MaveDB and popEVE heatmap tracks, # with Ter appended as a final row (MaveMD carries ~17k nonsense measurements). STANDARD_AAS = list('AVLIMFYWRHKDESTNQGCP') HEATMAP_ROWS = STANDARD_AAS + ['*'] THREE_TO_ONE = { 'Ala': 'A', 'Arg': 'R', 'Asn': 'N', 'Asp': 'D', 'Cys': 'C', 'Gln': 'Q', 'Glu': 'E', 'Gly': 'G', 'His': 'H', 'Ile': 'I', 'Leu': 'L', 'Lys': 'K', 'Met': 'M', 'Phe': 'F', 'Pro': 'P', 'Ser': 'S', 'Thr': 'T', 'Trp': 'W', 'Tyr': 'Y', 'Val': 'V', 'Ter': '*', } @@ -165,30 +166,31 @@ class CodonMap(object): """Genomic coordinates of every codon of one transcript. codons[n] is the list of three 0-based genomic positions of codon n (1-based protein numbering), in transcription order. For a codon split across an intron those three positions are not contiguous, which is why callers take min() and clamp block widths rather than assuming a 3bp run. """ def __init__(self, tx, chrom, strand, cdsBases): self.tx = tx self.chrom = chrom self.strand = strand self.cdsBases = cdsBases self.protLen = len(cdsBases) // 3 + self.protein = None # filled in by addProteinSequence() def codon(self, protPos): """Genomic positions of codon protPos (1-based), or None if out of range.""" i = (protPos - 1) * 3 if protPos < 1 or i + 3 > len(self.cdsBases): return None return self.cdsBases[i:i + 3] def codonSpan(self, protPos): """Full genomic span of a codon as (start0, end), or None. A codon at an exon junction is split, so its three bases are not a 3bp run; the span then covers the intervening intron. That is the honest extent of the codon, and it is what a protein-level measurement actually pins down. """ @@ -255,46 +257,157 @@ cand = CodonMap(tx, chrom, strand, bases) if best is None: best = cand elif isMainChrom(cand.chrom) and not isMainChrom(best.chrom): best = cand elif isMainChrom(cand.chrom) == isMainChrom(best.chrom) and \ len(cand.cdsBases) > len(best.cdsBases): best = cand if best is None: missing.append(tx) else: maps[tx] = best return maps, missing +DASHES = {0x2010: '-', 0x2011: '-', 0x2012: '-', 0x2013: '-', 0x2014: '-', 0x2015: '-', + 0x2212: '-', 0x00ad: '-'} + + +def asciiText(value): + """Flatten text to ASCII, for fields hgTracks draws as raster rather than HTML. + + bedField() turns non-ASCII into numeric HTML entities, which is right for the details + page and the mouseovers but wrong for the heatmap legend: that is drawn with a graphics + library, so an entity appears literally as "–" on the image. MaveDB score set + titles carry en dashes ("BRCA2 exons 15-26"), so anything bound for the legend comes + through here first. + """ + text = '' if value is None else str(value) + out = [] + for ch in text: + code = ord(ch) + if code < 128: + out.append(ch) + elif code in DASHES: + out.append('-') + else: + import unicodedata + folded = unicodedata.normalize('NFKD', ch).encode('ascii', 'ignore').decode() + out.append(folded) + return ''.join(out) + + def bedField(value): """Render one BED field: no tabs, no newlines, no non-ASCII. MaveDB free text (score set titles, formatted citations) carries all three. Tabs and newlines would split the row. Non-ASCII is subtler: the browser does not transcode UTF-8, so an en dash or an author name like Gr\u00f8nb\u00e6k-Thygesen reaches the details page as mojibake. Numeric HTML entities render correctly instead. """ text = '' if value is None else str(value) for ch in ('\t', '\n', '\r'): text = text.replace(ch, ' ') while ' ' in text: text = text.replace(' ', ' ') return ''.join(c if ord(c) < 128 else '&#%d;' % ord(c) for c in text.strip()) +CODON_TABLE = {} +for _bases, _aas in ( + ('TTT TTC', 'F'), ('TTA TTG CTT CTC CTA CTG', 'L'), ('ATT ATC ATA', 'I'), + ('ATG', 'M'), ('GTT GTC GTA GTG', 'V'), ('TCT TCC TCA TCG AGT AGC', 'S'), + ('CCT CCC CCA CCG', 'P'), ('ACT ACC ACA ACG', 'T'), ('GCT GCC GCA GCG', 'A'), + ('TAT TAC', 'Y'), ('TAA TAG TGA', '*'), ('CAT CAC', 'H'), ('CAA CAG', 'Q'), + ('AAT AAC', 'N'), ('AAA AAG', 'K'), ('GAT GAC', 'D'), ('GAA GAG', 'E'), + ('TGT TGC', 'C'), ('TGG', 'W'), ('CGT CGC CGA CGG AGA AGG', 'R'), + ('GGT GGC GGA GGG', 'G')): + for _b in _bases.split(): + CODON_TABLE[_b] = _aas + +COMPLEMENT = str.maketrans('ACGTacgtNn', 'TGCAtgcaNn') + + +def addProteinSequence(db, maps, twoBit, workDir): + """Translate each transcript's CDS from the genome and hang it on its CodonMap. + + This is what lets makeMaveMdVariants.py check a projected codon against the reference + residue the HGVS term asserts. The cross-check against MaveDB's own genomic mapping + cannot do that job: 18 of the 40 protein accessions have no genomic-route variants at + all, and they hold 97% of the projected items, so a shifted CDS in any of them would + move every item and still pass. + """ + bedPath = os.path.join(workDir, 'cdsExons.bed') + faPath = os.path.join(workDir, 'cdsExons.fa') + runs = [] + with open(bedPath, 'w') as fh: + for cm in maps.values(): + start = None + prev = None + for pos in sorted(cm.cdsBases): + if start is None: + start = prev = pos + elif pos == prev + 1: + prev = pos + else: + fh.write('%s\t%d\t%d\t%s:%d-%d\n' + % (cm.chrom, start, prev + 1, cm.chrom, start, prev + 1)) + runs.append((cm.chrom, start, prev + 1)) + start = prev = pos + if start is not None: + fh.write('%s\t%d\t%d\t%s:%d-%d\n' + % (cm.chrom, start, prev + 1, cm.chrom, start, prev + 1)) + runs.append((cm.chrom, start, prev + 1)) + + # -bedPos makes the fasta id chrom:start-end; the bed needs a name column either way + subprocess.run(['twoBitToFa', '-bed=' + bedPath, '-bedPos', twoBit, faPath], + check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + + seq = {} + name = None + chunks = [] + with open(faPath) as fh: + for line in fh: + if line.startswith('>'): + if name: + seq[name] = ''.join(chunks) + name = line[1:].strip() + chunks = [] + else: + chunks.append(line.strip()) + if name: + seq[name] = ''.join(chunks) + + base = {} + for chrom, start, end in runs: + s = seq.get('%s:%d-%d' % (chrom, start, end)) + if s is None: + continue + for offset, ch in enumerate(s): + base[(chrom, start + offset)] = ch + + for cm in maps.values(): + letters = [] + for pos in cm.cdsBases: + ch = base.get((cm.chrom, pos), 'N') + letters.append(ch.translate(COMPLEMENT) if cm.strand == '-' else ch) + cds = ''.join(letters).upper() + cm.protein = ''.join(CODON_TABLE.get(cds[i:i + 3], 'X') + for i in range(0, len(cds) - 2, 3)) + + def fmtScore(value, places=4): """Format a functional score for display, trimming trailing zeros.""" try: text = ('%.*f' % (places, float(value))).rstrip('0').rstrip('.') except (TypeError, ValueError): return '' return '0' if text in ('', '-0') else text def pickDisplayCall(calls, primaryTitle=None): """Choose which calibration's call drives the color and the filters. MaveDB curates a `primary` flag on score calibrations, with its own promote and demote API endpoints, so where that primary actually classifies the variant it is their editorial choice and we use it as-is.