b1fea807e3a0f6c6b0ebb6b1244b3291aeaaad28 max Mon Sep 14 07:51:13 2026 -0700 hg38 Fiber-seq makeDoc: drop the git attribution note from 8ceea19db11, refs #36210 Which commit carries a change is session state, not something a reader needs in order to rebuild the track. It belongs in the per-ticket notes under /hive/groups/browser/redmineNotes/36210/claude/, where it now is. diff --git src/hg/makeDb/doc/hg38/fiberSeq.txt src/hg/makeDb/doc/hg38/fiberSeq.txt index e9ba3dc578c..4a2bc77533c 100644 --- src/hg/makeDb/doc/hg38/fiberSeq.txt +++ src/hg/makeDb/doc/hg38/fiberSeq.txt @@ -1,370 +1,361 @@ # 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 4.1 million bases (PM00008, # the near-haploid Hap1 line, which has little phaseable heterozygosity). # 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. -# -# Where to find this change in git: the commits are -# d07356d7640 the PM00001 reissue and the nucleosome track -# 75e82896028 this hold-back -# The second one reads "uniprot otto: the miniprot cluster job needs absolute -# paths". Two sessions were committing in the same working copy at the same -# time, and the other one picked up these four files out of the shared index -# before its own commit. It was pushed that way, so the message stands; go by -# the file list, not the subject, when reading back through the log here. ~/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: 467 MB over the 41 files before rounding, 313 MB after. 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 3 values, all count > 1. Derived from the lab's free-text # cell type, and worth replacing when they give us real HPRC # metadata (donor sex, population) that would facet properly. # # 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. # # 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.