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/scripts/popEve/vcfToPopEveHeatmap.py src/hg/makeDb/scripts/popEve/vcfToPopEveHeatmap.py index 16f06cdacba..1bc2aab01c4 100644 --- src/hg/makeDb/scripts/popEve/vcfToPopEveHeatmap.py +++ src/hg/makeDb/scripts/popEve/vcfToPopEveHeatmap.py @@ -95,30 +95,36 @@ with open(path) as fh: fh.readline() # header for line in fh: f = line.rstrip('\n').split(',') if len(f) < 7: continue m = f[0] # e.g. G1042A wt, var = m[0], m[-1] if wt not in STANDARD_SET or var not in STANDARD_SET: continue try: pos = int(m[1:-1]) except ValueError: continue gap, pe, paEve, paEsm, eve, esm = f[1], f[2], f[3], f[4], f[5], f[6] + try: # skip rows with no usable popEVE score + fpe = float(pe) + except ValueError: + continue + if fpe != fpe or fpe in (float('inf'), float('-inf')): # NaN / inf + continue d = data.get(pos) if d is None: d = data[pos] = {'wt': wt, 'vars': {}} d['vars'][var] = (pe, eve, esm, paEve, paEsm, gap) return data 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: @@ -153,33 +159,39 @@ 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)) # Densify: where the full per-amino-acid CSV covers a position, replace the sparse # single-nucleotide-reachable scores with the complete set of 19 substitutions. Codon # coordinates (bases) and the wildtype residue stay from the genomic records. if csvData is not None: stats['dense'] += 1 for protPos, d in byPos.items(): cd = csvData.get(protPos) if cd is None: + # Position has genomic coordinates but no CSV row: leave it sparse. + stats['posSparse'] += 1 continue if cd['wt'] != d['wt']: + # CSV keyed to a different residue at this position (e.g. a different + # transcript version): keep the correct sparse genomic data rather than + # attaching CSV scores computed for the wrong wildtype. stats['wtMismatch'] += 1 + continue d['scores'] = {var: v[0] for var, v in cd['vars'].items()} d['extra'] = {var: v[1:] for var, v in cd['vars'].items()} else: stats['sparse'] += 1 # 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) @@ -289,48 +301,50 @@ sys.exit("Usage: %s <sorted_tsv> <strandMap> <output_bed> <loAnchor> <hiAnchor> [csvDir]" % sys.argv[0]) sortedTsv = sys.argv[1] strandMapPath = sys.argv[2] outputBed = sys.argv[3] loAnchor = float(sys.argv[4]) hiAnchor = float(sys.argv[5]) csvDir = sys.argv[6] if len(sys.argv) > 6 else None strandMap = loadStrandMap(strandMapPath) sys.stderr.write("Loaded %d NP_->strand mappings.\n" % len(strandMap)) if csvDir: sys.stderr.write("Dense mode: full per-amino-acid matrices from %s\n" % csvDir) stats = {'ok': 0, 'multiChrom': 0, 'overlap': 0, 'strandMismatch': 0, - 'strandDefault': 0, 'trailingFix': 0, 'dense': 0, 'sparse': 0, 'wtMismatch': 0} + 'strandDefault': 0, 'trailingFix': 0, 'dense': 0, 'sparse': 0, 'wtMismatch': 0, + 'posSparse': 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, loadCsv(csvDir, curProt), 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, loadCsv(csvDir, curProt), strandMap, loAnchor, hiAnchor, out, stats) sys.stderr.write("Done: %(ok)d proteins written; dense %(dense)d, sparse %(sparse)d, " "multiChrom %(multiChrom)d, overlap %(overlap)d, " "strandMismatch %(strandMismatch)d, strandDefault %(strandDefault)d, " - "wtMismatch %(wtMismatch)d, trailingFix %(trailingFix)d\n" % stats) + "wtMismatch %(wtMismatch)d, posSparse %(posSparse)d, " + "trailingFix %(trailingFix)d\n" % stats) if __name__ == '__main__': main()