a37d4ffb473b6ebf6c66545f2b7d5f7eef35ef4b max Fri Apr 10 03:00:33 2026 -0700 Color strVar subtracks by expected heterozygosity instead of motif period, fix hgTrackUi filter label truncation, refs #36652 Change all four strVar subtracks (webstr, tommoStr, trexplorer, viennaVntr) from motif-period-based coloring to expected heterozygosity (het = 1 - sum(p_i^2)), using a blue-to-red heat map: dark blue (het<0.1) through medium blue, light purple, salmon, to dark red (het>=0.7). Add het as a filterable bigBed field with scoreFilter and filterByRange on each track. Update mouseOver, track docs, and makedoc. Also fix hgTrackUi to strip the "|..." suffix from autoSql comments when displaying numeric filter labels. Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com> diff --git src/hg/makeDb/scripts/webstr/webstrToBed.py src/hg/makeDb/scripts/webstr/webstrToBed.py index 7754047adb4..f26bab2b5d8 100644 --- src/hg/makeDb/scripts/webstr/webstrToBed.py +++ src/hg/makeDb/scripts/webstr/webstrToBed.py @@ -1,126 +1,160 @@ #!/usr/bin/env python3 """Convert WebSTR CSV data to BED9+ format for bigBed conversion. Reads hg38_repeats_withlinks.csv.gz and hg38_afreqs.csv.gz from the input directory and writes a tab-separated BED file to stdout. Usage: webstrToBed.py <inputDir> > webstr.bed """ import csv import gzip import sys from collections import defaultdict -PERIOD_COLORS = { - 1: "255,0,0", # mono: red - 2: "0,0,255", # di: blue - 3: "0,128,0", # tri: green - 4: "255,165,0", # tetra: orange - 5: "128,0,128", # penta: purple - 6: "70,130,180", # hexa: steel blue -} -DEFAULT_COLOR = "128,128,128" # gray for period > 6 +HET_BINS = [ + (0.1, "0,0,180"), # het < 0.1: dark blue (nearly monomorphic) + (0.3, "70,130,230"), # het 0.1-0.3: medium blue (low diversity) + (0.5, "180,130,200"), # het 0.3-0.5: light purple (moderate diversity) + (0.7, "230,100,80"), # het 0.5-0.7: salmon (high diversity) + (999, "180,0,0"), # het >= 0.7: dark red (very high diversity) +] +NO_DATA_COLOR = "128,128,128" # gray when no allele freq data def truncateMotif(motif, maxLen=25): """Truncate motif to maxLen characters with '..' in the middle.""" if len(motif) <= maxLen: return motif keepLen = maxLen - 2 leftLen = (keepLen + 1) // 2 rightLen = keepLen - leftLen return motif[:leftLen] + ".." + motif[-rightLen:] +def hetToColor(het): + """Map heterozygosity value to an RGB color string.""" + if het < 0: + return NO_DATA_COLOR + for threshold, color in HET_BINS: + if het < threshold: + return color + return HET_BINS[-1][1] + + +def computeHet(af): + """Compute expected heterozygosity from allele freqs pooled across populations. + + Pools allele frequencies weighted by sample count, then computes + het = 1 - sum(p_i^2). + """ + if af is None: + return -1.0 + totalN = 0 + weightedFreqs = defaultdict(float) + for cohort in COHORT_ORDER: + entry = af[cohort] + n = entry["n"] + if n == 0: + continue + totalN += n + for allele, freq in zip(entry["alleles"], entry["freqs"]): + weightedFreqs[allele] += float(freq) * n + if totalN == 0: + return -1.0 + sumPSq = sum((wf / totalN) ** 2 for wf in weightedFreqs.values()) + return round(1.0 - sumPSq, 3) + + COHORT_ORDER = ["AFR", "AMR", "EAS", "EUR", "SAS"] COHORT_MAP = { "1000 Genomes AFR": "AFR", "1000 Genomes AMR": "AMR", "1000 Genomes EAS": "EAS", "1000 Genomes EUR": "EUR", "1000 Genomes SAS": "SAS", } def loadAlleleFreqs(inDir): """Load allele frequency data, grouped by repeatid and cohort.""" freqs = defaultdict(lambda: {c: {"alleles": [], "freqs": [], "n": 0} for c in COHORT_ORDER}) path = f"{inDir}/hg38_afreqs.csv.gz" with gzip.open(path, "rt") as f: reader = csv.reader(f) header = next(reader) # skip header for row in reader: cohort_raw, allele, freq, n, repeatid = row cohort = COHORT_MAP.get(cohort_raw) if cohort is None: continue entry = freqs[repeatid][cohort] entry["alleles"].append(allele) entry["freqs"].append(freq) entry["n"] = int(n) return freqs def main(): if len(sys.argv) != 2: print(__doc__, file=sys.stderr) sys.exit(1) inDir = sys.argv[1] print("Loading allele frequencies...", file=sys.stderr) afreqs = loadAlleleFreqs(inDir) print(f" Loaded frequencies for {len(afreqs)} repeats", file=sys.stderr) print("Processing repeats...", file=sys.stderr) repeatsPath = f"{inDir}/hg38_repeats_withlinks.csv.gz" count = 0 with gzip.open(repeatsPath, "rt") as f: reader = csv.reader(f) header = next(reader) # skip header for row in reader: repeatid, panel, chrom, motif, start, end, period, numcopies, _webstr_link = row - period_int = int(period) - color = PERIOD_COLORS.get(period_int, DEFAULT_COLOR) - # Source coordinates are 1-based; convert start to 0-based for BED start = str(int(start) - 1) + # Compute heterozygosity pooled across populations + af = afreqs.get(repeatid) + het = computeHet(af) + color = hetToColor(het) + score = max(0, int(het * 1000)) if het >= 0 else 0 + # BED9 fields fields = [ chrom, start, end, truncateMotif(motif) + "x" + numcopies, # name - "0", # score + str(score), # score (het * 1000) ".", # strand start, # thickStart end, # thickEnd color, # itemRgb motif, period, numcopies, repeatid, + str(het), # het field ] - - # Allele frequency fields for each cohort (logfmt: allele=freq pairs) - af = afreqs.get(repeatid) for cohort in COHORT_ORDER: if af and af[cohort]["alleles"]: entry = af[cohort] pairs = [a + "=" + f for a, f in zip(entry["alleles"], entry["freqs"])] fields.append(" ".join(pairs)) fields.append(str(entry["n"])) else: fields.append("") fields.append("0") print("\t".join(fields)) count += 1 if count % 500000 == 0: print(f" Processed {count} repeats...", file=sys.stderr) print(f"Done. Wrote {count} records.", file=sys.stderr) if __name__ == "__main__": main()