e9a2b5a28cb4e3970977c76d84649bdd339a8423 lrnassar Thu Jul 30 17:00:51 2026 -0700 Address popEVE code-review feedback (v502). refs #37950 refs #37791 - popEve.ra dataVersion now names the score source (March per-transcript release) and the July VCF as the coordinate/strand source, rather than only the VCF date. - vcfToPopEveHeatmap.py: when the CSV wildtype disagrees with the genomic wildtype, skip the position and keep the correct sparse data instead of attaching CSV scores computed for a different residue; skip CSV rows with a nan/empty popEVE; add a posSparse counter for positions with genomic coordinates but no CSV row. All three are 0 on the current data, so the bigBed output is unchanged (verified byte-identical), but they make the converter fail safe for future per-transcript releases. - Add the build drivers runBuild.sh and runBuildDense.sh to the tree (the anchor computation previously lived only in the hive build directory), and add a makedoc forward-pointer so the intermediate sparse section is not mistaken for the final dense build. diff --git src/hg/makeDb/doc/hg38/popEve.txt src/hg/makeDb/doc/hg38/popEve.txt index 94c34e3c581..8fdcd8b45aa 100644 --- src/hg/makeDb/doc/hg38/popEve.txt +++ src/hg/makeDb/doc/hg38/popEve.txt @@ -1,118 +1,128 @@ # [Claude/lrnassar] popEVE - proteome-wide missense deleteriousness scores (2026-06-30) # popEVE places missense variants on a single proteome-wide, human-specific spectrum of # deleteriousness by calibrating EVE + ESM-1v evolutionary scores against UK Biobank # population variation with a Gaussian process. # Reference: Orenbuch et al. (2025) Nat Genet 57:3165-3174, PMID 41286104. # https://doi.org/10.1038/s41588-025-02400-1 # Data: the complete UKBB-trained GRCh38 VCF (popEVE v1.1, dated 2025-07-15) from the # pop.evemodel.org bulk downloads. The VCF carries, per missense substitution, genomic # coordinates plus INFO fields protein (RefSeq NP_), gene, mutant (e.g. E773D), # gap_frequency, popEVE, pop-adjusted_EVE, pop-adjusted_ESM1v, EVE, ESM1v. # Source size 1476124568 bytes, Last-Modified 2025-10-29. mkdir -p /hive/data/outside/popEve mkdir -p /hive/data/genomes/hg38/bed/popEve/sorttmp cd /hive/data/genomes/hg38/bed/popEve ln -s /hive/data/outside/popEve input cd /hive/data/outside/popEve wget -nv https://data.evemodel.org/popeve/v1.1/downloads/grch38_popEVE_ukbb_20250715.vcf.gz +# NOTE: the steps in this section produce the intermediate SPARSE heatmap. The final track is +# DENSE (all 19 substitutions per position); see the DENSE REBUILD section below. These steps +# still run because they produce popEve_sorted.tsv, the codon-coordinate and strand source +# reused by the dense build. +# # The combined VCF is one position-sorted file with all proteins interleaved (overlapping # genes share positions), so records are grouped by protein on disk before conversion. # Each protein becomes one heatmap BED12+ entry: columns = amino acid positions at codon # genomic coordinates, rows = 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 heatmap track. Multiple codon changes encoding # the same amino acid substitution carry identical popEVE scores and are deduplicated. # Wildtype cells are empty. popEVE lists only positions carrying a missense alt, so a codon # may have only 2 of its 3 genomic positions; block sizes are clamped (min(3, gap-to-next)) # so adjacent codon blocks cannot overlap, which keeps the file valid for bedToBigBed. # Strand is taken from NCBI RefSeq (the VCF has no strand field), cross-checked against the # strand inferred from genomic-vs-protein-position direction. # Color: global cross-gene gradient keyed to raw popEVE; interior anchors at the published # severe (-5.056) and moderate (-4.617) cutoffs and the proteome median (~-3.5); outer # saturation anchors at the 0.5th / 99.5th percentiles of the proteome-wide distribution. # Records with popEVE=nan (e.g. start-codon M1 variants) are skipped. cd /hive/data/genomes/hg38/bed/popEve # Full build (download already done): strand map, extract, color anchors, sort, convert, # filter to chrom.sizes, bigBed. See runBuild.sh in this directory. bash runBuild.sh # runBuild.sh steps, in order: # 1. NP_->strand map: # hgsql hg38 -N -e "select distinct l.protAcc, g.chrom, g.strand from ncbiRefSeqLink l # join ncbiRefSeq g on g.name=l.mrnaAcc where l.protAcc like 'NP\_%'" # -> np_strand.tsv (prefer the strand on a primary chromosome) # 2. extract: zcat ...vcf.gz | extractPopEve.py > popEve_records.tsv # 3. color anchors: p0.5 / p99.5 of popEVE over all records -> anchors.txt # 4. sort -t$'\t' -k1,1 -k4,4n -S 4G -T sorttmp popEve_records.tsv > popEve_sorted.tsv # 5. vcfToPopEveHeatmap.py popEve_sorted.tsv np_strand.tsv popEve_raw.bed $LO $HI # 6. bedSort + awk filter to chrom.sizes -> popEve_filtered.bed # 7. bedToBigBed -type=bed12+ -tab -as=popEve_heatmap.as popEve_filtered.bed chrom.sizes popEve.bb -# Scripts: ~/kent/src/hg/makeDb/scripts/popEve/{extractPopEve.py,vcfToPopEveHeatmap.py,popEve_heatmap.as} +# Scripts (all in ~/kent/src/hg/makeDb/scripts/popEve/): extractPopEve.py, +# vcfToPopEveHeatmap.py, popEve_heatmap.as, and the build drivers runBuild.sh (sparse / +# coordinate foundation) and runBuildDense.sh (dense rebuild). # Notes: # - The genomic VCF contains only single-nucleotide-reachable substitutions (~6-9 of 19 per # codon). The DENSE REBUILD below supersedes this with the full per-amino-acid matrix; the # VCF-derived popEve_sorted.tsv remains the source of genomic codon coordinates and strand. # - The heatmap renderer parses the score array (chopCommas, keeps trailing empty) and the # label array (chopByCharRespectDoubleQuotesKeepEmpty, drops one trailing empty) # differently. When the last cell (last row, last column) is empty the counts disagree and # the track aborts; the converter sets a non-empty placeholder label on that one trailing # cell (its score stays empty so the cell is uncolored). trailingFix count below. # Build results: # records extracted: 66,400,085 (nan skipped: 15,156) # color anchors: loAnchor=-5.742 hiAnchor=-2.287 (popEVE p0.5 / p99.5; median -3.358) # proteins: 18,968 (18,343 distinct gene symbols; remainder are RefSeq isoforms) # amino-acid positions: 10,114,809; bases covered: 886,638,890 # strand: 0 inferred-vs-RefSeq mismatches; 25 proteins not in the RefSeq map (strand from # coordinate inference); 0 with no strand signal; 0 dropped by chrom.sizes filter # trailingFix (trailing empty cell relabeled): 15,555 proteins # mouseover: HTML multi-line labels (<br>/<b>); component scores rounded to 3 dp; # missing/nan components shown as NA # popEve.bb size: 1,592,586,266 bytes # ============================================================================ # DENSE REBUILD (2026-07-23) - full per-amino-acid matrix, per collaborator request # ============================================================================ # The popEVE authors confirmed the per-transcript CSV download carries the COMPLETE # per-amino-acid matrix (all 19 substitutions per position), unlike the genomic VCF above # (single-nucleotide-reachable only, ~6 of 19 per codon). This rebuild takes the SCORES from # the per-transcript CSVs and REUSES the genomic codon coordinates, strand, and blocks # already derived from the VCF (popEve_sorted.tsv), producing dense heatmaps. Proteins present # in the VCF but absent from the (older, March 2025) CSV release fall back to sparse VCF scores. cd /hive/data/outside/popEve wget -nv https://data.evemodel.org/popeve/v1.1/downloads/popEVE_ukbb_20250312.zip # 2915714283 bytes, Last-Modified 2025-10-29 mkdir -p popeve_csv unzip -q popEVE_ukbb_20250312.zip -d popeve_csv/ # 18,332 per-transcript CSVs at popeve_csv/popEVE_ukbb_20250312/<NP_acc>.csv # columns: mutant,gap frequency,popEVE,popped EVE,popped ESM-1v,EVE,ESM-1v cd /hive/data/genomes/hg38/bed/popEve bash runBuildDense.sh # runBuildDense.sh steps: # 1. proteins not in the CSV (sparse fallback): 637 # 2. color anchors: p0.5 / p99.5 of popEVE over the full matrix (all CSV popEVE values plus # VCF popEVE for the 637 non-CSV proteins) -> anchorsDense.txt # 3. vcfToPopEveHeatmap.py popEve_sorted.tsv np_strand.tsv popEve_dense_raw.bed $LO $HI <csvDir> # the optional <csvDir> arg enables dense mode; loadCsv() fills all 19 substitutions # per position from the CSV, keeping the VCF codon coordinates and wildtype residue # 4. bedSort + filter to chrom.sizes + bedToBigBed -> popEve.bb # Dense build results: # full-matrix popEVE values: 187,547,662 # color anchors: loAnchor=-6.035 hiAnchor=-2.407 (popEVE p0.5 / p99.5; median -3.489) # proteins: 18,968 (dense 18,331 from CSV; sparse 637 VCF fallback) # NP_065910.3 (SHROOM3) is in the CSV but not the VCF coordinate set (all its VCF records # were nan-skipped start codons), so its dense matrix is unused -> 18,331 dense, not 18,332 -# wtMismatch 0 (CSV and VCF wildtype residues agree); trailingFix 1043 +# wtMismatch 0 (CSV and VCF wildtype residues agree); posSparse 0 (every dense position +# also had CSV coverage); trailingFix 1043. The converter skips CSV rows with a nan/empty +# popEVE and skips positions where the CSV wildtype disagrees with the genomic one; both +# are 0 in this build (empty-popEVE rows already render as empty cells). # popEve.bb size: 4,726,263,622 bytes mkdir -p /gbdb/hg38/popEve ln -sf /hive/data/genomes/hg38/bed/popEve/popEve.bb /gbdb/hg38/popEve/popEve.bb