d07356d76401059fbe6c2d0890b1490546a62a4e max Mon Sep 14 02:54:50 2026 -0700 hg38 Fiber-seq: reissued GM12878 data and a nucleosome density track, refs #36210 The lab reprocessed GM12878 (PM00001) and replaced the files in place under the same hash directory. Ten of its twelve files changed; a size sweep over all 41 samples confirmed no other sample is affected. This fixes the two placeholder haplotype accessibility bigWigs that covered a single base, so that overlay now draws real data for the sample that comes up by default. Its peak calls changed substantially as well, 429,883 source peaks before and 196,742 now, which is noted on the description page since figures made from the first version of the track will not reproduce for GM12878. The downloader now fetches into <file>.part and moves it into place when complete. It used curl -C - straight onto the final file, which is right for an interrupted transfer and silently corrupting when the server has replaced the file: it would have appended the tail of the new 5.2 GB hap1 file to the 512-byte stub, and the size check afterwards would have passed. It also takes an optional list of accessions now, to refresh one sample without walking all 41. Nucleosome density (all.nucleosome.coverage.bw) was sent separately and is not in the lab's own hub. It is on the server for all 41 samples and is added as a seventh data type in the compendium. Unlike every other wiggle here it is a read depth rather than a percentage, so it cannot take fixed viewLimits: the genome-wide mean runs from 25 to 142 across samples with sequencing depth and single loci reach 1.7e5. It is drawn with autoScale, which the description page explains, and reads as the complement of the accessibility signal. diff --git src/hg/makeDb/doc/hg38/fiberSeq.txt src/hg/makeDb/doc/hg38/fiberSeq.txt index c22ae9d3be3..ce6d30291ec 100644 --- src/hg/makeDb/doc/hg38/fiberSeq.txt +++ src/hg/makeDb/doc/hg38/fiberSeq.txt @@ -45,58 +45,108 @@ 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 # --------------------------------------------------------------------------- -# 11 files per sample, 451 files, 397 GB: +# 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. +# 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 -# This reports two placeholder files at the source: +# 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 are valid bigWigs that cover a single base with a value of zero, so -# GM12878's haplotype accessibility overlay draws nothing. Note that they are +# 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). -# Every other sample has real haplotype files, and PM00001's CpG haplotype -# files are fine. Reported to the lab; noted on fiberSeqCompendium.html. -# Nothing else was dropped: all 451 files downloaded, and all 449 non-empty -# ones parse as bigWig or bigBed. +# 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 +# --------------------------------------------------------------------------- + +# 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 + +# 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. Drawn with +# autoScale on, so the axis is comparable within a sample and not between them. + +~/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. @@ -110,82 +160,89 @@ # 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,487,043 -# peaks in track 9,486,622 -# dropped chrEBV 421 in 20 of the 41 samples, 2 to 166 each +# 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 166 peaks, which is 0.04 percent of its 429,883. +# 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 9486622 source 9487043 chrEBV 421 mismatched 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 585 stanzas stay consistent: +# metadata and color files, so the 625 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 +# fiberSeqCompendium faceted composite, 41 samples, seven data types: +# acc, peaks, hap, nuc, cpg, cpgHap, cpgDiff # # 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 +# data type as the inner one, so a single composite keeps a sample's seven 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