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/fiberSeqTrackDb.py src/hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py
new file mode 100755
index 00000000000..cd38fe438b0
--- /dev/null
+++ src/hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py
@@ -0,0 +1,404 @@
+#!/usr/bin/env python3
+"""Generate the hg38 Fiber-seq trackDb stanza plus the faceted-composite
+metadata and color files.
+
+Writes, given the sample list in fiberSeqSamples.tsv:
+
+  <trackDbDir>/fiberSeq.ra                      the track stanzas
+  <dataDir>/fiberSeqCompendium_metadata.tsv     facet table for accessibility
+  <dataDir>/fiberSeqCompendium_colors.json      facet swatches
+  <dataDir>/fiberSeqMeth_metadata.tsv           facet table for methylation
+  <dataDir>/fiberSeqMeth_colors.json            facet swatches
+
+Subtrack names are deliberately "<composite>_<accession>_<dataType>" with the
+accession as the only middle component.  facetedCompositeUi() in
+hg/hgTrackUi/hgTrackUi.c derives the data element by cutting at the first
+underscore after the composite name, and cartDump.c rebuilds the subtrack name
+as composite_element_type, so an accession holding an underscore (or a middle
+part like "GM12878_PM00001") would make the browser look for tracks that do not
+exist.  Sample name and cell type live in the metadata table instead.
+
+Usage: fiberSeqTrackDb.py [--data-dir DIR] [--gbdb-dir DIR] [--trackdb-dir DIR]
+"""
+
+import argparse
+import json
+import os
+import sys
+
+SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+SAMPLE_LIST = os.path.join(SCRIPT_DIR, "fiberSeqSamples.tsv")
+
+# The seven cell lines the lab wants in the always-on overlay, in their order,
+# with the Okabe-Ito colors they chose in their own hub.
+DEFAULT_OVERLAY = [
+    ("PM00001", "230,159,0"),    # GM12878, orange
+    ("PM00004", "86,180,233"),   # K562, sky blue
+    ("PM00005", "0,158,115"),    # HepG2, bluish green
+    ("PM00010", "240,228,66"),   # H1, yellow
+    ("PM00008", "0,114,178"),    # Hap1, blue
+    ("PM00012", "213,94,0"),     # Hek293T, vermillion
+    ("PM00009", "204,121,167"),  # Jurkat, reddish purple
+]
+
+# Samples selected in the two faceted composites on a first visit.  A clean
+# cross-product, so the facet table comes up as a tidy grid rather than the
+# ragged per-sample mix the source hub had.
+DEFAULT_SELECTED = ["PM00001", "PM00004", "PM00005"]
+
+HAP1_COLOR = "0,114,178"    # Okabe-Ito blue
+HAP2_COLOR = "213,94,0"     # Okabe-Ito vermillion
+
+# CpG haplotype-difference thresholds: file suffix, label, color.  A sequential
+# yellow-to-red ramp over nested significance cutoffs, monotonic in lightness.
+DIFF_LEVELS = [
+    ("cpg.diffs_all.bw", "All", "137,143,143"),
+    ("cpg.diffs_p0.01.bw", "p < 0.01", "235,229,52"),
+    ("cpg.diffs_p0.001.bw", "p < 0.001", "245,148,22"),
+    ("cpg.diffs_p0.0001.bw", "p < 0.0001", "255,0,0"),
+]
+
+# Sample classes, derived from the lab's own free-text cell type.  27 of the 41
+# samples are lymphoblastoid, so cell type alone gives one useless bucket; this
+# splits them into three facet values.  Okabe-Ito colors for the swatches.
+SAMPLE_CLASS_COLORS = {
+    "Lymphoblastoid cell line": "#0072B2",
+    "Stem cell": "#009E73",
+    "Cancer or immortalized cell line": "#D55E00",
+}
+
+
+def sampleClass(cellType):
+    """Group a free-text cell type into one of three facet values."""
+    low = cellType.lower()
+    if "lymphoblastoid" in low:
+        return "Lymphoblastoid cell line"
+    if "stem cell" in low or "ipsc" in low:
+        return "Stem cell"
+    return "Cancer or immortalized cell line"
+
+
+def readSamples(path):
+    """Read fiberSeqSamples.tsv into a list of dicts, in file order."""
+    samples = []
+    with open(path) as f:
+        for line in f:
+            if line.startswith("#") or not line.strip():
+                continue
+            fields = line.rstrip("\n").split("\t")
+            if len(fields) < 4:
+                sys.exit("bad sample line, want 4 fields: %s" % line.rstrip())
+            acc, sample, cellType, _hash = fields[:4]
+            samples.append({
+                "accession": acc,
+                "sample": sample,
+                "cellType": cellType,
+                "sampleClass": sampleClass(cellType),
+            })
+    if not samples:
+        sys.exit("no samples read from %s" % path)
+    return samples
+
+
+def writeMetadata(path, samples):
+    """Facet table.  The first column is the primaryKey; plain column names get
+    facet checkboxes, a leading underscore means searchable but not faceted."""
+    with open(path, "w") as f:
+        f.write("accession\tsampleClass\tcellType\t_sample\n")
+        for s in samples:
+            f.write("%s\t%s\t%s\t%s\n" % (s["accession"], s["sampleClass"],
+                                          s["cellType"], s["sample"]))
+
+
+def writeColors(path):
+    with open(path, "w") as f:
+        json.dump({"sampleClass": SAMPLE_CLASS_COLORS}, f, indent=4)
+        f.write("\n")
+
+
+def stanza(indent, lines):
+    """Render one trackDb stanza at the given indent, with a trailing blank."""
+    pad = " " * indent
+    return "".join("%s%s\n" % (pad, l) for l in lines) + "\n"
+
+
+def accOverlay(gbdb, samples):
+    """The always-on overlay of the seven common cell lines."""
+    byAcc = {s["accession"]: s for s in samples}
+    out = stanza(4, [
+        "track fiberSeqAcc",
+        "parent fiberSeq",
+        "container multiWig",
+        "aggregate transparentOverlay",
+        "showSubtrackColorOnUi on",
+        "type bigWig 0 100",
+        "viewLimits 0:100",
+        "autoScale off",
+        "alwaysZero on",
+        "graphTypeDefault bar",
+        "windowingFunction maximum",
+        "maxHeightPixels 100:50:8",
+        "visibility full",
+        "priority 1",
+        "shortLabel Fiber-seq Acc",
+        "longLabel Fiber-seq percent-accessible chromatin in seven common cell lines",
+    ])
+    for acc, color in DEFAULT_OVERLAY:
+        s = byAcc[acc]
+        out += stanza(8, [
+            "track fiberSeqAcc_%s" % acc,
+            "parent fiberSeqAcc",
+            "type bigWig",
+            "bigDataUrl %s/%s/all.percent.accessible.bw" % (gbdb, acc),
+            "color %s" % color,
+            "shortLabel %s" % s["sample"],
+            "longLabel %s Fiber-seq percent-accessible chromatin, both haplotypes"
+            % s["sample"],
+        ])
+    return out
+
+
+def compendium(gbdb, dataUrlDir, samples):
+    """Faceted composite: percent accessible, FIRE peaks, haplotype overlay."""
+    out = stanza(4, [
+        "track fiberSeqCompendium",
+        "parent fiberSeq",
+        "compositeTrack faceted",
+        "type bigWig",
+        "shortLabel Fiber-seq Compendium",
+        "longLabel Fiber-seq percent accessible, FIRE peaks and haplotype overlays "
+        "in %d samples" % len(samples),
+        "metaDataUrl %s/fiberSeqCompendium_metadata.tsv" % dataUrlDir,
+        "colorSettingsUrl %s/fiberSeqCompendium_colors.json" % dataUrlDir,
+        "primaryKey accession",
+        'dataTypes acc|"Percent accessible" peaks|"FIRE peaks" '
+        'hap|"Haplotype accessibility"',
+        "defaultSortField accession",
+        "maxCheckboxes 50",
+        "noInherit on",
+        "visibility hide",
+        "priority 2",
+    ])
+    for s in samples:
+        acc, name = s["accession"], s["sample"]
+        accOn = "on" if acc in DEFAULT_SELECTED else "off"
+        peaksOn = "on" if acc in DEFAULT_SELECTED else "off"
+
+        out += stanza(8, [
+            "track fiberSeqCompendium_%s_acc" % acc,
+            "parent fiberSeqCompendium %s" % accOn,
+            "type bigWig",
+            "bigDataUrl %s/%s/all.percent.accessible.bw" % (gbdb, acc),
+            "shortLabel %s Acc" % name,
+            "longLabel %s Fiber-seq percent accessible, both haplotypes" % name,
+            "color 0,0,0",
+            "viewLimits 0:100",
+            "autoScale off",
+            "alwaysZero on",
+            "graphTypeDefault bar",
+            "windowingFunction maximum",
+            "maxHeightPixels 100:40:8",
+            "onlyVisibility full",
+        ])
+
+        # bigNarrowPeak, so the point-source offset in the tenth column is drawn
+        # as a tick inside the peak (lfFromEncodePeak() sets tallStart/tallEnd
+        # from it).  Its filters are the ENCODE peak settings, signalFilter and
+        # qValueFilter, which encodePeakCfgUi() already knows how to draw; the
+        # bigBed-generic filter.<field> settings are not read by this type.
+        # Both defaults are the full range, so nothing is hidden until the user
+        # narrows it.  signalValue tops out near 90 and the FIRE pipeline caps
+        # qValue at 100.  pValue is -1 throughout, so no filter is offered.
+        out += stanza(8, [
+            "track fiberSeqCompendium_%s_peaks" % acc,
+            "parent fiberSeqCompendium %s" % peaksOn,
+            "type bigNarrowPeak",
+            "bigDataUrl %s/%s/fire-peaks.ucsc.bb" % (gbdb, acc),
+            "shortLabel %s Peaks" % name,
+            "longLabel %s Fiber-seq FIRE peaks" % name,
+            "signalFilter 0",
+            "signalFilterLimits 0:100",
+            "qValueFilter 0",
+            "qValueFilterLimits 0:100",
+            "scoreFilter 0",
+            "scoreFilterLimits 0:1000",
+            "mouseOver <b>%s FIRE peak</b><br>FIRE score: ${signalValue}"
+            "<br>-log10 FDR: ${qValue}<br>Score: ${score}" % name,
+            "onlyVisibility dense",
+        ])
+
+        out += stanza(8, [
+            "track fiberSeqCompendium_%s_hap" % acc,
+            "parent fiberSeqCompendium off",
+            "container multiWig",
+            "aggregate transparentOverlay",
+            "showSubtrackColorOnUi on",
+            "type bigWig 0 100",
+            "viewLimits 0:100",
+            "autoScale off",
+            "alwaysZero on",
+            "windowingFunction maximum",
+            "maxHeightPixels 100:40:8",
+            "shortLabel %s Hap1/2" % name,
+            "longLabel %s Fiber-seq percent accessible, haplotype 1 (blue) and 2 (orange)"
+            % name,
+            "onlyVisibility full",
+        ])
+        for hap, color in (("h1", HAP1_COLOR), ("h2", HAP2_COLOR)):
+            n = hap[1]
+            out += stanza(12, [
+                "track fiberSeqCompendium_%s_hap_%s" % (acc, hap),
+                "parent fiberSeqCompendium_%s_hap" % acc,
+                "type bigWig",
+                "bigDataUrl %s/%s/hap%s.percent.accessible.bw" % (gbdb, acc, n),
+                "color %s" % color,
+                "shortLabel %s Hap%s" % (name, n),
+                "longLabel %s Fiber-seq percent accessible, haplotype %s" % (name, n),
+            ])
+    return out
+
+
+def methylation(gbdb, dataUrlDir, samples):
+    """Faceted composite: CpG methylation, combined and per haplotype."""
+    out = stanza(4, [
+        "track fiberSeqMeth",
+        "parent fiberSeq",
+        "compositeTrack faceted",
+        "type bigWig",
+        "shortLabel Methylation",
+        "longLabel CpG methylation from Fiber-seq reads, combined and by haplotype, "
+        "in %d samples" % len(samples),
+        "metaDataUrl %s/fiberSeqMeth_metadata.tsv" % dataUrlDir,
+        "colorSettingsUrl %s/fiberSeqMeth_colors.json" % dataUrlDir,
+        "primaryKey accession",
+        'dataTypes comb|"Combined CpG" hap|"Hap1/Hap2 CpG" '
+        'diffs|"Haplotype differences"',
+        "defaultSortField accession",
+        "maxCheckboxes 50",
+        "noInherit on",
+        "visibility hide",
+        "priority 3",
+    ])
+    for s in samples:
+        acc, name = s["accession"], s["sample"]
+        combOn = "on" if acc in DEFAULT_SELECTED else "off"
+
+        out += stanza(8, [
+            "track fiberSeqMeth_%s_comb" % acc,
+            "parent fiberSeqMeth %s" % combOn,
+            "type bigWig 0 100",
+            "bigDataUrl %s/%s/cpg.combined.bw" % (gbdb, acc),
+            "shortLabel %s CpG" % name,
+            "longLabel %s CpG methylation, both haplotypes" % name,
+            "color 0,0,0",
+            "viewLimits 0:100",
+            "autoScale off",
+            "windowingFunction mean",
+            "maxHeightPixels 100:40:8",
+            "onlyVisibility full",
+        ])
+
+        out += stanza(8, [
+            "track fiberSeqMeth_%s_hap" % acc,
+            "parent fiberSeqMeth off",
+            "container multiWig",
+            "aggregate transparentOverlay",
+            "showSubtrackColorOnUi on",
+            "type bigWig 0 100",
+            "viewLimits 0:100",
+            "autoScale off",
+            "windowingFunction mean",
+            "maxHeightPixels 100:40:8",
+            "shortLabel %s CpG Hap1/2" % name,
+            "longLabel %s CpG methylation, haplotype 1 (blue) and 2 (orange)" % name,
+            "onlyVisibility full",
+        ])
+        for hap, color in (("h1", HAP1_COLOR), ("h2", HAP2_COLOR)):
+            n = hap[1]
+            out += stanza(12, [
+                "track fiberSeqMeth_%s_hap_%s" % (acc, hap),
+                "parent fiberSeqMeth_%s_hap" % acc,
+                "type bigWig",
+                "bigDataUrl %s/%s/cpg.hap%s.bw" % (gbdb, acc, n),
+                "color %s" % color,
+                "shortLabel %s CpG Hap%s" % (name, n),
+                "longLabel %s CpG methylation, haplotype %s" % (name, n),
+            ])
+
+        out += stanza(8, [
+            "track fiberSeqMeth_%s_diffs" % acc,
+            "parent fiberSeqMeth off",
+            "container multiWig",
+            "aggregate solidOverlay",
+            "showSubtrackColorOnUi on",
+            "type bigWig -100 100",
+            "viewLimits -100:100",
+            "autoScale off",
+            "windowingFunction mean",
+            "maxHeightPixels 100:50:8",
+            "shortLabel %s CpG diffs" % name,
+            "longLabel %s CpG methylation difference between haplotypes, "
+            "by significance threshold" % name,
+            "onlyVisibility full",
+        ])
+        # Least significant first, so the more significant levels draw on top.
+        for i, (fname, label, color) in enumerate(DIFF_LEVELS):
+            out += stanza(12, [
+                "track fiberSeqMeth_%s_diffs_l%d" % (acc, i),
+                "parent fiberSeqMeth_%s_diffs" % acc,
+                "type bigWig",
+                "bigDataUrl %s/%s/%s" % (gbdb, acc, fname),
+                "color %s" % color,
+                "shortLabel %s %s" % (name, label),
+                "longLabel %s CpG haplotype difference, %s" % (name, label),
+            ])
+    return out
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__,
+                                 formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--data-dir", default="/hive/data/genomes/hg38/bed/fiberSeq",
+                    help="where the metadata and color files are written")
+    ap.add_argument("--gbdb-dir", default="/gbdb/hg38/fiberSeq",
+                    help="bigDataUrl prefix, i.e. the symlink directory")
+    ap.add_argument("--trackdb-dir",
+                    default=os.path.expanduser(
+                        "~/kent/src/hg/makeDb/trackDb/human/hg38"),
+                    help="where fiberSeq.ra is written")
+    args = ap.parse_args()
+
+    samples = readSamples(SAMPLE_LIST)
+
+    # The faceted composite fetches these two over http, from the same /gbdb
+    # path the browser serves, so trackDb refers to them the same way.
+    for name in ("fiberSeqCompendium", "fiberSeqMeth"):
+        writeMetadata(os.path.join(args.data_dir, "%s_metadata.tsv" % name), samples)
+        writeColors(os.path.join(args.data_dir, "%s_colors.json" % name))
+
+    raPath = os.path.join(args.trackdb_dir, "fiberSeq.ra")
+    with open(raPath, "w") as f:
+        f.write("# Fiber-seq: chromatin accessibility, FIRE regulatory elements and CpG\n"
+                "# methylation from PacBio HiFi Fiber-seq, Stergachis and Vollger labs.\n"
+                "# Generated by hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py.\n"
+                "# Do not edit by hand, edit the script and regenerate.\n\n")
+        f.write(stanza(0, [
+            "track fiberSeq",
+            "superTrack on show",
+            "shortLabel Fiber-seq",
+            "longLabel Fiber-seq chromatin accessibility, regulatory elements and CpG methylation",
+            "group regulation",
+            "priority 2.5",
+        ]))
+        f.write(accOverlay(args.gbdb_dir, samples))
+        f.write(compendium(args.gbdb_dir, args.gbdb_dir, samples))
+        f.write(methylation(args.gbdb_dir, args.gbdb_dir, samples))
+
+    nSub = len(DEFAULT_OVERLAY) + len(samples) * (3 + 2) + len(samples) * (3 + 2 + 4)
+    print("wrote %s" % raPath)
+    print("  %d samples, %d track stanzas" % (len(samples), nSub + 4))
+    print("  metadata and colors in %s" % args.data_dir)
+
+
+if __name__ == "__main__":
+    main()