# 2026-09-07 Claude max: Fiber-seq accessibility, FIRE peaks and CpG methylation, refs #36210

# The Stergachis and Vollger labs (UW / Utah) built their own hg38 track hub for
# Fiber-seq and asked for it as a native track:
#   https://fiberseq.github.io/UCSC-Fiber-seq-hub/hub.txt
# Shane Neph then sent 41 per-sample CpG methylation hubs to include alongside it
# (attachment cpg-hprc.txt on the "Fiber-seq in the genome browser" thread,
# 2026-09-07).  The CpG sample list is the same 41 samples as the accessibility
# compendium, so both became faceted composites under one container.

# The data lives on the UW Kopah S3 (Ceph) server under a per-sample hash
# directory.  We mirror it rather than pointing bigDataUrl at UW, since that
# server has had availability problems (see MLQ #37961) and a native track
# should not depend on it.

# ---------------------------------------------------------------------------
# Sample list
# ---------------------------------------------------------------------------

# The accession, sample name, cell type and per-sample S3 hash directory are
# checked in as
#   ~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqSamples.tsv
# The accessions and hashes were extracted from the lab's hub trackDb, and the
# sample names and cell types from the metadata TSV that hub points at:

cd /hive/data/genomes/hg38/bed/fiberSeq
wget -O hubTrackDb.txt https://fiberseq.github.io/UCSC-Fiber-seq-hub/hg38/trackDb.txt
wget -O hubMeta.tsv https://fiberseq.github.io/UCSC-Fiber-seq-hub/hg38/fireCompendium_metadata.tsv

# accession -> hash, from the bigDataUrls in the hub
grep -o "hashed.PacBio-Fiber-seq/[A-Z0-9]*/[0-9a-f]*" hubTrackDb.txt \
  | sed 's|hashed.PacBio-Fiber-seq/||' | tr '/' '\t' | sort -u > acc_hash.tsv
wc -l acc_hash.tsv
# 41 acc_hash.tsv

# accession -> sample, cellType, one row per sample
tail -n +2 hubMeta.tsv | awk -F'\t' '!seen[$4]++ {print $4"\t"$1"\t"$2}' | sort > acc_sample_type.tsv

join -t$'\t' acc_sample_type.tsv <(sort acc_hash.tsv) > manifest.tsv
wc -l manifest.tsv
# 41 manifest.tsv

# Cross-check that the CpG hubs Shane sent cover exactly these 41 accessions and
# agree on the sample names.  cpg-hprc.txt is his attachment, sample<TAB>hubUrl.
sed 's|.*ucsc-hprc-cpg/\([A-Z0-9]*\)/.*|\1|' cpg-hprc.txt | grep . | sort -u > cpg_acc.txt
comm -3 <(cut -f1 manifest.tsv) cpg_acc.txt
# no output: identical accession sets

join -t$'\t' -1 1 -2 1 <(cut -f1,2 manifest.tsv) \
  <(awk -F'\t' 'NF==2{sub(/.*ucsc-hprc-cpg\//,"",$2); split($2,a,"/"); print a[1]"\t"$1}' cpg-hprc.txt | sort) \
  | awk -F'\t' '$2!=$3'
# no output: all 41 sample labels agree

# manifest.tsv was then committed as fiberSeqSamples.tsv with a header line.

# ---------------------------------------------------------------------------
# Mirror the data
# ---------------------------------------------------------------------------

# 12 files per sample, 492 files, 576 GB:
#   all.percent.accessible.bw, hap1/hap2.percent.accessible.bw, fire-peaks.bb
#   all.nucleosome.coverage.bw
#   cpg.combined.bw, cpg.hap1.bw, cpg.hap2.bw
#   cpg.diffs_all.bw, cpg.diffs_p0.01.bw, cpg.diffs_p0.001.bw, cpg.diffs_p0.0001.bw
# The script resumes, verifies each file against the size the server reports and
# skips files that are already complete, so it can be re-run safely.  It takes an
# optional list of accessions, to refresh one sample without walking all 41.

~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqDownload.sh \
    /hive/data/genomes/hg38/bed/fiberSeq 12 > download.log 2>&1

grep -c '^got' download.log
grep -c FAIL download.log

# Check every file parses and holds data:
~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqCheck.sh /hive/data/genomes/hg38/bed/fiberSeq

# On the first mirror this reported two placeholder files at the source:
#   PM00001/hap1.percent.accessible.bw  (512 bytes, basesCovered 1)
#   PM00001/hap2.percent.accessible.bw  (512 bytes, basesCovered 1)
# They were valid bigWigs covering a single base with a value of zero, so
# GM12878's haplotype accessibility overlay drew nothing.  Note that they were
# NOT basesCovered 0, which is why the check tests for less than 1000; the
# smallest legitimate haplotype file here covers 399,383,234 bases (PM00008
# hap1, the near-haploid Hap1 line, which has little phaseable heterozygosity;
# every other sample's haplotype files cover 2.1 billion bases or more).
# Reported to the lab, and fixed by them - see the PM00001 section below.

# ---------------------------------------------------------------------------
# 2026-09-14  PM00001 (GM12878) reissued
# ---------------------------------------------------------------------------

# The lab reprocessed GM12878 and replaced the files IN PLACE, under the same
# per-sample hash directory, so fiberSeqSamples.tsv did not change.  Ten of its
# twelve files differ from the first mirror; only cpg.combined.bw is unchanged.
# A size sweep over all 41 samples confirmed PM00001 is the only one affected.

~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqDownload.sh \
    /hive/data/genomes/hg38/bed/fiberSeq 10 PM00001 > download.PM00001.log 2>&1

# The two haplotype files are now real: basesCovered 2,894,029,052 (hap1) and
# 2,382,899,132 (hap2), against 2,897,956,157 for the combined signal.
#
# This is what forced the change to how the downloader writes.  It used
# "curl -C - -o <finalFile>", which resumes from whatever is already on disk.
# That is correct for an interrupted transfer and silently corrupting when the
# server has REPLACED the file: curl would have appended the tail of the new
# 5.2 GB hap1 file to the 512-byte stub and reported success, and the size check
# afterwards would have passed because the total matches.  The script now
# fetches into <file>.part and moves it into place only when complete, so a
# resume can never land on a finished file, and the final mv is atomic, so
# hgTracks never reads a half-written file out of /gbdb.
#
# The peak file was reissued too but still has the fieldCount 3 header, so the
# rebuild below is still required.  Its peak calls changed substantially,
# 429,883 source peaks before and 196,742 now, which is a different call set
# rather than a reformatting.  Noted on fiberSeqCompendium.html, since figures
# made from the first version of the track will not reproduce for this sample.

# ---------------------------------------------------------------------------
# Nucleosome density: mirrored, NOT displayed
# ---------------------------------------------------------------------------

# Sent separately in September 2026, as all.nucleosome.coverage.bw under the
# same per-sample trackHub/bw directory.  It is NOT in the lab's own hub.
# Present for all 41 samples, 141 GB, added to the file list in the download
# and check scripts and mirrored with a plain re-run:

~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqDownload.sh \
    /hive/data/genomes/hg38/bed/fiberSeq 12 > download.nucleosome.log 2>&1

# The lab then asked us to hold the display back while they settle internally
# on how they want nucleosomes shown, and to keep the track set on FIRE peaks
# and CpG methylation.  So the files are mirrored and checked but no subtrack
# is generated: INCLUDE_NUC in fiberSeqTrackDb.py is False.  Setting it True and
# regenerating is the whole of turning it on, no download needed.
#
# What that flag builds, when the lab comes back: it is read depth, not a
# percentage, which makes it the one wiggle here that cannot take fixed
# viewLimits.  The genome-wide mean runs from 25 (PS00971) to 142 (PM00001) with
# sequencing depth, and single loci reach 1.7e5, so it is drawn with autoScale
# on and the axis is comparable within a sample but not between samples.
#
# The flag also moves the per-sample priorities, which is why they are assigned
# from a running counter rather than a fixed index per data type.  With the flag
# off the generated fiberSeq.ra is byte-identical to the version before the
# nucleosome work, and with it on to the version that added it, which is the
# check that the counter renumbers correctly either way.

~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqCheck.sh /hive/data/genomes/hg38/bed/fiberSeq
# no lines: all 492 files present, parsing and non-empty

# ---------------------------------------------------------------------------
# Rebuild the FIRE peak bigBeds
# ---------------------------------------------------------------------------

# The peak files carry narrowPeak data and even embed the narrowPeak autoSql,
# but their bigBed header records a field count of 3:
bigBedInfo -as PM00001/fire-peaks.bb | head -25
#   fieldCount: 3
#   definedFieldCount: 3
#   extraFieldCount: 0
#   as: table fire_peaks ... 10 columns ...
bigBedToBed PM00001/fire-peaks.bb -chrom=chr1 -start=1000000 -end=1100000 stdout | head -1
#   chr1 1000104 1000268 peak-247 224 . 22.421965 -1 31.751720958940915 94
# i.e. the data really has all 10 columns.  All 41 files are like this.
#
# With fieldCount 3, bbFieldIndex() cannot find signalValue or qValue, so no
# filter, no mouseOver and no details-page column can work.  Worse for this
# type: bigNarrowPeakLoadItems() calls bigBedMakeNumberFilter() for score,
# signalValue, pValue and qValue on every draw, and that calls getFieldNum()
# before it checks whether a filter was even requested, so hgTracks would
# errAbort with "error building filter with field signalValue.  Field not
# found." on any view of these files.  The rebuild is therefore required, not
# cosmetic.  It changes the header only, no data.

~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqFixPeaks.sh \
    /hive/data/genomes/hg38/bed/fiberSeq 8

# Peaks were called against the GRCh38 analysis set, which carries the
# Epstein-Barr virus decoy.  chrEBV is not in UCSC hg38 (not in chrom.sizes and
# not in chromInfo), so bedToBigBed refuses those rows and the script drops
# them, counting each one.  chrEBV is the only contig in the peak files that
# hg38 does not have; the bigWigs also contain chrEBV but need no rebuild,
# since hgTracks simply never queries a chrom it does not know.
#
# Reconciliation, from peakCounts.tsv (accession, source peaks, peaks in the
# track, peaks dropped on chrEBV):
#   source peaks     9,253,902
#   peaks in track   9,253,593
#   dropped chrEBV         309   in 20 of the 41 samples, 2 to 54 each
# Every sample balances exactly: source - chrEBV == in track.  The largest
# single loss is GM12878 with 54 peaks, which is 0.03 percent of its 196,742.
# This is mentioned in the methods section of fiberSeqCompendium.html.
# (Totals as of the PM00001 reissue; before it they were 9,487,043 source,
# 9,486,622 in track and 421 on chrEBV, all of the difference in that sample.)

# The rebuild also rounds signalValue and qValue to 3 decimals.  They arrive
# with full double precision (qValues like 22.807437495448326), which reads
# badly in a mouseover and on the details page and costs a third of the file
# size: 500 MB over the 41 files before rounding, 320 MB after (recounted after
# the PM00001 reissue).  pValue is the sentinel -1 throughout and is untouched.
#
# The rebuilt files are named fire-peaks.ucsc.bb and are what trackDb points
# at; the originals are kept next to them for comparison.

# Reconciliation is re-derived by the script on every run, into peakFixReport.txt:
sed -E 's/^([A-Z0-9]+) ([0-9]+) peaks \(source ([0-9]+), dropped ([0-9]+) on chrEBV\)$/\1\t\2\t\3\t\4/' \
    peakFixReport.txt | awk -F'\t' '{kept+=$2; src+=$3; ebv+=$4; if($3-$4!=$2) mism++}
    END{printf "kept %d source %d chrEBV %d mismatched %d\n", kept, src, ebv, mism+0}'
# kept 9253593 source 9253902 chrEBV 309 mismatched 0
#
# fixOne() skips a sample whose fire-peaks.ucsc.bb is newer than its source, so
# after the PM00001 reissue only that one rebuilt and the other 40 said "have".
# That also means the run only prints one "got" line, and peakFixReport.txt has
# to have PM00001's row replaced rather than being rewritten wholesale.

# Chromosome naming is otherwise already UCSC style, no chromAlias needed:
bigWigInfo -chroms PM00001/all.percent.accessible.bw | head -12
# chromCount 25, chr1 ... chrY (primary chromosomes only)
bigWigInfo -chroms PM00001/cpg.combined.bw | head -12
# chromCount 195, primary plus alts and patches

# ---------------------------------------------------------------------------
# trackDb, metadata and colors
# ---------------------------------------------------------------------------

# One script writes the whole track stanza plus the two faceted-composite
# metadata and color files, so the 584 stanzas stay consistent:

~/kent/src/hg/makeDb/scripts/fiberSeq/fiberSeqTrackDb.py

# It writes
#   ~/kent/src/hg/makeDb/trackDb/human/hg38/fiberSeq.ra
#   /hive/data/genomes/hg38/bed/fiberSeq/fiberSeqCompendium_metadata.tsv
#   /hive/data/genomes/hg38/bed/fiberSeq/fiberSeqCompendium_colors.json

# Structure:
#   fiberSeq            container, group regulation
#     fiberSeqAcc         multiWig overlay of 7 common cell lines, shown by default
#     fiberSeqCompendium  faceted composite, 41 samples, six data types:
#                         acc, peaks, hap, cpg, cpgHap, cpgDiff
#                         (a seventh, nuc, is built but held back: see below)
#
# Accessibility and methylation started as two composites, fiberSeqCompendium and
# fiberSeqMeth, and were merged into one.  They cover the identical 41 samples,
# and cartDump.c assigns priority with the data element as the OUTER loop and the
# data type as the inner one, so a single composite keeps a sample's six subtracks
# contiguous in the image.  As two composites the display was an accessibility
# block followed by a methylation block, so comparing the two assays for one
# sample meant reading across every other sample - which is the whole point,
# since both come off the same molecules in the same experiment.
#
# The composite uses the same faceted-composite machinery as Methbase
# (methbase2.ra): metaDataUrl for the sample table, colorSettingsUrl for the
# facet swatches, primaryKey for the row id, and dataTypes for the per-sample
# data kinds.  Two things about that machinery drove the naming:
#
# 1. Subtracks must be named <composite>_<accession>_<dataType>, with the
#    accession as the ONLY middle component.  facetedCompositeUi() in
#    hg/hgTrackUi/hgTrackUi.c cuts the data element at the first underscore
#    after the composite name, and cartDump.c rebuilds subtrack names as
#    <composite>_<element>_<type>.  The lab's hub used
#    fireCompendium_GM12878_PM00001_all, which would resolve to element
#    "GM12878" and then look for a track named fireCompendium_GM12878_all that
#    does not exist.  Sample name and cell type live in the metadata TSV.
#
# 2. dataTypes brings onlyVisibility with it, which is what lets the bigWig
#    subtracks default to full while the peak subtracks default to dense.  A
#    plain composite shares one visibility across all its subtracks, which is
#    the problem the lab ran into when they asked for peaks in dense and signal
#    in full.
#
# 3. The peaks are declared bigNarrowPeak, so the point-source offset in the
#    tenth column is drawn as a tick inside the peak (lfFromEncodePeak() sets
#    tallStart/tallEnd from it).  Getting there needed two small additions to
#    the tree, since bigNarrowPeak had been left behind by the bigBed-like
#    types:
#      hg/hgTracks/encode.c   bigNarrowPeakLoadItems() now honours mouseOver
#                             and mouseOverField, using the existing
#                             mouseOverSetupForBbi() / mouseOverGetBbiText()
#                             helpers in hg/hgTracks/mouseOver.c
#      hg/makeDb/trackDb/tagTypes.tab
#                             bigNarrowPeak added to mouseOver, scoreFilter,
#                             scoreFilterLimits, scoreMin, scoreMax,
#                             signalFilter and signalFilterLimits, and
#                             pValueFilter/qValueFilter (plus their Limits)
#                             registered for the first time.  encodePeakCfgUi()
#                             in hg/lib/hui.c and bigNarrowPeakLoadItems() both
#                             already handled the q/pValue filters; they had
#                             just never been declared here, so tdbQuery
#                             -strict rejected them.
#
# 4. The peak filters are the ENCODE peak settings, signalFilter, qValueFilter
#    and scoreFilter with their *Limits, which encodePeakCfgUi() draws with its
#    own labels ("Minimum Signal value", "Minimum Q-Value (-log10)").  The
#    bigBed-generic filter.<field> settings are NOT read by this type.  All
#    three defaults are the full range, so nothing is hidden until the user
#    narrows one.  Verified in the rendered UI:
#      hgTrackUi?db=hg38&g=fiberSeqCompendium_PM00004_peaks
#    draws min/max boxes for all three with the right limit hints.
#
# 5. Facet columns: only "sampleClass" is faceted.  facetedComposite.js offers
#    a facet value only when it occurs more than once (a checkbox matching a
#    single row is just a slow search box), so:
#      accession  41 distinct, all count 1, and excluded anyway as primaryKey
#      _sample    41 distinct, all count 1 - can never be a facet
#      _cellType  14 distinct but only 2 with count > 1 (Lymphoblastoid 27,
#                 Embryonic stem cell 2), so as a facet it drew two checkboxes
#                 and left 12 samples unreachable.  Underscored, so it is a
#                 searchable and sortable column instead.
#      sampleClass 2 values, HPRC (20) and Common Cell Line (21), both count > 1.
#                 Read from the fifth column of fiberSeqSamples.tsv, which
#                 carries the lab's own classification (the Note column of
#                 Mitchell's sample sheet).  It used to be derived from the
#                 free-text cell type, which misfiled five lymphoblastoid lines
#                 that are not HPRC - see the note below.  Still worth extending
#                 when the lab sends real donor metadata (sex, population), which
#                 would give more than one useful facet.
#
# 6. Subtracks need an explicit priority.  Without one they fall back to a label
#    sort, which showed a sample's data types as Peaks, CpG, Acc on a first
#    visit.  The script now numbers them sample-outer, declared-data-type-inner
#    (i*10 + j + 1), which matches the row of data type checkboxes across the
#    top of the table.  cartDump.c clears "<mdid>_*.priority" and writes its own
#    on every submit, so this only sets the starting order.
#
#    The children INSIDE each multiWig need their own priority too, and there it
#    is not cosmetic.  makeContainerTrack() in hg/hgTracks/container.c sorts a
#    container's children with trackPriCmp, which compares priority alone and
#    returns 0 on a tie, and slSort is not stable - so children left on the
#    parent's inherited priority draw in an arbitrary order.  For cpgDiff that
#    decides the picture: it is a solidOverlay and all four files hold the same
#    value at a shared base, so the level drawn last is the colour the user sees.
#    Left unset, p<0.01 painted over p<0.001 and p<0.0001, i.e. the least
#    stringent tier won and the track said the opposite of what it meant.
#    Measured at chr20:29,300,000-29,305,000, where 107 of the 133 p<0.01 bases
#    are also p<0.0001: before, 98 image columns showed only yellow and 2 showed
#    red; after giving the four levels priority 1-4, 88 columns show red and 51
#    show grey.  The two haplotype children carry priority 1 and 2 for the same
#    reason, though transparentOverlay makes their order much less visible.
#
#    The haplotype children also take their shortLabel prefix from the
#    container's own shortLabel, so the CpG pair reads "<sample> CpG Hap1/2"
#    rather than colliding with the accessibility pair's "<sample> Hap1/2".
#
# 7. mouseOver text only appears in pack or full.  Dense makes no per-item map
#    boxes at all, so no bigBed-like track has a per-item hover there.  Since
#    Andrew asked for peaks in dense by default, the mouseover is there for
#    whoever switches a peak track to pack.  Checked by temporarily dropping
#    onlyVisibility and reading the image map: in pack the AREA tags carry
#      data-tooltip='<b>K562 FIRE peak</b><br>FIRE score: 17.76015
#                    <br>-log10 FDR: 22.807437495448326<br>Score: 177'
#    and in dense there are no per-peak AREA tags.  Note the tooltip is HTML
#    entity encoded in the page, so grep for it with the entities in mind.

# ---------------------------------------------------------------------------
# Rendering checks
# ---------------------------------------------------------------------------

# ACTB promoter, the accessibility overlay, all seven cell lines in color:
#   hgRenderTracks?db=hg38&position=chr7:5,527,000-5,533,000&hideTracks=1
#     &fiberSeq=show&fiberSeqAcc=full
# Note that a child of a container needs the container shown as well; without
# fiberSeq=show only the ruler comes back.
#
# The compendium at the same locus, showing that onlyVisibility works: the
# three default samples come up with peaks in dense and signal in full, in one
# composite.
#   ...&fiberSeqCompendium=full
#
# GNAS, an imprinted locus, for the methylation composite.  GM12878's
# haplotype overlay switches from orange (haplotype 2 methylated) to blue
# (haplotype 1 methylated) across the locus, and the difference track follows:
#   hgRenderTracks?db=hg38&position=chr20:58,838,000-58,860,000&hideTracks=1
#     &fiberSeq=show&fiberSeqCompendium=full
#     &fiberSeqCompendium_PM00001_cpgHap_sel=1
#     &fiberSeqCompendium_PM00001_cpgDiff_sel=1
#
# The merged composite at the same locus, with the three default data types on,
# renders each sample's tracks together and in the declared order:
#   GM12878 Acc / GM12878 Peaks / GM12878 CpG / K562 Acc / K562 Peaks / ...
#
# The faceted table's embedded JSON on the composite UI page confirms the
# naming works: dataElements comes back as bare accessions ("PM00001",
# "PM00004", "PM00005") and dataTypes as acc/peaks active, hap inactive.
