75e828960283291546d2c1a27845e2cf3823adcd max Mon Sep 14 05:29:02 2026 -0700 uniprot otto: the miniprot cluster job needs absolute paths GRCz12ab failed with the parasol job crashing four times, return 1, no output. The wrapper I wrote ran miniprot -t 16 --gff protToGenome/GRCz12ab/.../genome.fa fasta/7955.fa > $1 and a parasol job runs with its working directory set to the batch directory, not to the directory the pipeline runs in, so neither input existed from the job's point of view. The BLAST batch next door gets away with relative paths because it cds into its own workdir and its jobList is written relative to that; this batch directory sits a level deeper and its paths were relative to the otto root. Every path in the wrapper, the jobList command and the output check is now absolute. Verified on the cluster against the real 1.48 Gb zebrafish genome: successful batch, a 195 MB GFF with 93518 mRNA records. refs #38300 diff --git src/hg/makeDb/doc/hg38/fiberSeq.txt src/hg/makeDb/doc/hg38/fiberSeq.txt index ce6d30291ec..7abb25d631e 100644 --- src/hg/makeDb/doc/hg38/fiberSeq.txt +++ src/hg/makeDb/doc/hg38/fiberSeq.txt @@ -1,347 +1,360 @@ # 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 +# 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 -# 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. +# 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, which is the check that the counter renumbers correctly. ~/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 625 stanzas stay consistent: +# 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, seven data types: -# acc, peaks, hap, nuc, cpg, cpgHap, cpgDiff +# 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 seven subtracks +# 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.