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/makeMaveMdVariants.py src/hg/makeDb/scripts/mavemd/makeMaveMdVariants.py index 7fb55163d01..8a56a426f06 100755 --- src/hg/makeDb/scripts/mavemd/makeMaveMdVariants.py +++ src/hg/makeDb/scripts/mavemd/makeMaveMdVariants.py @@ -33,30 +33,35 @@ sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import mavemdLib as lib csv.field_size_limit(10 ** 7) # MaveDB states protein terms against RefSeq for most score sets and Ensembl for a few. PROTEIN_TERM = re.compile(r'^(?P<acc>(?:N[PM]_|ENSP)[\w.]+):p\.' r'(?P<wt>[A-Z][a-z]{2})(?P<pos>\d+)' r'(?P<var>[A-Z][a-z]{2}|Ter|=)$') NA = ('NA', '', '-', None) # Above this fraction of disagreements between our codon projection and MaveDB's own # genomic mapping, stop rather than publish coordinates we no longer trust. MAX_PROJECTION_MISMATCH = 0.005 +# Per-accession ceiling for the wild-type residue check. Individual variants that fail are +# dropped rather than placed at a position the genome contradicts, so this threshold only +# has to catch a transcript that has moved wholesale, which fails at close to 100%. +MAX_RESIDUE_MISMATCH = 0.25 + def clean(value): """Normalise the CSV's several spellings of "no value" to an empty string.""" return '' if value in NA else value.strip() def basesToBlocks(bases, chromStart): """Contiguous runs of a codon's genomic bases, as BED block starts and sizes.""" ordered = sorted(bases) starts, sizes = [ordered[0]], [1] for prev, cur in zip(ordered, ordered[1:]): if cur == prev + 1: sizes[-1] += 1 else: starts.append(cur) @@ -252,72 +257,78 @@ # the bases that actually change. if len(ref) > 1 and len(alt) >= 1 and ref[0] == alt[0]: start += 1 if end <= start: end = start + 1 placements[term] = (chrom, start, end) return placements, failures 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('--classPalette', default='purple', choices=sorted(lib.CLASS_PALETTES), help='palette for measurements with no ACMG code') parser.add_argument('--raFragment', help='write generated trackDb filterValues here') parser.add_argument('--workDir', help='scratch directory (default: alongside outBed)') args = parser.parse_args() lib.setClassPalette(args.classPalette) workDir = args.workDir or os.path.dirname(os.path.abspath(args.outBed)) or '.' os.makedirs(workDir, exist_ok=True) scoreSets = loadScoreSets(args.downloadDir) sys.stderr.write("Loaded %d score set records\n" % len(scoreSets)) scoreSetTx = loadScoreSetTranscripts(args.downloadDir) sys.stderr.write("Resolved a submission transcript for %d of %d score sets\n" % (len(scoreSetTx), len(scoreSets))) sys.stderr.write("Pass 1: collecting HGVS terms\n") terms, protAccs = collectTerms(args.downloadDir, scoreSetTx) sys.stderr.write(" %d distinct HGVS terms, %d protein accessions\n" % (len(terms), len(protAccs))) sys.stderr.write("Resolving protein accessions to transcripts\n") protToTx, unresolvedProt = lib.loadProteinToTranscript(args.db, sorted(protAccs)) if unresolvedProt: sys.stderr.write(" WARNING: no transcript for %s\n" % ', '.join(unresolvedProt)) codonMaps, missingTx = lib.loadCodonMaps(args.db, sorted(set(protToTx.values()))) + lib.addProteinSequence(args.db, codonMaps, args.twoBit, workDir) if missingTx: sys.stderr.write(" WARNING: no ncbiRefSeqCurated entry for %s\n" % ', '.join(missingTx)) for acc, tx in sorted(protToTx.items()): cm = codonMaps.get(tx) if cm: sys.stderr.write(" %-16s -> %-16s %s %s %d codons\n" % (acc, tx, cm.chrom, cm.strand, cm.protLen)) sys.stderr.write("Running hgvsToVcf on %d terms\n" % len(terms)) placements, hgvsFailures = runHgvsToVcf(args.db, terms, workDir) sys.stderr.write(" placed %d, failures %s\n" % (len(placements), dict(hgvsFailures))) stats = collections.Counter() filterValues = collections.defaultdict(set) + residueChecked = collections.Counter() + residueBad = collections.Counter() + residueExamples = [] projectionChecked = 0 projectionMismatch = 0 mismatchExamples = [] rows = [] sys.stderr.write("Pass 2: building items\n") 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'): clinvarRelease = name.split('.')[1] for row in reader: @@ -326,38 +337,61 @@ meta = scoreSets.get(urn, {}) gene = meta.get('gene') or clean(row.get('score_set.target_gene')) genomic = clean(row.get('mavedb.post_mapped_hgvs_g')) protein = clean(row.get('mavedb.post_mapped_hgvs_p')) original = qualifyTerm(clean(row.get('hgvs_nt')), scoreSetTx.get(urn)) protMatch = PROTEIN_TERM.match(protein) if protein else None # Codon projection, computed whenever there is a usable protein term so it # can be cross-checked even when a genomic term is also present. A codon # at an exon junction is split, so the item spans the intervening intron: # that is the real extent of the codon, and a protein-level measurement # pins down nothing narrower. projected = None codonBases = None + residueMismatch = False strand = '.' if protMatch: tx = protToTx.get(protMatch.group('acc')) codonMap = codonMaps.get(tx) if tx else None if codonMap: codonBases = codonMap.codon(int(protMatch.group('pos'))) strand = codonMap.strand - if codonBases: + # Does the genome actually hold the residue the HGVS term asserts? + # This covers every projected codon, unlike the comparison against + # MaveDB's own genomic mapping, which can only see the accessions + # that happen to have genomic-route variants too. + pos1 = int(protMatch.group('pos')) + wtOne = lib.THREE_TO_ONE.get(protMatch.group('wt')) + if codonMap.protein and wtOne and pos1 <= len(codonMap.protein): + acc = protMatch.group('acc') + residueChecked[acc] += 1 + if codonMap.protein[pos1 - 1] != wtOne: + # The codon we would place does not hold the residue the + # term names, so the projection would put this variant at + # the wrong position. Refuse to place it. In this build + # every case is a nonsense term whose numbering is one + # codon downstream of its own reference residue. + residueBad[acc] += 1 + residueMismatch = True + if len(residueExamples) < 10: + residueExamples.append( + '%s %s: term says %s, genome has %s' + % (urn, protein, wtOne, + codonMap.protein[pos1 - 1])) + if codonBases and not residueMismatch: start, end = codonMap.codonSpan(int(protMatch.group('pos'))) projected = (codonMap.chrom, start, end) place = None method = '' if genomic and genomic in placements: place = placements[genomic] method = 'MaveDB genomic mapping' if projected and codonBases: # Compare against the codon's actual three bases, not its span: # a split codon's span includes an intron that the genomic # coordinate should never fall in. projectionChecked += 1 overlap = (projected[0] == place[0] and any(place[1] <= b < place[2] for b in codonBases)) @@ -365,35 +399,43 @@ projectionMismatch += 1 if len(mismatchExamples) < 10: mismatchExamples.append( '%s %s g=%s:%d-%d codon bases=%s' % (urn, protein, place[0], place[1], place[2], ','.join(str(b + 1) for b in sorted(codonBases)))) elif projected: place = projected method = 'codon projection' elif original and original in placements: place = placements[original] method = 'hgvsToVcf on submitted term' if place is None: submittedProtein = clean(row.get('hgvs_pro')) - if submittedProtein.startswith('p.[') or ';' in submittedProtein: + submittedNt = clean(row.get('hgvs_nt')) + # A haplotype can be stated on either axis. PTEN 00000054-a-1 states + # 1,236 of them as nucleotides only, e.g. c.[1207G>T;1209C>T], with no + # protein term at all; testing hgvs_pro alone files them under + # "submitted term rejected" and understates the haplotype count. + if (submittedProtein.startswith('p.[') or ';' in submittedProtein + or submittedNt.startswith('c.[') or ';' in submittedNt): # A haplotype is several substitutions measured as one unit, so it # has no single position. The earlier MaveDB track excluded these # for the same reason. stats['skipHaplotype'] += 1 + elif residueMismatch: + stats['skipResidueMismatch'] += 1 elif protein and not protMatch: stats['skipProteinTermNotASubstitution'] += 1 elif genomic: stats['skipGenomicTermRejected'] += 1 elif original: stats['skipSubmittedTermRejected'] += 1 else: stats['skipNoUsableTerm'] += 1 continue chrom, start, end = place # Draw the real bases as blocks. For a codon split across an exon junction # that is two blocks with a thin connector, rather than one solid bar across # the intron: the intron is not part of the codon, and a solid bar there # shows a measurement where none was made and pollutes any range query. @@ -504,30 +546,51 @@ if funcClass: filterValues['funcClass'].add(funcClass) if gene: filterValues['gene'].add(gene) # Harvest from the named values, never from positions in `item`: the row # layout has changed several times and positional indices drift silently # into the wrong column. for field, value in (('clinvarSig', clinvarSig), ('vepConsequence', vepConsequence), ('assayMethod', assayMethod), ('assayModel', assayModel), ('libraryMethod', libraryMethod)): if value: filterValues[field].add(value) + if residueChecked: + sys.stderr.write("Wild-type residue check, per protein accession:\n") + failed = [] + for acc in sorted(residueChecked): + n, bad = residueChecked[acc], residueBad[acc] + rate = bad / n + flag = (' <-- FAILS' if rate > MAX_RESIDUE_MISMATCH + else (' (dropped)' if bad else '')) + sys.stderr.write(" %-18s %6d checked, %5d mismatched (%.3f%%)%s\n" + % (acc, n, bad, 100 * rate, flag)) + if rate > MAX_RESIDUE_MISMATCH: + failed.append(acc) + for example in residueExamples: + sys.stderr.write(" %s\n" % example) + if failed: + sys.exit("ERROR: the reference residue in the genome disagrees with the HGVS " + "term for more than %.0f%% of the projected codons in: %s. That is a " + "moved transcript or assembly annotation, not stray upstream rows; " + "do not publish these coordinates." + % (100 * MAX_RESIDUE_MISMATCH, ', '.join(failed))) + if projectionChecked: rate = projectionMismatch / projectionChecked sys.stderr.write("Codon projection cross-check: %d of %d disagreed with MaveDB's " "genomic mapping (%.4f%%)\n" % (projectionMismatch, projectionChecked, 100 * rate)) for example in mismatchExamples: sys.stderr.write(" %s\n" % example) if rate > MAX_PROJECTION_MISMATCH: sys.exit("ERROR: codon projection disagrees with MaveDB on %.2f%% of the " "variants that carry both coordinate systems, above the %.2f%% " "threshold. RefSeq or the MaveDB mapper has moved; investigate before " "publishing." % (100 * rate, 100 * MAX_PROJECTION_MISMATCH)) rows.sort(key=lambda r: (r[0], r[1], r[2])) with open(args.outBed, 'w') as out: @@ -552,21 +615,24 @@ values = sorted(filterValues[field]) if not values: continue # filterValues is read with slNameListFromCommaEscaped, which takes a doubled # comma as a literal one, so ClinVar terms like "Pathogenic, low penetrance" # are escaped to build the menu correctly. # # No filterType is emitted on purpose. The list types (multipleListOr and # friends) 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 in it could never match its own menu entry. Every one of these fields # holds a single value, so the default FILTERBY_MULTIPLE is both correct and # still multi-select. escaped = [v.replace(',', ',,') for v in values] fh.write('filterValues.%s %s\n' % (field, ','.join(escaped))) - fh.write('filterLabel.%s %s\n\n' % (field, label)) + # One newline, not two: a blank line ends a trackDb stanza, so a + # fragment with blank separators silently orphans every setting after + # the first when it is pasted in. + fh.write('filterLabel.%s %s\n' % (field, label)) sys.stderr.write("Wrote trackDb filter fragment to %s\n" % args.raFragment) if __name__ == '__main__': main()