7485753cb34c1c3fd4a3f7b8894d993110db50ff lrnassar Thu Jul 23 14:03:39 2026 -0700 Rebuild popEVE heatmap with the full per-amino-acid matrix. refs #37791 The track was built from the genomic single-nucleotide VCF, so each column showed only the ~6 of 19 substitutions reachable by a single nucleotide change. The popEVE authors confirmed the per-transcript CSV download carries the complete matrix, so the converter now takes the all-19 scores from those CSVs while reusing the genomic codon coordinates, strand, and blocks derived from the VCF. Proteins absent from the (older) CSV release fall back to the sparse VCF scores. Color anchors recomputed over the full-matrix distribution, and the mouseover's raw EVE field is relabeled "EVE index" to distinguish it from the EVE track's 0-1 score. Updates the description page and makedoc to match. diff --git src/hg/makeDb/doc/hg38/popEve.txt src/hg/makeDb/doc/hg38/popEve.txt index a4d743dfb7b..94c34e3c581 100644 --- src/hg/makeDb/doc/hg38/popEve.txt +++ src/hg/makeDb/doc/hg38/popEve.txt @@ -43,38 +43,76 @@ # 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} # Notes: -# - popEVE is distributed as genomic SNVs, so only single-nucleotide-reachable missense -# substitutions are scored; each heatmap column (codon) therefore has ~6-9 of 19 rows -# filled. This is expected and sparser than the EVE track (which scores all 19). +# - 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 +# 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