Oncofusion

start / graph / concept-q154-a-genome-wide-gene-body-matrix-costs-one-fragment-bed-of-peak-disk

Explainerstrong confidence 4 of 5: Replicated in peer-reviewed work by more than one group.Q154-gene-scale-occupancy-null

A genome-wide gene-body occupancy matrix over 133 libraries costs one fragment BED of peak disk, not 133 - stream, bin in flight, delete

also called bin-and-delete over HTTPS · peak disk is not total transfer

Working out how much of a protein sits on every gene in the genome, for 133 sequencing runs, sounds like it needs a lot of disk space. It does not: you can download each file, count as it arrives, throw it away, and keep only a tiny table - so the most disk you ever use at once is one file.

The same thing, for a scientist

A 19,402 gene-body x 133 library occupancy matrix for GSE252378 was built with **peak disk equal to one fragment BED** (~100 MB) plus the output TSVs (~190 kB each, ~25 MB total), on a machine with 9.3 GB free on C: and 12.5 GB on D: and **zero cached bigWigs** [E1]. Each `GSM*_hg38.bed.gz` is streamed from the GEO FTP over HTTPS into a temporary file, piped through `zcat` into a single awk pass that assigns each fragment midpoint to a 500 bp bin and increments the intervals that bin carries, and deleted before the next download starts. **The whole 133-library pass took about 45 minutes at a median of roughly 19 s per library**, dominated by transfer, not by compute [E1]. The bin map holds 1,223 Mb of interval span in ~2.45 M bins; **39,192 bins carry more than one interval**, and the map stores a list at those bins so nested and overlapping intervals - here, ten nested length-matched control windows per negative-control locus - are counted exactly rather than first-come-first-served [E1].

Go deeper
Why this is a concept and not a note

Two beats in this project have deferred a genome-wide re-analysis on the grounds that the input data would not fit. The relevant quantity is peak disk, not total transfer. 133 fragment BEDs are perhaps 15 GB in total and do not fit; one at a time is 100 MB and fits on a full disk with room to spare. The same pattern applies to bigWigs, BAMs and any per-sample archive: if the reduction is a fixed-size summary, the archive never needs to land.

The q144 beat wrote the pattern first, for promoters; this beat swapped the interval file for gene bodies and reused everything else 1. That is the second point: an extractor whose only sample-specific input is an interval table is reusable at any scale, and building the next matrix costs an interval file rather than a pipeline.

The parts that are easy to get wrong
  • A bin can hold more than one interval. The promoter version of this awk kept the first interval to claim a bin and counted collisions; that is safe for disjoint 2 kb windows and wrong for nested windows. This version stores a list at multi-interval bins - 39,192 of them - and checks containment against every one.
  • Verify a failed download rather than a low count. The wrapper treats a curl non-zero exit or a total fragment count below 100,000 as a failure, deletes the partial output, and reports it. A partial matrix that looks like a complete one is the failure mode this project has paid for repeatedly.
  • A background job's notification exit code is the wrapper's, not the command's. The run was checked by opening the log, counting the output files and reading the tail, not by trusting a completion notice.
  • Fragment midpoints are not coverage. A midpoint count loses signal at an interval's edges in proportion to fragment length over interval width, so it understates short intervals relative to long ones. Measured against q147's bigWig mean coverage over 21 named regions, the two agree at Spearman 0.9932 with a per-region median ratio of 1.038 (range 0.867-1.246) 1 - but the direction of the residual bias is known and must be declared, not discovered afterwards.
Facts and figures
scriptanalyses/q154-gene-scale-occupancy-null/extract.sh
statusrun
result summary133 libraries, peak disk one fragment BED plus ~25 MB of TSVs, ~45 minutes, Spearman 0.9932 against the bigWig route.
n intervals19503
total interval span mb1223
n multi interval bins39192
bin size bp500
inputs
[
  "analyses/q154-gene-scale-occupancy-null/intervals.tsv (19,503 intervals, 1,223 Mb)",
  "GSE252378 GSM*_hg38.bed.gz over the GEO FTP"
]
measurements
[
  {
    "quantity": "peak_disk",
    "subject": "building a 19,402 x 133 gene-body occupancy matrix from GEO fragment BEDs by streaming",
    "value": 100,
    "units": "MB (one fragment BED) plus ~25 MB of output TSVs",
    "n": 133,
    "universe": "GSE252378, 133 AutoCUT&RUN libraries, hg38 fragment BEDs on the GEO FTP",
    "null": 15000,
    "note": "the null is the approximate total transfer if every BED were landed; the whole scratch directory ended at 26 MB, and the free-space reading on C: is not a usable bound because other processes moved it"
  },
  {
    "quantity": "instrument_agreement_spearman",
    "subject": "fragment-midpoint density shares versus q147's bigWig mean-coverage shares",
    "value": 0.9932,
    "units": "Spearman rho",
    "n": 1134,
    "universe": "21 named regions x 54 KMT2A-N libraries; per-region median ratio 1.038, range 0.867-1.246",
    "null": null,
    "note": "the pooled rho is inflated by the regions spanning three orders of magnitude; the per-region ratio is the informative figure"
  }
]
Still open
  1. The same extractor with an enhancer or CpG-island interval table would give this project its first non-promoter, non-gene-body occupancy matrix for the same cost. Which interval set would most change the current conclusions?
  2. Fragment midpoints understate short intervals relative to long ones. What is the exact size of that bias as a function of interval width, measured by re-running one library with full-fragment overlap counting instead of midpoints?
  3. GSE159608 (CD34+ controls, hg19) was partially extracted by the q144 beat against an hg19 promoter table. Would a gene-body pass over the normal-HSPC libraries give the missing not-a-leukaemia baseline for the pass rate?
Every source behind this entry (1)
  1. E1strong confidence 4 of 5: Replicated in peer-reviewed work by more than one group.dataset

    A 19,402 gene-body plus 101 control-interval matrix over all 133 GSE252378 libraries was built by streaming each fragment BED from the GEO FTP, counting midpoints in a single awk pass with a 500 bp bin map holding lists at the 39,192 multi-interval bins, and deleting each stream; peak disk was one BED (~100 MB) plus ~25 MB of output, and the pass took about 45 minutes. Agreement with q147's bigWig-derived shares over 21 named regions is Spearman 0.9932 with a per-region median ratio of 1.038.

    Q154 streaming gene-body extractor and instrument gate · Iqbal Lab oncofusion knowledge graph, beat Q154 · this project 2026
    What this evidence cannot show

    The 100 MB peak-disk figure is the size of the one stream measured directly (GSM8000495, 100,545,363 bytes); other libraries in the series are of similar order but were not individually measured, and a larger one would raise peak disk proportionally. Free space on C: was 9.3 GB before the run and 8.6 GB after it, but the whole scratch directory is only 26 MB, so most of that 0.7 GB drop belongs to other processes on this machine and the free-space delta does NOT bound the peak. The peak is bounded by the one-file-at-a-time structure of the loop and by the measured size of one stream, not by the disk reading.

    Evidence that argues against this (1)
    • Streaming trades disk for network and for re-runnability. Every re-analysis at a new interval set costs the full transfer again - here about 15 GB - and a GEO outage or a rate limit stops the run where a local copy would not. On a machine with room, caching the archive once is the better trade. NCBI Gene Expression Omnibus 2024
Connected to