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/hprcPclaiMakeBb.sh src/hg/makeDb/scripts/hprcPclai/hprcPclaiMakeBb.sh new file mode 100755 index 00000000000..626e0c86aa1 --- /dev/null +++ src/hg/makeDb/scripts/hprcPclai/hprcPclaiMakeBb.sh @@ -0,0 +1,46 @@ +#!/bin/bash +# Convert one pcLAI GRCh38-coordinate BED into a UCSC bigBed. +# +# The source is a bed9+1: the name column packs three things +# ("SAMPLE/hN/<window>_(PC1,PC2)") and column 10 is the source's "centroid" +# column, the discretized ancestry of the window written as the PCA centroid of +# that ancestry cluster (hence only four distinct values across the whole +# release). Those are split into their own named fields (haplotype, window, pca, +# centroid) so they can be shown on mouseover and fed to the details-page +# scatterplot, and the name column is blanked: the +# packed string is too long to draw as an item label and carries no extra +# information once the fields are split out. +# +# thickStart/thickEnd are rewritten to the item bounds. The source format +# specifies thickStart == chromStart and thickEnd == chromEnd; thickEnd always +# holds, and thickStart holds except in +# 54811 of the 11936603 windows (0.46%), where it is chromStart-1 and so outside +# the item -- an off-by-one in the projection that bedToBigBed rejects. The two +# columns carry no information for this annotation, so they are set to the item +# bounds rather than dropping those windows. +# +# Usage: hprcPclaiMakeBb.sh <in.bed> <out.bb> <chrom.sizes> +set -u -o pipefail +inBed=$1; outBb=$2; sizes=$3 +SCR=$HOME/kent/src/hg/makeDb/scripts/hprcPclai +tmp=$(mktemp "${TMPDIR:-/data/tmp}/hprcPclai.XXXXXX.bed"); trap 'rm -f "$tmp"' EXIT + +awk -F'\t' 'BEGIN{OFS="\t"} + { n=split($4,a,"/"); + if (n!=3) { bad++; next } + hap=a[1]"/"a[2]; rest=a[3]; + k=split(rest,b,"_"); pca=b[k]; + win=(pca!="") ? substr(rest,1,length(rest)-length(pca)-1) : rest; + if ($7!=$2) fixedThick++; + print $1,$2,$3,"",$5,$6,$7=$2,$8=$3,$9,hap,win,pca,$10 } + END{ if (bad) print "UNPARSED_NAMES "bad" in '"$inBed"'" > "/dev/stderr"; + if (fixedThick) print "FIXED_THICKSTART "fixedThick" in '"$inBed"'" > "/dev/stderr" }' "$inBed" \ + | LC_COLLATE=C sort -k1,1 -k2,2n > "$tmp" + +inCount=$(wc -l < "$inBed") +outCount=$(wc -l < "$tmp") +[ "$inCount" = "$outCount" ] || { echo "COUNT_MISMATCH $inBed $inCount -> $outCount" >&2; exit 4; } + +bedToBigBed -type=bed9+4 -tab -as="$SCR/hprcPclai.as" "$tmp" "$sizes" "$outBb" \ + || { echo "BB_FAIL $inBed" >&2; exit 6; } +echo "OK $outBb $outCount"