# [Claude/lrnassar] MaveMD - clinically calibrated multiplexed variant effect measurements (2026-09-17)
# refs #37800

# MaveMD ("MAVEs for MeDicine") is not a separate database. It is a curated collection inside
# MaveDB holding the score sets that carry clinical evidence calibrations, so everything comes
# from the ordinary MaveDB API. The collection URN is
#   urn:mavedb:collection-603dafbf-4a3f-4d70-ab8c-aafb226fbff4
# which is not documented; it was read out of the MaveDB front end, which fetches it to render
# the mavedb.org/mavemd page.
# Reference: McEwen et al. (2025) medRxiv, PMID 41332838. A local copy of the paper is in
# paper/ alongside this build:
#   mkdir -p paper && cd paper
#   curl -s "https://www.ebi.ac.uk/europepmc/webservices/rest/PMC12668102/fullTextXML" -o paper.xml
#   # paper.md is that XML flattened to markdown (title, abstract, sections) for grepping
#
# Two things from the paper shaped the build:
#  1. The paper's counts are a snapshot the collection has since grown past: it reports 74
#     datasets over 32 genes and 438,318 measurements, against 84 score sets drawn from 81
#     experiments over 39 genes and 481,717 measurements here. The paper does not define how
#     a "dataset" maps onto a score set or an experiment, and neither 84 nor 81 is 74, so do
#     not try to reconcile them. Say "score set" in the docs, not "dataset".
#  2. Counts are measurement-level. A variant measured by several assays contributes several
#     measurements. Our 462,486 items are 275,791 distinct ClinGen alleles, which is why a
#     codon studied by several groups stacks over a hundred items deep.

# Request pacing is what MaveDB asked us for directly (Benjamin Capodanno, 21 Aug 2026 email):
# metadata at no more than two concurrent requests, score and variant endpoints sequential,
# page at 100k variants, and back off rather than retry immediately on a 504, because the
# gateway gives up while their query keeps running. fetchMaveMd.py implements all of this.

mkdir -p /hive/data/outside/mavemd/2026-09-17
mkdir -p /hive/data/genomes/hg38/bed/mavemd/2026-09-17

# Download: the collection, one metadata record per score set, and one flat CSV per score set
# from score-sets/{urn}/variants/data. The CSV joins the assay score to the post-mapped HGVS
# terms, ClinGen allele IDs, gnomAD frequencies, the most recent ClinVar release, and every
# calibration's functional class and ACMG call. All calibrations are requested, research-use-only
# ones included; only the newest of the twelve available ClinVar releases is requested.
~/kent/src/hg/makeDb/scripts/mavemd/fetchMaveMd.py /hive/data/outside/mavemd/2026-09-17

#   MaveDB API 2026.2.7.2
#   Collection MaveMD: 84 score sets, modified 2026-04-30
#   Downloaded 242.1 MB in 257 seconds

# Build both tracks: the per-variant bigBed and the per-score-set heatmap.
~/kent/src/hg/makeDb/scripts/mavemd/runBuild.sh \
    /hive/data/outside/mavemd/2026-09-17 \
    /hive/data/genomes/hg38/bed/mavemd/2026-09-17

# Scripts (all in ~/kent/src/hg/makeDb/scripts/mavemd/):
#   fetchMaveMd.py         download the collection from the API
#
#   mavemdLib.py           calibration choice and color tables; also puts
#                          utils/otto/mavedb on sys.path for mavedbPlace
#   makeMaveMdVariants.py  per-variant bed12+35 -> mavemdVar.bb
#   makeMaveMdHeatmap.py   per-score-set bed12+20 -> mavemdMap.bb
#   mavemdVariants.as, mavemdHeatmap.as
#   runBuild.sh            driver for the two converters plus bedToBigBed

# COORDINATES. Placement is not in this directory: it is mavedbPlace.py, in
# ~/kent/src/hg/utils/otto/mavedb, shared with the MaveDB track (#38396) so that a variant
# lands at the same base in both. mavemdLib.py puts that directory on sys.path, so the two
# converters can "import mavedbPlace as mp"; nothing else is needed to run them by hand.
#
# The rule is one offset per score set target: reference position = experiment position +
# offset, chosen by testing the wild-type residues the experiment itself states against the
# reference protein. MaveDB supplies an offset of its own and it wins wherever it fits those
# residues exactly, or at every position but one (a construct carrying a known change, see
# below). Otherwise our own search over the gene's RefSeq proteins finds the offset that does
# fit. Every placed row is then checked against the reference protein and against the genome
# codon, and MaveDB's genomic term, where there is one, must fall inside that codon.
#
# Placed, 462,486 of 481,717 measurements (96.0%), by where the offset came from:
#   430,133  MaveDB, at the position it mapped
#    18,513  MaveDB's offset did not fit the experiment's residues; ours did
#     6,519  MaveDB named another isoform of the same protein
#     7,321  MaveDB did not map it at all; the experiment's own term is on a reference
#            sequence (NM_/NC_), placed through kent's hgvsToVcf with a reference-base check
#
# THREE SCORE SETS MAVEDB MAPS AT THE WRONG CODON. This is why the offset is checked rather
# than taken. In each the experiment's stated residues fit the reference protein at exactly
# one offset, and it is not the one in MaveDB's mapped terms. The mapped terms are not a
# clue: MaveDB rewrites the reference residue in them to agree with the codon it chose, so
# they are self-consistent and agree with the genome, and the error is only visible from the
# experiment's own residues. Until 2026-10 this track drew all of these at MaveDB's position.
#   00000663-a-1  NDUFAF6, NP_689629.2. All 5,714 stated residues fit at offset 0; MaveDB
#                 mapped 3,668 rows one codon downstream. p.Ser45Phe was drawn as V46S.
#   00001251-a-1  TPK1, NP_071890.2. All fit at offset 0; MaveDB mapped 2,789 rows 49 codons
#                 downstream, 49 being the length of the segment missing from the isoform
#                 deposited as the target sequence. p.Val195Thr was drawn as *244T, at the
#                 stop codon.
#   00000060-a-1  TARDBP, NP_031401.1, and -a-2. The constructs are the TDP-43 prion-like
#                 domain, AA 290-331 and AA 332-373, numbered 1-42; MaveDB read 1-42 as
#                 full-protein positions, so both score sets drew at codons 1-42. Offsets
#                 +289 and +331.
# Reported to the MaveDB authors 2026-10; see the draft in /hive/users/lrnassar/claude/RM37800.
#
# Our own search also wins on 00000665-a/b/c/d-1 (HBB, NP_000271.1), at offset 0, where
# MaveDB mapped nothing to disagree with.
#
# Score sets numbered on a domain rather than the whole protein, where MaveDB's offset fits
# the stated residues exactly: 00000098-a-2 (+1619), 00001222-a-2 (+1672), 00001229-a-2
# (+544), 00001267-a-2 (+553) and -b-1 (+1511).
#
# CONSTRUCTS THAT DIFFER FROM THE REFERENCE. Eleven score sets use a construct carrying one
# known change of its own, and the offset is accepted with that one position excused. The
# measurements at that position are of a different starting residue from the reference, so
# they are dropped rather than drawn as though they were the reference variant: 221 rows in
# all. TP53 00000068-0/a/b/c-1, 00001235-a-1 and 00001236-0-1 at residue 72 (the common P72R
# form, which is in the target name of three of them); KCNE1 00000674-a-2/b-1/c-1 at 38;
# BRCA1 00001222-b-2 at 1613 (offset +1576); TPK1 00001251-a-1 at 193.
# The other 68,095 items of those eleven score sets are placed normally, since the
# construct matches the reference everywhere else, and each carries the construct field:
# "differs from NP_000537.3 at position 72 (Pro in the reference, Arg in the construct)"
# and its three counterparts. It is OffsetChoice.constructNote() in mavedbPlace.py, the
# same string the MaveDB track prints, so a reader can see which backbone a score was
# measured on. 48,742 of those items are TP53, 9,518 TPK1, 7,574 KCNE1, 2,261 BRCA1.
#
# 4,474 items sit in a codon split across an exon junction, so its three bases are not a 3 bp
# run. Those are written as BED12 with the two real blocks, so the intron draws as a thin
# connector rather than as a solid item up to 133 kb wide. The heatmap column is drawn on the
# longest contiguous run of the codon's bases (CodonMap.codonBlock), keeping the block off
# the intron.
#
# An item placed on a protein covers its codon, three bases, even where MaveDB also gave a
# genomic term naming the one changed base. Before the switch to the shared code those rows
# were drawn one base wide, which is why basesCovered went from 667,558 to 1,272,893 with
# almost no change in item count. The codon is what the measurement is of, and it is what
# the heatmap column covers, so the two tracks now agree.
#
# CHECKS. All of these are in mavedbPlace.py and run on every build:
#   - every placed row's wild-type residue must match the genome translation of its codon;
#     an accession failing on more than 25% of its rows is dropped whole (0 accessions)
#   - MaveDB's genomic term, where the row also has a protein placement, must fall inside
#     that codon. 124,335 rows carry both; 0 disagreed. crossCheckGate() stops the build
#     above 0.5%, which is what a drift in RefSeq or in MaveDB's mapper would look like.
#   - a reference-term placement must match the reference base hgvsToVcf reports
#
# GENCODE: the six ENSP-stated score sets (33,201 items) resolve through the newest
# wgEncodeGencodeAttrsV* table, found at run time by mavedbPlace.newestGencode(), so a new
# GENCODE on hg38 needs no edit here. The RefSeq route (ncbiRefSeqLink, ncbiRefSeqCurated,
# and ncbiRefSeqLinkHistorical/ncbiRefSeqHistorical for a retired NP_ version) is unversioned.
#
# NOT PLACED (19,231 of 481,717):
#   18,963  more than one residue or site changed: several substitutions measured as one
#           unit, with no single position. Almost all from urn:mavedb:00000665-a/b/c/d-1.
#           A bracketed nucleotide term whose bases all fall in one codon, such as the PTEN
#           c.[1207G>T;1209C>T], is one amino acid substitution and is placed; the judgement
#           is made on the protein term where there is one (isMultiVariant).
#      221  at the one position where the construct differs from the reference, see above
#       28  MaveDB did not map it and its own term is not on a reference sequence
#       13  frameshift: the protein term names the codon where the stop appears, not the
#           site of the change, so it is not a position we can draw
#        4  multi-residue delins
#        2  protein term could not be parsed

# HEATMAP. All 84 score sets produce a map. Rows are the 20 standard amino acids ordered by
# class (A V L I M F Y W R H K D E S T N Q G C P), matching the MaveDB and popEVE heatmap
# tracks, plus a final row for nonsense. A synonymous measurement fills the reference
# residue's cell.
# The columns come from the same placements as the variant track's items, so the two cannot
# disagree about where a measurement belongs; a cell exists wherever a protein-placed item
# does and the substitution is one of the 21 rows.
#   432,754  protein-placed measurements (462,486 items less the 29,732 placed at genomic
#            level only, which have no residue to put in a row)
#      460  dropped: the substitution is not one of the 21 rows (Sec, Xaa and friends)
#   391,758  cells drawn. 40,536 measurements are a later nucleotide change encoding a
#            substitution already measured at that position. The strongest evidence wins the
#            cell (CSV order is not a defined rule and flipped hundreds of cells between
#            pathogenic and benign); the cell is marked '+', or '!' where the measurements
#            disagree in direction. 1,057 cells carry '!' and 27,932 carry '+'. The
#            cellConflictingDirection counter in the log reads 1,281 because it counts
#            pairwise disagreements, so a cell with three conflicting measurements scores
#            more than once; it is not a cell count.
#
# One calibration colors a whole map, or its cells are not comparable to each other. Of the
# score sets with no usable MaveDB primary, 16 have rows that fall to more than one
# calibration; for those the one that supplied the strongest call on the most rows wins and
# every cell is recolored to it. That is a vote over rows, not a count of how many variants
# each calibration classifies: in the competing score sets the rivals usually classify an
# identical number, and 00001234-c-1 picks the opposite calibration from its seven siblings
# on the strength tiebreak alone. (55 score sets reach the fallback in total: 20 with no primary at all and
# 35 whose primary classifies nothing. One further score set, urn:mavedb:00001267-a-2, carries
# no calibration of any kind.)
#
# The renderer splits the score array with chopCommas (keeps a trailing empty field) but the
# label array with chopByCharRespectDoubleQuotesKeepEmpty (drops one). If the last cell of the
# last row is empty the counts disagree and the track aborts, so the converter puts a
# placeholder label on that one cell. 42 score sets needed it. Same trap as popEVE.

# ASSAY METADATA. The experiment record's controlled keywords carry the fields McEwen et al.
# argue a clinician judges an assay by, so four of them ride on every item: Phenotypic Assay
# Method, Phenotypic Assay Model System, Molecular Mechanism Assessed, Variant Library Creation
# Method. The last matters most: an in vitro construct library introduces a synthetic copy of
# the target, so it cannot see an effect on splicing or NMD and can read falsely normal for
# one; an endogenous locus library can. 61 score sets use construct libraries, 16 endogenous,
# 3 state "N/A (meta-analysis)" and 4 have no keyword record at all. Molecular Mechanism
# Assessed is the sparsest of the four, present on 45 of 84. The keyword text is in
# keyword.label; there is no keyword.value in the payload.
#
# Not exposed by the API: the paper's per-assay "can this assay detect splicing / NMD
# variants" curation, and the rest of the ~180 curated fields. Only the controlled keywords
# above come through, so the library route is the closest available proxy.
#
# NON-ASCII. MaveDB free text carries en dashes and accented author names (Gronbaek-Thygesen
# and Olvera-Leon are spelled with o-slash, ae and o-acute upstream). The browser does not
# transcode UTF-8, so these reach the details page as
# mojibake. mavemdLib.bedField() converts every non-ASCII character to a numeric HTML entity;
# 35,695 variant rows and 6 heatmap rows needed it. Re-check after any rebuild with:
#   LC_ALL=C grep -cP '[^\x00-\x7F]' mavemdVar.bed mavemdMap.bed     # must be 0

# EMPTY FIELDS IN THE MOUSEOVER. The trackDb mouseOver pattern substitutes ${field}
# unconditionally, so an empty value leaves a dangling label ("ACMG evidence:" with nothing
# after it). Three fields can legitimately be empty, so the converter writes a readable
# placeholder rather than an empty string: acmgOutcome "none assigned" (63,338 items, score
# sets that classify variants without assigning an ACMG code), funcClass "not classified"
# (3,896 items, all from urn:mavedb:00001267-a-2, which carries no calibration at all) and
# funcScore "not reported" (5,124 items MaveDB reports as NA, all of them indeterminate).
# The placeholders are set after the color and the BED score are computed from the raw
# values, so neither changes, and both label fields go into their filter menus as real
# categories. Re-check after any rebuild by looking for a label with nothing after it:
#   awk -F'\t' '$4==""||$15==""||$17==""||$18==""||$35==""' mavemdVar.bed | wc -l   # must be 0

# SCORE DIRECTION. Functional scores have no common scale and no common direction across score
# sets; the paper calls out MSH2, where LOWER scores mean normal activity. Nothing in the build
# compares raw scores across score sets, and the colors come from the functional class and the
# ACMG code rather than from the score, so this is a documentation problem only. It is called
# out in the Description section of mavemdVar.html.

# COLORS. Cells and items are colored by the ACMG evidence code on an 11-class red-to-blue
# diverging ramp, red toward pathogenic and blue toward benign, darkening with strength. The
# strength ladder is MaveDB's own StrengthOfEvidenceProvided enum, taken from their OpenAPI
# spec, not assumed: VERY_STRONG, STRONG, MODERATE_PLUS, MODERATE, SUPPORTING. Score sets that
# classify variants without assigning an ACMG code use a separate purple-green palette, so a
# measurement is not mistaken for a clinical claim.

# WHICH CALIBRATION IS SHOWN. A score set can carry up to six calibrations. MaveDB curates a
# `primary` flag with its own promote and demote API endpoints, so where that primary actually
# classifies the variant it is their editorial choice and is used as-is: 171,467 items.
#
# Two traps. 21 score sets have no primary at all (one of them has no calibration whatsoever).
# A further 35 designate one that carries no
# functionalClassifications, so it says nothing about any variant and arrives as NA on every
# CSV row - it exists but cannot be displayed. Those are different situations and the item's
# "Calibration chosen by" field distinguishes them, because telling a user "no primary
# designated" about a score set that has one is simply false. 172,287 items are in that second
# case; an earlier build mislabelled all of them.
#
# Every calibration's call is listed on the details page regardless of which one is displayed.

# Build results:
#   mavemdVar.bb   462,486 items, 47 fields (bed12+35)
#
# PROVENANCE OF THE ANNOTATIONS. clinvarRelease and gnomadVersion ride on each item that has a
# ClinVar or gnomAD annotation (60,666 and 3,978 respectively), so a stale annotation is visible
# rather than silent. This build is mixed: 79 score sets served ClinVar 2026-01 and 4 served
# 2025-01 (57,182 and 3,555 items respectively), on gnomAD v4.1 throughout. A single release
# per build is not something to expect;
# fetchMaveMd.py always asks for the newest ClinVar namespace the score set offers, so the
# release moves on its own between builds and is not pinned anywhere.

# LINKOUTS. The urls setting linkifies publicationUrl, scoreSet (to its MaveDB page) and
# clinGenId (to the ClinGen Allele Registry, which serves both CA genomic and PA protein
# alleles). hgc percent-encodes the substituted value, so the score set URN arrives as
# urn%3Amavedb%3A...%2D0%2D1; mavedb.org resolves that fine. There is no separate mavedbUrl
# field: it used to store 462,486 copies of 84 distinct strings that are just a prefix plus
# scoreSet, so linkifying scoreSet replaces it.
#
# DEFAULT DISPLAY. mavemdMap pack at priority 1, mavemdVar dense at priority 2. Pack height
# scales with how many items the window holds: over TP53 it measures 3,600 px at 30 bp and
# 7,090 px at 100 bp. Past 200 kb of coverage maxWindowCoverage suppresses the track in pack
# altogether, so pack is not something to open onto at gene scale; dense gives a compact
# colored strip at any zoom and the reader switches to pack at a few codons.
#   mavemdMap.bb        84 items, 32 fields (bed12+20)
#
# Install: the bigDataUrl and dataVersion settings point at /gbdb, so the links have to be
# repointed at each new dated build directory. runBuild.sh does not do this.
mkdir -p /gbdb/hg38/mavemd
cd /gbdb/hg38/mavemd
ln -sf /hive/data/genomes/hg38/bed/mavemd/2026-09-17/mavemdVar.bb mavemdVar.bb
ln -sf /hive/data/genomes/hg38/bed/mavemd/2026-09-17/mavemdMap.bb mavemdMap.bb
ln -sf /hive/data/genomes/hg38/bed/mavemd/2026-09-17/version.txt  version.txt

# Load trackDb. "make DBS=hg38" writes the per-user sandbox tables (trackDb_$USER); the track
# is gated "include mavemd.ra alpha" so "make alpha DBS=hg38" is what puts it on dev.
cd ~/kent/src/hg/makeDb/trackDb
make DBS=hg38
#
# STRAND is the transcript's, so it is set on every item placed on a protein: 190,188 + and
# 242,566 -. The remaining 29,732 are '.', the items placed at genomic level only, where
# there is no transcript to take a strand from: 22,411 from a MaveDB genomic term whose row
# has no usable single-substitution protein term, and 7,321 from the experiment's own
# reference term.

# Version file for the dataVersion trackDb setting.
cd /hive/data/genomes/hg38/bed/mavemd/2026-09-17
python3 -c "
import json
c=json.load(open('/hive/data/outside/mavemd/2026-09-17/collection.json'))
api=open('/hive/data/outside/mavemd/2026-09-17/apiVersion.txt').read().split()[-1]
open('version.txt','w').write('MaveDB collection modified %s (MaveDB API %s)\n'%(c['modificationDate'],api))
"

# trackDb: human/hg38/mavemd.ra, included from trackDb.ra as "include mavemd.ra alpha".
# The filterValues blocks in that file are generated by makeMaveMdVariants.py --raFragment;
# paste mavemdFilters.ra in rather than maintaining the value lists by hand. Paste it INSIDE
# the stanza with no blank lines: a blank line ends a .ra stanza, so blank separators between
# the filter groups silently orphan every setting after the first one. hgTrackDb loads without
# complaint and the filters just never appear in hgTrackUi. Check after a make with:
#   hgsql hg38 -Ne "select settings from trackDb_$USER where tableName='mavemdVar'" \
#     | tr '\\' '\n' | grep -c filterValues
#
# ClinVar terms containing a comma are escaped by doubling it, which is what
# slNameListFromCommaEscaped wants.
# No filterType is set on those fields on purpose: the list types run COMPARE_HASH_LIST_OR in
# hgTracks/bigBedTrack.c, which splits the field VALUE on commas before looking it up, so a
# ClinVar term with a comma could never match its own menu entry. The default whole-field
# match is correct here and still multi-select.
# relatedTracks.ra carries reciprocal links between mavedb and mavemd, plus one-way links
# from mavemd to clinvar and gnomadVariants.
