2c0adaa48b2f9a14109c3f90713405f259d920bc max Wed Sep 9 06:08:46 2026 -0700 hg38: pcLAI local ancestry track for HPRC Release 2 haplotypes Point cloud local ancestry inference (pcLAI) for HPRC Release 2, projected onto GRCh38: a composite with one subtrack per haplotype, both haplotypes of 231 samples plus CHM13, 463 in all. 11,936,603 windows, autosomes only, no windows dropped from the source files. Each window carries the (PC1,PC2) coordinate pcLAI predicts for it, the discretized ancestry centroid, and a confidence score. The details page draws the window's position against the 1000 Genomes reference panel that defines the space, via detailsScript/scatterPlot; metaDataUrl is what lets hgTrackUi serve that panel file for a native (non-hub) track. thickStart is one base before chromStart in 54,811 of the windows (0.46%), against the format the pcLAI authors document, so bedToBigBed rejects it. Neither thick column carries information here, so both are set to the item bounds rather than dropping those windows. Testing only for now, alpha, no ticket yet; a ticket may follow if this becomes a real track. The makeDoc carries the detail in the meantime, including what was deliberately left undone. diff --git src/hg/makeDb/scripts/hprcPclai/hprcPclaiMakeTrackDb.py src/hg/makeDb/scripts/hprcPclai/hprcPclaiMakeTrackDb.py new file mode 100755 index 00000000000..15c2f8dde0f --- /dev/null +++ src/hg/makeDb/scripts/hprcPclai/hprcPclaiMakeTrackDb.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python3 +"""Write the trackDb stanzas for the hg38 hprcPclai composite track. + +One subtrack per HPRC Release 2 haplotype, read from the pcLAI GRCh38-coordinate +index CSV, so the track follows whatever the index lists rather than a hardcoded +sample list. Subtracks are grouped by sample and haplotype so hgTrackUi draws a +sample-by-haplotype matrix instead of a flat list of several hundred checkboxes. + +detailsScript and metaDataUrl are repeated on every subtrack rather than set once +on the composite. Both are read as local settings at the point they are used: +bigBedClick.c collects detailsScript.* with trackDbLocalSettingsWildMatch(), and +hgTrackUi's file fetch looks the path up on the tdb of the clicked track, so a +value on the parent alone would not be seen. + + hprcPclaiMakeTrackDb.py <index.csv> <defaultOnFile> <out.ra> + +defaultOnFile holds one "<sample>.<haplotype>" per line, the subtracks that start +out checked; blank lines and # comments are ignored. +""" +import sys +import csv + +TRACK = "hprcPclai" +GBDB = "/gbdb/$D/hprcPclai" +DATA_VERSION = "HPRC Release 2, pcLAI v1.1 (GRCh38 coordinates)" + +# Only four coordinates ever appear in the centroid field across all haplotypes: +# the four ancestry clusters pcLAI discretizes to. They would make a good +# categorical filter, but each value contains a comma and comma is the separator +# filterValues uses, so filtering on them needs a comma-free label field first. + +DETAILS_SCRIPT = ('{"dataUrl":"pclaiRefPanel.json","exportFields":["pca","centroid"],' + '"title":"Position in ancestry space","xLabel":"PC1","yLabel":"PC2"}') + + +def readIndex(fname): + """Return [(sample, haplotype)] from the HPRC index CSV, sorted by sample.""" + rows = [] + with open(fname, encoding="utf-8", newline="") as fh: + for row in csv.DictReader(fh): + samp = (row["sample_id"] or "").strip() + hap = (row["haplotype"] or "").strip() + if samp and hap: + rows.append((samp, hap)) + if not rows: + sys.exit("no rows in %s" % fname) + return sorted(set(rows)) + + +def readDefaultOn(fname): + on = set() + with open(fname, encoding="utf-8") as fh: + for line in fh: + line = line.split("#")[0].strip() + if line: + on.add(line) + return on + + +def main(): + if len(sys.argv) != 4: + sys.exit(__doc__) + idxFname, onFname, outFname = sys.argv[1:4] + haps = readIndex(idxFname) + wanted = readDefaultOn(onFname) + + samples = sorted({s for s, _h in haps}) + hapVals = sorted({h for _s, h in haps}) + + out = open(outFname, "w", encoding="utf-8") + w = out.write + w("# hg38 pcLAI local ancestry, one subtrack per HPRC Release 2 haplotype.\n") + w("# Generated by hg/makeDb/scripts/hprcPclai/hprcPclaiMakeTrackDb.py -- do not\n") + w("# hand-edit; see hg/makeDb/doc/hg38/hprcPclai.txt.\n\n") + + w("track %s\n" % TRACK) + w("compositeTrack on\n") + w("shortLabel pcLAI Ancestry\n") + w("longLabel Point cloud local ancestry inference (pcLAI) along HPRC assembly haplotypes\n") + w("group hprc\n") + w("type bigBed 9 +\n") + w("itemRgb on\n") + w("visibility dense\n") + w("priority 30\n") + # A whole chromosome holds a few thousand windows per haplotype, well over the + # 1000-item default at which pack mode gives up drawing; the block structure + # this track exists to show is only visible at that zoom. + w("maxItems 40000\n") + w("subGroup1 sample Sample %s\n" % " ".join("%s=%s" % (s, s) for s in samples)) + w("subGroup2 hap Haplotype %s\n" % " ".join("h%s=%s" % (h, h) for h in hapVals)) + w("sortOrder sample=+ hap=+\n") + w("dimensions dimensionY=sample dimensionX=hap\n") + w("dragAndDrop subTracks\n") + w("dataVersion %s\n" % DATA_VERSION) + w("\n") + + # Explicit priorities: without them the subtracks come out in an arbitrary + # (in practice reversed) order in the image, and with several hundred rows a + # stable sample order is what makes the display readable. + nOn = 0 + for prio, (samp, hap) in enumerate(haps, start=1): + key = "%s.%s" % (samp, hap) + sub = "%s%sH%s" % (TRACK, samp.replace("-", ""), hap) + on = key in wanted + nOn += on + w(" track %s\n" % sub) + w(" parent %s %s\n" % (TRACK, "on" if on else "off")) + w(" subGroups sample=%s hap=h%s\n" % (samp, hap)) + w(" priority %d\n" % prio) + w(" shortLabel %s\n" % key) + w(" longLabel pcLAI ancestry of HPRC assembly %s haplotype %s\n" % (samp, hap)) + w(" type bigBed 9 +\n") + w(" bigDataUrl %s/%s.bb\n" % (GBDB, key)) + w(" itemRgb on\n") + w(" mouseOver ${haplotype} ${window}<br>Window PC1,PC2: ${pca}" + "<br>Ancestry centroid: ${centroid}<br>Confidence: ${score}\n") + w(" detailsScript.scatterPlot.pca %s\n" % DETAILS_SCRIPT) + w(" metaDataUrl %s/pclaiRefPanel.json\n" % GBDB) + w("\n") + out.close() + + missing = wanted - {"%s.%s" % (s, h) for s, h in haps} + if missing: + sys.stderr.write("WARNING: default-on not in index: %s\n" % ", ".join(sorted(missing))) + sys.stderr.write("%s: %d subtracks, %d samples, %d on by default\n" + % (outFname, len(haps), len(samples), nOn)) + + +if __name__ == "__main__": + main()