fcf788d7357f6f207f2bb0b61acde9dc0507d89f lrnassar Mon Sep 21 04:04:10 2026 -0700 Fix two stale figures and guard the accession interpolation in the MaveMD build, per CR. refs #37800 The makeDoc script listing still said the variant converter writes bed12+31. It writes bed12+34, as the autoSql, mavemd.ra, runBuild.sh and the makeDoc's own build-results line all already said. mavemdLib.py's docstring claimed the codon projection is cross-checked against "~154k variants that carry both a genomic and a protein term". 154,895 is the number placed by the genomic route; the cross-check set is the smaller number that also has a resolvable protein term. The docstring now describes the set rather than quoting a figure that drifts with every build. Also adds checkAccession() and calls it before either query that interpolates an accession into SQL. Nothing can currently reach those queries with a quote in it, since the accessions come from PROTEIN_TERM whose character class excludes one, but the regex is a hundred lines from the query and a later edit to it should not be able to open this up silently. Output is byte-identical to the previous build. diff --git src/hg/makeDb/scripts/mavemd/mavemdLib.py src/hg/makeDb/scripts/mavemd/mavemdLib.py index 7c16157c30b..24b8102165a 100644 --- src/hg/makeDb/scripts/mavemd/mavemdLib.py +++ src/hg/makeDb/scripts/mavemd/mavemdLib.py @@ -1,30 +1,33 @@ #!/usr/bin/env python3 """Shared helpers for the MaveMD track build: coordinate projection and colors. 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. It is validated in makeMaveMdVariants.py against the ~154k variants -that carry both a genomic and a protein term, so a drift in either RefSeq or MaveDB's +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 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': '*', } # Two palettes, because the track carries two different kinds of statement. @@ -84,30 +87,45 @@ # Order used only when a score set has no MaveDB-designated primary calibration and we # have to say which of several calls to show. Strongest pathogenic first, then strongest # benign, with "not met" last because it is the absence of evidence either way. # # This list ranks every pathogenic call above every benign one, so if two calibrations ever # disagree in direction the pathogenic one wins regardless of strength. No variant in the # collection currently has calibrations that disagree in direction, so the bias is latent; # if one appears, this ordering is the thing to revisit. ACMG_SEVERITY = ['PS3_very_strong', 'PS3', 'PS3_moderate_plus', 'PS3_moderate', 'PS3_supporting', 'BS3_very_strong', 'BS3', 'BS3_moderate_plus', 'BS3_moderate', 'BS3_supporting', 'PS3_not_met', 'BS3_not_met'] +# Accessions reach hgsql() by string interpolation, so they are checked against this +# first. Today they can only arrive via PROTEIN_TERM in makeMaveMdVariants.py, whose +# character class already excludes quotes, but that regex is far from the query and a +# future edit to it should not be able to open this up silently. +SAFE_ACCESSION = re.compile(r'^[A-Za-z0-9_.]+$') + + +def checkAccession(acc): + """Fail loudly on an accession that has no business being pasted into SQL.""" + if not SAFE_ACCESSION.match(acc): + raise ValueError('refusing to query on accession %r: expected letters, digits, ' + 'underscore and dot only' % acc) + return acc + + def hgsql(db, query): """Run a query and return rows as lists of strings.""" out = subprocess.run(['hgsql', db, '-N', '-e', query], check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, universal_newlines=True) return [line.split('\t') for line in out.stdout.rstrip('\n').split('\n') if line] def isMainChrom(chrom): """True for chr1..chr22, chrX, chrY, chrM - not alts, randoms or patches.""" return '_' not in chrom GENCODE_ATTRS = 'wgEncodeGencodeAttrsV50' GENCODE_GENEPRED = 'wgEncodeGencodeCompV50' @@ -115,30 +133,31 @@ def loadProteinToTranscript(db, protAccs): """Map protein accessions to their transcripts. MaveDB states protein terms against RefSeq (NP_) for most score sets and against Ensembl (ENSP) for a handful, so both routes are needed: NP_ through ncbiRefSeqLink, ENSP through the GENCODE attributes table. Each tries the exact versioned accession first, then any version of the same base accession. Returns (mapping, unresolved) so the caller can report a whole score set going missing rather than silently dropping it. """ mapping = {} unresolved = [] for acc in protAccs: + checkAccession(acc) base = acc.split('.')[0] if acc.startswith('ENSP'): table, col, key = GENCODE_ATTRS, 'transcriptId', 'proteinId' else: table, col, key = 'ncbiRefSeqLink', 'mrnaAcc', 'protAcc' rows = hgsql(db, "select %s from %s where %s = '%s'" % (col, table, key, acc)) if not rows: rows = hgsql(db, "select %s from %s where %s like '%s.%%'" % (col, table, key, base)) if rows: mapping[acc] = rows[0][0] else: unresolved.append(acc) return mapping, unresolved @@ -202,30 +221,31 @@ runs.append((runStart, runLen)) return max(runs, key=lambda r: r[1]) def loadCodonMaps(db, transcripts): """Build a CodonMap for each transcript from ncbiRefSeqCurated. RefSeq transcripts come from ncbiRefSeqCurated and Ensembl ones from the GENCODE genePred, keyed off the accession prefix. A transcript can align to more than one place (alt haplotypes, fix patches); the alignment on a main chromosome wins, and among several the longest CDS wins. """ maps = {} missing = [] for tx in transcripts: + checkAccession(tx) table = GENCODE_GENEPRED if tx.startswith('ENST') else 'ncbiRefSeqCurated' rows = hgsql(db, "select chrom, strand, cdsStart, cdsEnd, exonStarts, exonEnds " "from %s where name = '%s'" % (table, tx)) best = None for chrom, strand, cdsStart, cdsEnd, exonStarts, exonEnds in rows: cdsStart, cdsEnd = int(cdsStart), int(cdsEnd) starts = [int(x) for x in exonStarts.rstrip(',').split(',')] ends = [int(x) for x in exonEnds.rstrip(',').split(',')] bases = [] for s, e in zip(starts, ends): s = max(s, cdsStart) e = min(e, cdsEnd) if s < e: bases.extend(range(s, e)) if not bases: