ce780dd2f1216ce728ab6bb69ac19a39ddc694fd max Tue Sep 8 00:26:39 2026 -0700 hg38: Fiber-seq container with accessibility, FIRE peaks and CpG methylation, 41 samples Native version of the Stergachis/Vollger lab hub at https://fiberseq.github.io/UCSC-Fiber-seq-hub/hub.txt, plus the per-sample CpG methylation Shane Neph asked to have alongside it. Both cover the same 41 samples: 14 cell lines and 27 lymphoblastoid lines from HPRC and GIAB individuals. fiberSeq container, group regulation fiberSeqAcc multiWig overlay of 7 common cell lines, on by default fiberSeqCompendium faceted composite, dataTypes acc/peaks/hap fiberSeqMeth faceted composite, dataTypes comb/hap/diffs, "Methylation" Both composites use the Methbase faceted-composite machinery. Subtracks are named <composite>_<accession>_<dataType> with the accession as the only middle component, because facetedCompositeUi() cuts the data element at the first underscore and cartDump.c reassembles the name from the pieces; the hub's <composite>_<sample>_<accession>_<type> names would have resolved to tracks that do not exist. Sample name and cell type live in the metadata TSV instead. Using dataTypes also brings onlyVisibility, which is what lets the peaks default to dense while the signal tracks default to full, the mixed-visibility default Andrew Stergachis asked for. 397 GB mirrored from the UW Kopah S3 server rather than pointed at over the network, since a native track should not depend on it. The FIRE peak bigBeds had to be rebuilt: they carry full narrowPeak data but their header records a field count of 3, which hides signalValue and qValue from the browser and would have made hgTracks errAbort in bigNarrowPeakLoadItems(). The rebuild fixes the header and rounds the two float columns to 3 decimals, 467 MB to 313 MB. It drops 421 of 9,487,043 peaks called on chrEBV, the EBV decoy of the GRCh38 analysis set, which hg38 does not have; 9,486,622 remain and every sample reconciles exactly. Reported upstream, along with GM12878's two haplotype accessibility bigWigs, which are one-base placeholders at the source. refs #36210 diff --git src/hg/makeDb/scripts/fiberSeq/fiberSeqFixPeaks.sh src/hg/makeDb/scripts/fiberSeq/fiberSeqFixPeaks.sh new file mode 100755 index 00000000000..05f5e521178 --- /dev/null +++ src/hg/makeDb/scripts/fiberSeq/fiberSeqFixPeaks.sh @@ -0,0 +1,108 @@ +#!/bin/bash +# Rewrite the FIRE peak bigBeds with a correct header. +# +# The files on the UW server carry narrowPeak data (10 columns: name, score, +# strand, signalValue, pValue, qValue and the point-source offset) and even +# embed the narrowPeak autoSql, but their bigBed header records a field count +# of 3. bbFieldIndex() therefore cannot find signalValue or qValue, so no +# filter, no mouseOver and no details-page column works on them, and asking +# hgTracks for a signalValue filter aborts with "Field not found". +# +# This converts each file back to BED and rebuilds it as a bed6+4 bigBed with +# the narrowPeak schema, which changes no data, only the header. Reported +# upstream; drop this step once the lab regenerates the files. +# +# Usage: fiberSeqFixPeaks.sh <dataDir> [jobs] + +set -o pipefail + +dataDir=$1 +jobs=${2:-8} +if [ -z "$dataDir" ]; then + echo "usage: $0 <dataDir> [jobs]" >&2 + exit 1 +fi + +scriptDir=$(dirname "$(readlink -f "$0")") +sampleList=$scriptDir/fiberSeqSamples.tsv +asFile=$(readlink -f ~/kent/src/hg/lib/encode/narrowPeak.as) +chromSizes=/hive/data/genomes/hg38/chrom.sizes + +for f in "$asFile" "$chromSizes"; do + [ -s "$f" ] || { echo "missing $f" >&2; exit 1; } +done + +fixOne() { + acc=$1; dataDir=$2; asFile=$3; chromSizes=$4 + src=$dataDir/$acc/fire-peaks.bb + out=$dataDir/$acc/fire-peaks.ucsc.bb + tmp=$dataDir/$acc/fire-peaks.$$.bed + + if [ ! -s "$src" ]; then + echo "MISSING $acc/fire-peaks.bb" >&2 + return 1 + fi + # Already converted and newer than the source: nothing to do. + if [ -s "$out" ] && [ "$out" -nt "$src" ]; then + echo "have $acc" + return 0 + fi + + if ! bigBedToBed "$src" "$tmp"; then + echo "FAIL $acc bigBedToBed" >&2 + rm -f "$tmp" + return 1 + fi + inRows=$(wc -l < "$tmp") + # Every row must have the 10 narrowPeak columns; stop rather than let + # bedToBigBed silently reinterpret a short row. + ragged=$(awk -F'\t' 'NF!=10' "$tmp" | wc -l) + if [ "$ragged" != "0" ]; then + echo "FAIL $acc $ragged of $inRows rows are not 10 columns" >&2 + rm -f "$tmp" + return 1 + fi + # The peaks were called against the GRCh38 analysis set, which carries the + # Epstein-Barr virus decoy. chrEBV is not part of UCSC hg38, so those + # peaks cannot be placed and are dropped here; the count is reported so it + # is never a silent loss. Any other unplaceable contig is an error. + unknown=$(cut -f1 "$tmp" | sort -u | grep -v -x -F -f <(cut -f1 "$chromSizes" | sort -u) | grep -v -x chrEBV) + if [ -n "$unknown" ]; then + echo "FAIL $acc unexpected contigs: $(echo $unknown | tr '\n' ' ')" >&2 + rm -f "$tmp" + return 1 + fi + ebv=$(awk -F'\t' '$1=="chrEBV"' "$tmp" | wc -l) + # signalValue and qValue come out of the pipeline with full double precision + # (qValue values like 22.807437495448326), which bloats the file and reads + # badly in a mouseover and on the details page. Round both to 3 decimals. + # pValue is the sentinel -1 throughout and is left exactly as it is. + awk -F'\t' -v OFS='\t' '$1!="chrEBV" { + $7 = sprintf("%.3f", $7); + $9 = sprintf("%.3f", $9); + print + }' "$tmp" | sort -k1,1 -k2,2n > "$tmp.s" && mv "$tmp.s" "$tmp" + if ! bedToBigBed -as="$asFile" -type=bed6+4 -tab "$tmp" "$chromSizes" "$out" 2>/dev/null; then + echo "FAIL $acc bedToBigBed" >&2 + rm -f "$tmp" "$out" + return 1 + fi + outRows=$(bigBedInfo "$out" | awk -F': ' '$1=="itemCount"{gsub(",","",$2); print $2}') + rm -f "$tmp" + want=$((inRows - ebv)) + if [ "$want" != "$outRows" ]; then + echo "FAIL $acc kept $outRows, expected $want (source $inRows, chrEBV $ebv)" >&2 + return 1 + fi + echo "got $acc $outRows peaks (source $inRows, dropped $ebv on chrEBV)" +} +export -f fixOne + +grep -v '^#' "$sampleList" | awk -F'\t' 'NF{print $1}' \ + | xargs -P "$jobs" -I{} bash -c 'fixOne "$1" "$2" "$3" "$4"' _ {} "$dataDir" "$asFile" "$chromSizes" + +echo +echo "field counts after conversion (all should be 10):" +for f in "$dataDir"/*/fire-peaks.ucsc.bb; do + bigBedInfo "$f" | awk -F': ' '$1=="fieldCount"{print $2}' +done | sort | uniq -c