15bdf707a4bde0a281109bb810499a2fa98ef6ac
lrnassar
Tue Jul 14 07:48:05 2026 -0700
Fix stale "row Y" reference in popEVE trailingFix comment per CR feedback. refs #37791
The amino acid row reorder made P the last row, not Y, but the trailingFix explanation in
the converter comment and the makedoc still said "row Y, last column". Reworded both to be
order-agnostic ("last row, last column"). Documentation only; the code already operates
positionally on labelParts[-1].
diff --git src/hg/makeDb/scripts/popEve/vcfToPopEveHeatmap.py src/hg/makeDb/scripts/popEve/vcfToPopEveHeatmap.py
index a40682bd110..ceb315681c7 100644
--- src/hg/makeDb/scripts/popEve/vcfToPopEveHeatmap.py
+++ src/hg/makeDb/scripts/popEve/vcfToPopEveHeatmap.py
@@ -1,268 +1,268 @@
#!/usr/bin/env python3
"""Convert popEVE records (extracted, sorted by protein then position) to heatmap bigBed.
Usage: vcfToPopEveHeatmap.py <sorted_tsv> <strandMap> <output_bed> [loAnchor] [hiAnchor]
Input <sorted_tsv> is the output of extractPopEve.py sorted by:
sort -t$'\\t' -k1,1 -k4,4n
so every protein's records are contiguous and in ascending genomic order.
<strandMap> is a TSV "protAcc<TAB>strand" (one preferred mapping per NP_ accession,
from hg38 ncbiRefSeq) used to validate / override the strand inferred from coordinates.
Produces one heatmap BED12+ line per protein (see popEve_heatmap.as):
Columns = protein positions ordered by ascending genomic coordinate; reverse-strand
genes therefore read C-terminus -> N-terminus, matching genomic orientation.
Rows = 20 standard amino acids (A C D E F G H I K L M N P Q R S T V W Y).
Each codon is a block. popEVE lists only positions carrying a missense alt, so a codon
may have 2 of 3 genomic positions; block sizes are therefore clamped so adjacent blocks
cannot overlap (min(3, gap-to-next)), keeping the file valid for bedToBigBed.
Wildtype cells are left empty. Mouseover carries popEVE + severity class plus the
component EVE / ESM-1v / pop-adjusted scores and gap frequency.
"""
import sys
STANDARD_AAS = list('AVLIMFYWRHKDESTNQGCP') # 20 standard AAs, by class, to match MaveDB
# Published popEVE severity cutoffs (fixed interior color anchors).
SEVERE_MAX = -5.056 # popEVE < SEVERE_MAX -> severe
MODERATE_MAX = -4.617 # SEVERE_MAX <= popEVE < this -> moderate; >= this -> tolerated
MEDIAN_ANCHOR = -3.5 # interior white anchor (~ proteome-wide median)
COLOR_VALUES = "#b2182b,#d6604d,#f4a582,#f7f7f7,#2166ac" # deleterious red -> tolerated blue
def r3(s):
"""Round a numeric string to 3 decimals, trimming trailing zeros."""
try:
v = float(s)
except (ValueError, TypeError):
return s
out = ('%.3f' % v).rstrip('0').rstrip('.')
return '0' if out in ('', '-0') else out
def comp(s):
"""Format a component score for the mouseover: 'NA' for missing/nan, else 3 dp."""
if s == '' or s == 'nan':
return 'NA'
return r3(s)
def classify(score):
"""Severity class for a popEVE score."""
if score < SEVERE_MAX:
return "severe"
if score < MODERATE_MAX:
return "moderate"
return "tolerated"
def loadStrandMap(path):
"""Read protAcc<TAB>strand into a dict."""
d = {}
with open(path) as fh:
for line in fh:
f = line.rstrip('\n').split('\t')
if len(f) >= 2 and f[1] in ('+', '-'):
d[f[0]] = f[1]
return d
def inferStrand(sortedPos):
"""Infer strand from whether protein position decreases as genomic coordinate increases.
sortedPos is a list of (genomic_start0, prot_pos) sorted by genomic_start0.
Returns '+', '-', or None (no signal, e.g. single codon)."""
first = None
last = None
for _, pp in sortedPos:
if first is None:
first = pp
last = pp
if first is None or first == last:
return None
return '-' if last < first else '+'
def buildEntry(protein, gene, records, strandMap, loAnchor, hiAnchor, out, stats):
"""Write one heatmap BED12+ line for a protein.
records: list of (chrom, pos0, wt, protPos, var, popEVE, EVE, ESM1v, popAdjEVE,
popAdjESM1v, gapFreq) where pos0 is the 0-based genomic base.
"""
# Group by protein position: wt aa, codon's genomic base set, and {var: (score, extras)}.
byPos = {}
chroms = set()
for (chrom, pos0, wt, protPos, var, pe, eve, esm, paEve, paEsm, gap) in records:
chroms.add(chrom)
d = byPos.get(protPos)
if d is None:
d = byPos[protPos] = {'wt': wt, 'bases': set(), 'scores': {}, 'extra': {}}
d['bases'].add(pos0)
if var not in d['scores']: # dedup: identical substitution carries identical score
d['scores'][var] = pe
d['extra'][var] = (eve, esm, paEve, paEsm, gap)
if len(chroms) > 1:
sys.stderr.write(" WARNING: %s (%s) spans multiple chromosomes %s - skipping\n" %
(gene, protein, sorted(chroms)))
stats['multiChrom'] += 1
return
chrom = next(iter(chroms))
# One column per protein position, ordered by ascending genomic coordinate.
# Column start = min(codon bases) (0-based); size clamped below.
cols = sorted(((min(d['bases']), protPos) for protPos, d in byPos.items()))
colStarts0 = [c[0] for c in cols]
colProtPos = [c[1] for c in cols]
nCols = len(cols)
# Strand: infer from coordinates, then validate / override against ncbiRefSeq.
inferred = inferStrand(cols)
refStrand = strandMap.get(protein)
if refStrand is not None:
if inferred is not None and inferred != refStrand:
sys.stderr.write(" WARNING: %s (%s) strand %s != ncbiRefSeq %s; using RefSeq\n"
% (gene, protein, inferred, refStrand))
stats['strandMismatch'] += 1
strand = refStrand
elif inferred is not None:
strand = inferred
else:
strand = '+'
stats['strandDefault'] += 1
sys.stderr.write(" NOTE: %s (%s) no strand signal and not in ncbiRefSeq - defaulting +\n" %
(gene, protein))
chromStart = colStarts0[0]
# Clamp block sizes so adjacent blocks cannot overlap (popEVE codons may have only
# 2 of 3 positions, so naive size-3 blocks would collide -> bedToBigBed aborts).
blockSizes = []
for i in range(nCols):
if i < nCols - 1:
gap = colStarts0[i + 1] - colStarts0[i]
blockSizes.append(max(1, min(3, gap)))
else:
blockSizes.append(3) # last block keeps full codon width
relStarts = [s - chromStart for s in colStarts0]
chromEnd = colStarts0[-1] + blockSizes[-1]
# Build assertion: blocks strictly non-overlapping, last ends at chromEnd.
for i in range(1, nCols):
if relStarts[i] < relStarts[i - 1] + blockSizes[i - 1]:
sys.stderr.write(" ERROR: %s (%s) overlapping blocks at col %d\n" % (gene, protein, i))
stats['overlap'] += 1
return
if chromStart + relStarts[-1] + blockSizes[-1] != chromEnd:
sys.stderr.write(" ERROR: %s (%s) last block does not end at chromEnd\n" % (gene, protein))
return
# BED score (0-1000): most-deleterious residue -> high, spread across loAnchor..hiAnchor.
minPe = min(float(d2) for protPos in byPos for d2 in byPos[protPos]['scores'].values())
span = hiAnchor - loAnchor
bedScore = int(round((hiAnchor - minPe) / span * 1000)) if span else 0
bedScore = max(0, min(1000, bedScore))
rowCount = len(STANDARD_AAS)
labels = ','.join(STANDARD_AAS)
colorBounds = "%s,%s,%s,%s,%s" % (loAnchor, SEVERE_MAX, MODERATE_MAX, MEDIAN_ANCHOR, hiAnchor)
# Row-major score and label arrays (all columns for row 0, then row 1, ...).
scoreParts = []
labelParts = []
for aa in STANDARD_AAS:
for protPos in colProtPos:
d = byPos[protPos]
wt = d['wt']
if aa == wt:
scoreParts.append('') # wildtype cell -> empty (background)
labelParts.append('')
elif aa in d['scores']:
pe = d['scores'][aa]
eve, esm, paEve, paEsm, gap = d['extra'][aa]
scoreParts.append(r3(pe))
cls = classify(float(pe))
# Mouseover rendered as HTML (<br>/<b>). No commas inside labels (the field is
# comma-split by the renderer); values use '.' and '/' only, so labels are unquoted.
lbl = ("%s%s→%s<br><b>popEVE:</b> %s (%s)<br><b>EVE:</b> %s<br>"
"<b>ESM1v:</b> %s<br><b>popAdj:</b> EVE %s / ESM1v %s<br><b>gap:</b> %s" %
(wt, protPos, aa, r3(pe), cls, comp(eve), comp(esm), comp(paEve),
comp(paEsm), comp(gap)))
labelParts.append(lbl)
else:
scoreParts.append('')
labelParts.append('')
# The heatmap renderer parses the score array with chopCommas (keeps a trailing empty
# field) but the label array with chopByCharRespectDoubleQuotesKeepEmpty (drops one
- # trailing empty field). When the very last cell (row Y, last column) is empty - common on
- # minus-strand genes whose last column is the start codon, whose M1 substitutions are
- # 'nan' and were skipped - the two field counts disagree and the track aborts. Guarantee a
+ # trailing empty field). When the very last cell (last row, last column) has no scored
+ # substitution - common, since each codon has only a few single-nucleotide-reachable
+ # substitutions - the two field counts disagree and the track aborts. Guarantee a
# non-empty final label cell; the empty score keeps that cell uncolored (background).
if labelParts[-1] == '':
labelParts[-1] = '(no popEVE score)'
stats['trailingFix'] += 1
fields = [
chrom, chromStart, chromEnd,
gene, bedScore, strand,
chromStart, chromEnd, 0,
nCols,
','.join(str(s) for s in blockSizes) + ',',
','.join(str(s) for s in relStarts) + ',',
rowCount, labels,
colorBounds, COLOR_VALUES,
','.join(scoreParts), ','.join(labelParts),
"popEVE: red <= -5.056 severe; orange moderate; white ~ -3.5; blue tolerated",
protein,
]
out.write('\t'.join(str(f) for f in fields) + '\n')
stats['ok'] += 1
def main():
if len(sys.argv) < 4:
sys.exit("Usage: %s <sorted_tsv> <strandMap> <output_bed> [loAnchor] [hiAnchor]"
% sys.argv[0])
sortedTsv = sys.argv[1]
strandMapPath = sys.argv[2]
outputBed = sys.argv[3]
loAnchor = float(sys.argv[4]) if len(sys.argv) > 4 else -7.0
hiAnchor = float(sys.argv[5]) if len(sys.argv) > 5 else -2.0
strandMap = loadStrandMap(strandMapPath)
sys.stderr.write("Loaded %d NP_->strand mappings.\n" % len(strandMap))
stats = {'ok': 0, 'multiChrom': 0, 'overlap': 0, 'strandMismatch': 0,
'strandDefault': 0, 'trailingFix': 0}
with open(sortedTsv) as fh, open(outputBed, 'w') as out:
curProt = None
curGene = None
recs = []
for line in fh:
f = line.rstrip('\n').split('\t')
if len(f) < 13:
continue
(protein, gene, chrom, pos, wt, protPos, var, pe,
eve, esm, paEve, paEsm, gap) = f[:13]
if protein != curProt:
if curProt is not None:
buildEntry(curProt, curGene, recs, strandMap, loAnchor, hiAnchor, out, stats)
curProt = protein
curGene = gene
recs = []
recs.append((chrom, int(pos) - 1, wt, int(protPos), var, pe,
eve, esm, paEve, paEsm, gap))
if curProt is not None:
buildEntry(curProt, curGene, recs, strandMap, loAnchor, hiAnchor, out, stats)
sys.stderr.write("Done: %(ok)d proteins written; multiChrom %(multiChrom)d, "
"overlap %(overlap)d, strandMismatch %(strandMismatch)d, "
"strandDefault %(strandDefault)d, trailingFix %(trailingFix)d\n" % stats)
if __name__ == '__main__':
main()