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/makeMaveMdHeatmap.py src/hg/makeDb/scripts/mavemd/makeMaveMdHeatmap.py index 2cc579d17de..e78fc42ace9 100755 --- src/hg/makeDb/scripts/mavemd/makeMaveMdHeatmap.py +++ src/hg/makeDb/scripts/mavemd/makeMaveMdHeatmap.py @@ -34,38 +34,40 @@ FALLBACK_BOUNDS = '0,1' FALLBACK_COLORS = '#f7f7f7,#b2182b' def assayLine(meta): """One line naming what a score set measured, for the legend and the cell mouseovers. Two maps of the same gene routinely disagree because they measured different things: PTEN abundance against PTEN lipid phosphatase activity, GCK activity against GCK abundance, KCNE1 trafficking with and without KCNQ1. A reader cannot make sense of that without knowing which assay they are looking at, so the assay travels with the map rather than sitting a click away on the details page. Method and model system come first because they are short and always present; the score set title can be long and is the part that gets truncated. - The separator is a plain hyphen, not a middot: the legend is drawn as raster text by - hgTracks, so an HTML entity from bedField() would appear literally as "·". + The whole line is flattened to ASCII, separator included: the legend is drawn as raster + text by hgTracks, so anything bedField() turns into an HTML entity appears literally on + the image. MaveDB titles carry en dashes ("BRCA2 exons 15-26", "CARD11 exons 3-5"), so + the dash in the title matters as much as the one between the fields. """ method = meta.get('assayMethod') or '' model = meta.get('assayModel') or '' title = meta.get('title') or '' head = '%s in %s' % (method, model) if method and model else (method or model) - return ' - '.join(p for p in (head, title) if p) + return lib.asciiText(' - '.join(p for p in (head, title) if p)) def severityRank(cell): """Rank a cell's call so the strongest evidence wins a tie. Lower is stronger.""" code = cell.get('outcome') if code in lib.ACMG_SEVERITY: return lib.ACMG_SEVERITY.index(code) order = {'abnormal': 0, 'normal': 1, 'indeterminate': 2} return len(lib.ACMG_SEVERITY) + order.get(cell.get('funcClass'), 3) def direction(cell): """'path', 'benign' or '' for a cell's call, ignoring strength.""" code = cell.get('outcome') or '' if code and not code.endswith('_not_met'): @@ -91,48 +93,53 @@ call = (cell.get('calls') or {}).get(calUrn) if call is None: stats['cellOutsideChosenCalibration'] += 1 continue cell['outcome'] = call['acmgOutcome'] cell['funcClass'] = call['funcClass'] cell['color'] = lib.cellColor(call['acmgOutcome'], call['funcClass']) def main(): parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) parser.add_argument('downloadDir') parser.add_argument('outBed') parser.add_argument('--db', default='hg38') + parser.add_argument('--twoBit', default='/hive/data/genomes/hg38/hg38.2bit', + help='genome sequence, for the wild-type residue check') + parser.add_argument('--workDir', default='.', + help='scratch directory for the sequence fetch') parser.add_argument('--classPalette', default='purple', choices=sorted(lib.CLASS_PALETTES), help='palette for measurements with no ACMG code') args = parser.parse_args() lib.setClassPalette(args.classPalette) scoreSets = loadScoreSets(args.downloadDir) # Resolve every protein accession once. protAccs = set() for path in sorted(glob.glob(os.path.join(args.downloadDir, 'variants', '*.csv'))): with open(path, newline='') as fh: for row in csv.DictReader(fh): match = PROTEIN_TERM.match(clean(row.get('mavedb.post_mapped_hgvs_p')) or '') if match: protAccs.add(match.group('acc')) protToTx, unresolved = lib.loadProteinToTranscript(args.db, sorted(protAccs)) codonMaps, missing = lib.loadCodonMaps(args.db, sorted(set(protToTx.values()))) + lib.addProteinSequence(args.db, codonMaps, args.twoBit, args.workDir) if unresolved: sys.stderr.write(" WARNING: no transcript for %s\n" % ', '.join(unresolved)) if missing: sys.stderr.write(" WARNING: no genePred for %s\n" % ', '.join(missing)) stats = collections.Counter() entries = [] for path in sorted(glob.glob(os.path.join(args.downloadDir, 'variants', '*.csv'))): with open(path, newline='') as fh: reader = csv.DictReader(fh) calCols = calibrationColumns(reader.fieldnames) clinvarRelease = '' for name in reader.fieldnames: if name.startswith('clinvar.') and name.endswith('.clinical_significance'): @@ -154,30 +161,39 @@ match = PROTEIN_TERM.match(protein) if protein else None if not match: stats['skipNoProteinTerm'] += 1 continue tx = protToTx.get(match.group('acc')) thisMap = codonMaps.get(tx) if tx else None if thisMap is None: stats['skipNoCodonMap'] += 1 continue if codonMap is None: codonMap = thisMap elif thisMap.tx != codonMap.tx: stats['skipOtherTranscript'] += 1 continue protPos = int(match.group('pos')) + # The column is placed from the protein term, so a term whose numbering + # disagrees with the genome would put the cell one codon off. The variant + # track can fall back on MaveDB's genomic mapping for these; a map has no + # such fallback, so the cell is dropped instead of drawn in the wrong place. + wtOne = lib.THREE_TO_ONE.get(match.group('wt')) + if (codonMap.protein and wtOne and protPos <= len(codonMap.protein) + and codonMap.protein[protPos - 1] != wtOne): + stats['skipResidueMismatch'] += 1 + continue block = codonMap.codonBlock(protPos) if block is None: stats['skipPositionPastCds'] += 1 continue # The column is drawn on the longest contiguous run of the codon's bases, so # a codon split across an exon junction keeps its block on coding sequence # instead of starting inside the intron. blockStart, blockLen = block bases = range(blockStart, blockStart + blockLen) wt = lib.THREE_TO_ONE.get(match.group('wt'), match.group('wt')) rawVar = match.group('var') aa = wt if rawVar == '=' else lib.THREE_TO_ONE.get(rawVar, rawVar) if aa not in lib.HEATMAP_ROWS: stats['skipNonStandardResidue'] += 1