INFERNO — Containerized & Performance-Optimized

Reproducible, version-pinned Docker images that run the INFERNO pipeline — INFERring the molecular mechanisms of NOncoding genetic variants (Amlie-Wolf et al., Nucleic Acids Research, 2018) — for annotating noncoding GWAS variants and serving its website, with the historical command-injection vulnerabilities fixed.

This repository (a) faithfully freezes the published pipeline in a reproducible container (exact toolchain + data layout, validated bit-for-bit against the reference outputs) and (b) adds a substantial performance rebuild of the advanced analyses (GTEx colocalization, lncRNA correlation) and the annotation phase — every optimization gated to reproduce the published results to the byte. See Performance below.

Images

Image	Built on	Purpose
inferno-pipeline:dev	ubuntu:18.04 (pinned digest)	the reproducibility artifact — runs the core pipeline
inferno-website:dev	FROM inferno-pipeline:dev	Apache/PHP site + safe job submission
inferno-advanced:dev	FROM inferno-pipeline:dev	advanced steps — GTEx COLOC → lncRNA → pathway; see below

Annotation data is not baked into the image; it is mounted at /data at runtime.

Pinned toolchain

Component	Version	Source
Base OS	Ubuntu 18.04	digest sha256:152dc042452c496007f07ca9127571cb9c29697f42acbfad72324b2bb2e43c98
Python	2.7.17	bionic apt
R	3.6.3	CRAN bionic-cran35 apt repo (key E298A3A825C0D65DFD57CBB651716619E084DAB9)
R packages	—	date-pinned Posit PPM snapshot 2020-04-23 (just before R 4.0)
bedtools	2.25.0	exact production binary, vendored (pipeline/vendor/bedtools)
plink	1.90b6.24	vendored (pipeline/vendor/plink)
INFERNO code	commit 2e61f2948bb62cb17a876758f6f815ce2a3ce47e (branch active-website)	vendored git archive + one prod patch

Core R packages: data.table ggplot2 cowplot gplots gtools plyr psych RColorBrewer reshape reshape2 scales bitops caTools.

1. Stage the reference data (once)

The pipeline reads its reference data from /data at runtime — it is not baked into the images. Two ways to obtain it:

• Base annotation data (public): download the processed INFERNO annotations bundle, full_INFERNO_annotations.tar.gz, linked from the INFERNO web server, and extract it (~55 GB).
• From a host you control: scripts/fetch-data.sh / scripts/fetch-ibd-oracle.sh rsync the config-referenced directories from a remote you set via REMOTE=user@host. The advanced GTEx v6p reference (~110 GB total) must be sourced from GTEx separately.

REMOTE=user@your-data-host bash scripts/fetch-data.sh   # rsync into /mnt/ssd1/inferno/data
bash scripts/verify-data.sh                             # key files present

Expected layout (under the data root): HOMER_motifs, GTEx, 1000_genomes, roadmap, targetscan, factorbook, FANTOM5, refgenomes, INFERNO/input (rsID DB).

2. Build

docker build -t inferno-pipeline:dev pipeline      # base/reproducibility image
docker build -t inferno-website:dev  website       # FROM inferno-pipeline:dev

3. Run the pipeline (CLI)

bash scripts/run-example.sh        # runs the default SCZ2_128 example, output in _analysis/example/out
bash scripts/verify-config-paths.sh # (optional) assert every rendered-config path resolves

Direct invocation:

docker run --rm --user "$(id -u):$(id -g)" -e HOME=/tmp \
  -v /mnt/ssd1/inferno/data:/data:ro \
  -v /path/to/input.txt:/work/input.txt:ro -v /path/to/out:/work/out \
  -e KG_POP=EUR -e LD_THRESH=0.7 -e LD_AREA=500000 \
  inferno-pipeline:dev /work/input.txt /work/out MYPREFIX

The entrypoint renders the config from INFERNO_DATA/KG_POP/LD_THRESH/LD_AREA and runs INFERNO.py ... --cluster_system shell (sequential; no LSF/bsub).

4. Validate against the reference

bash scripts/validate.sh           # compares our run's metric lines to the oracle's

See Validation notes below.

5. Run the website

docker run --rm -d --name inferno-web -p 8080:80 \
  -v /mnt/ssd1/inferno/data:/data:ro inferno-website:dev
# http://localhost:8080/   (informational site + job-submission form)
bash scripts/test-website-submit.sh         # injection + form tests (fast)
RUN_FULL=1 bash scripts/test-website-submit.sh   # also drives one real submission (slow)

Validation notes

The oracle (example_output_scz2_128_top_variants.tar.gz) was produced with the reference pipeline using the same parameters (EUR / R² 0.7 / 500 kb). Validation compares the canonical metric lines (the LD-expansion SNP count and per-annotation overlap counts), which are path/run-id free.

Oracle summary-file caveat: the oracle's on-disk *_annotation_overlap_summary.txt is truncated to one line. On the reference system that file is written onto a goofys S3 mount with tee -a, and goofys does not support append (Operation not supported), so only the first metric persisted there. The complete reference set survives in the oracle's *_full_pipeline_log.txt; scripts/validate.sh reads from there. Our container writes to ext4, so our summary file is complete.

Security fixes (website)

The three historical command-injection vectors are fixed in website/html/analyze_top_snps.php: 1kg_pop (sed) → whitelist {EUR,AFR,AMR,SAS,EAS}; rsIDs (bash) → ^rs[0-9]+$ validation + escapeshellarg; file_upload (arbitrary copy) → default_inputs/ whitelist + realpath containment. outprefix and uploaded filenames are charset-constrained; all shell arguments are escapeshellarg'd; the pipeline runs via the containerized entrypoint (no sed-built config).

Advanced steps

inferno-advanced:dev (FROM inferno-pipeline:dev) adds INFERNO's GTEx colocalization chain — COLOC → lncRNA correlation → pathway — gated behind --run_* flags that require GWAS summary statistics. The base and inferno-website:dev images are never rebuilt.

Adds
R packages	coloc 3.2.1, WebGestaltR 0.4.3 + Bioconductor snpStats (same pinned PPM 2020-04-23)
Reference data	~110 GB GTEx v6p all-associations + RNA-seq under /data/advanced
Entrypoint	entrypoint-advanced.sh — env-driven config, --cluster_system shell

docker build -t inferno-advanced:dev -f pipeline/Dockerfile.advanced pipeline
bash scripts/fetch-advanced-data.sh                 # extract GTEx coloc/RNA-seq subtrees -> /data/advanced
bash scripts/fetch-ibd-oracle.sh                    # stage IBD inputs + published oracle
COLOC_CORES=8 bash scripts/validate-advanced.sh     # run COLOC+lncRNA, diff vs oracle

Validation. Both gates reproduce the published IBD oracle exactly — COLOC 202 colocalized genes (primary), lncRNA 4750 target genes (secondary) — using float-free hit-set comparison.

Performance (Tier-0 shipped). COLOC was I/O-bound — the read_eqtl pass re-decompressed each per-chr·tissue GTEx file once per region (~77% of wall-clock). The Tier-0 refactor reads each chr·tissue file once (data.table::fread, loop inversion) and parallelizes over (chromosome, tissue) tasks (COLOC_CORES, parallel::mclapply), dissolving the old chr6/MHC serial long-pole. Result on the 60-region IBD run (24 cores): COLOC phase ~30 min, down from ~9.3 h — 18.7× faster, reproducing the same 202 colocalized genes bit-exact, at ~3.2 GB/fork peak RSS (below the old 49.5 GB spike). The COLOC refactor is gated bit-exact by profiling/equiv-check.sh (per-gene coloc-summary diff) and profiling/find_genes_unit.R (the in-memory find_genes reproduces the awk early-exit exactly).

COLOC summary aggregation (Tier-0.5 shipped). After both compute steps were optimized, the post-COLOC "Reading in all summaries" phase — which rbind-grows one data.frame over 85,609 per-gene summary files — became the single largest phase (~74 min, O(n²): each iteration copies the growing 85k-row frame). Fix: collect the rows into a pre-sized list and do.call(rbind, …) once (byte-identical by construction; serial reads kept for Tier-0.5). Result: 4,444 s → 588 s (7.6×), 202 + 4750 bit-exact, gated by the (extended) profiling/equiv-check.sh. Full advanced pipeline ~3.4 h → ~2.5 h.

COLOC summary reads (Tier-0.6 shipped). Tier-0.6 then parallelized that residual read (for→mclapply over COLOC_CORES, mc.preschedule=TRUE; byte-identical by mclapply's order-preserving contract — read.table + do.call(rbind) kept verbatim, NULL/try-error fail-loud guard added). The read dropped 587.9 s → 353.1 s (1.7×) at 24 workers, 202 + 4750 still bit-exact, gated by profiling/summ-parallel-unit.R (invariants) + the existing equiv-check.sh aggregated-output diff. Honest finding: this is far below the profiler's ~12× projection because the 85,609-tiny-file read is metadata/IOPS-bound, not CPU-bound — per-file cost rose with worker count (~2.5 ms/file at 4 workers on the scoped gate vs ~4.1 ms/file at 24 on the full run: saturated random-read IOPS), so more cores don't help. The real remaining lever for this phase is reducing the number of files (consolidate the per-gene summaries into fewer, larger files), a separate larger change — not more parallelism.

COLOC summary consolidation (Variant A shipped). Acting on that lever: the per-gene *_summary.txt files were a pure write-then-reread-once round-trip inside one R process tree (the compute worker holds the coloc.abf summary in memory when it writes the file; the only reader was the aggregation). Variant A has each (chr,tissue) compute task return its summary rows in memory (do.call(rbind, …) per task), drops the L433 file writes entirely, and replaces the read loop with do.call(rbind, .coloc_results) + an order(.sortkey, method="radix") on each row's would-be on-disk path — which reproduces list.files(recursive=TRUE) order exactly (the image is C/POSIX, so list.files byte-orders == radix), keeping _gtex_coloc_summaries.txt byte-identical. One subtlety the equiv-check caught: the old path re-read each value from its file, and that write.table→read.table round-trip re-rounds the double (write.table renders ~15 sig figs here, so the raw in-memory double and the re-read double differ on ~7 of 85,609 rows); Variant A reproduces the round-trip in memory (read.table(text=capture.output( write.table(.))), no disk) so values stay bit-exact. Result: the aggregation phase collapses 587.9 s → 0.359 s and 85,609 per-gene summary writes are eliminated per run (the *_full_results.txt files are untouched). 202 COLOC + 4750 lncRNA bit-exact, and the full-scale 85,609-row _gtex_coloc_summaries.txt is byte-identical to the pre-change output (a full-scale byte-diff was added because the order-insensitive gene-set oracle can't catch a row-ordering regression). Full advanced pipeline ~2h13m → ~2h04m (the attributable win is the eliminated read phase + file I/O; compute-step wall-clock is run-to-run noisy). Gated by profiling/summ-consolidation-unit.R (radix==list.files, field-derivation parse, NULL-drop rbind, in-memory==on-disk round-trip) + equiv-check.sh (aggregated byte-diff + radix-order + C/POSIX assert) + the full oracle. This is the correct lever for the phase — removing the I/O, not parallelizing it.

COLOC eQTL reads → tabix (Tier-1 shipped). The last COLOC lever was the read itself: each (chr,tissue) worker fread(zcat)'d the whole chromosome file (~15.2 M rows, ~4 GB RSS) just to serve a handful of ±500 kb tag windows. A data-plane benchmark (profiling/dataplane-bench-*) measured bgzip+tabix range queries at 119× sparse / 14× dense cold-read and 15–70× less RSS vs the whole-file read, in-image (tabix/bgzip already ship in inferno-advanced:dev — zero new dependency, base image unchanged). Tier-1 productionizes it. An order-preserving re-encode tool (scripts/reencode-gtex-tabix.sh) prepends bare-chr+pos parsed from variant_id and bgzip+tabix -s1 -b2 -e2 indexes each file without sorting (the source is already position-sorted), asserting cmp(strip 2 cols) == original byte-for-byte per file so within-position row order — load-bearing for find_genes' first-contiguous-run — is provably preserved. The only production change (read_eqtl_for_task) swaps the whole-file read for one multi-region tabix query over the chromosome's interval-merged ±500 kb tag windows (disjoint + ascending ⇒ file-order, duplicate-free rows), selected per file by COLOC_READ=auto|tabix|zcat (auto = tabix iff the .txt.tbx.gz sibling exists, else the verbatim zcat path) — so deployed behavior is byte-identical until a file is re-encoded and rollout is incremental. find_genes, the per-gene ±500 kb GWAS intersection, coloc.abf, and aggregation are untouched. Validated on the IBD oracle's chromosomes (chr1+chr6, all 44 tissues re-encoded, 0 fidelity failures): a new row-order read-equivalence gate (profiling/tabix-read-equiv.sh) confirms the tabix read == the zcat whole-file subset byte+order identical on 4.7 M rows (incl. the MHC); the scoped equiv-check.sh in COLOC_READ=tabix is byte-identical; and a full mixed oracle run (tabix on chr1+chr6, zcat elsewhere) reproduces 202 COLOC + 4750 lncRNA bit-exact with the full-scale 85,609-row _gtex_coloc_summaries.txt byte-identical (same SHA-256) to the pre-change golden. Also gated by profiling/tabix-windows-unit.R (interval-merge / window-math / path / COLOC_READ-selection invariants). Full corpus since re-encoded (all 1012 (chr,tissue) files / ~104 GB, 0 fidelity failures via scripts/reencode-gtex-tabix.sh <corpus root>, idempotent + parallel) and re-validated all-tabix end-to-end: a full oracle with COLOC_READ=auto resolving to tabix on every chromosome reproduces 202 + 4750 and the 85,609-row summary byte-identical (same SHA-256).

Honest performance footnote. The benchmarked 14–119× / 15–70× RSS numbers are per-read figures and do not translate to a headline pipeline win. Post-Tier-0 each file was already read once, so tabix removes a cost that was no longer dominant: at COLOC_CORES=12 the all-tabix COLOC phase is 2407 s and the full pipeline ~1h51m (vs the mixed run's ~2h08m — directional, totals are run-to-run noisy). COLOC is now compute-bound (coloc.abf + find_genes + the serial pre-pass), not read-bound. And per-fork peak RSS is unchanged: median 2.96 GB, peak 3.93 GB across all 836 forks (≈ the Tier-0 zcat ~3.2 GB), because per-fork VmHWM is floored by the parent footprint inherited via fork() copy-on-write (the loaded GWAS stats + all-chromosome chr_ctx

• R/data.table/coloc libraries) — not the per-task eqtl_dt tabix shrank (the near-flat 2.88→2.96 GB band across 836 tasks is the tell). So the real lever to raise COLOC_CORES is shrinking that parent footprint (e.g. load per-chr context inside the worker rather than fork-inheriting all of chr_ctx), not the read. Tier-1's banked value is correctness-preserving infra + cold-read robustness + eliminated per-read decompress CPU/latency, not a wall-clock or RSS headline at this core count.

COLOC full-results toggle (COLOC_FULLRESULTS=all|hits shipped). With COLOC now compute-bound, the next discrete lever was disk: the (chr,tissue) compute workers write 85,609 per-gene *_full_results.txt files — one per gene — even though ~99 % of them are non-colocalized pairs that nothing downstream ever reads (top_coloc_hits re-reads only entries with PP.H4 ≥ coloc_h4_thresh). A new COLOC_FULLRESULTS=all|hits env toggle (default all) lets low-resource deployments opt in to hits mode, which gates each write on the round-tripped H4 (.rt[6,2] — the exact value top_coloc_hits re-reads from disk, so the written set equals the downstream-read set with no missing-file risk). Default all is unchanged and byte-identical — validated by profiling/equiv-check.sh in all mode (aggregated + all 7,382 full_results byte-identical) and the full 202/4750 oracle. In hits mode, the run writes 688 full_results (vs 85,609 in all), shrinking gtex_coloc/ disk from 62 GB → 351 MB (~61.6 GB saved); the COLOC step runs 1,415 s → 982 s (1.44×) at 12 cores (the proportional win is larger on low-core machines, where the per-gene file writes were a bigger fraction of per-core wall-clock). The aggregated _gtex_coloc_summaries.txt is byte-identical in both modes (the summaries come from the Variant-A in-memory path, untouched). 202 COLOC + 4750 lncRNA reproduced exactly in hits mode. Zero new dependencies. Gated by a mode-aware profiling/equiv-check.sh (all = full byte-diff; hits = strict-subset check on full_results) + the full hits-mode oracle.

COLOC parent-RAM trim (COLOC_TRIM=auto|keep shipped). This realizes the parent-footprint lever predicted in the honest-performance footnote above (shrink what fork() COW-inherits, not the read). Profiling the pre-fork parent found ~2.43 GB of objects the canonical path never reads: the per-chromosome GTEx-ID→rsID match tables (chr_ctx[[*]]$rsid_match, 2,469 MB over 19 chrs, read only by the exists("locuszoom")-gated locuszoom branch — and locuszoom is rm()'d to "NONE" in the canonical config) plus ensembl_xref (24 MB, dead once gene_name_lookup is built). Every worker COW-inherits them, and because R's GC writes SEXP headers during its mark phase, each forked worker that GCs copies the shared pages private (fork+gc COW-breakage) — so the dead weight is the bulk of the "~3 GB per fork" floor, not a shared one-off. A new COLOC_TRIM=auto|keep env toggle (default auto) drops both before the fork when locuszoom is off: it skips reading/holding the rsID tables (stores NULL in chr_ctx), rm()s ensembl_xref, and gc()s before mclapply; keep reproduces the old footprint (for the A/B + equiv gate) and locuszoom-on runs are unaffected (the read+hold path still runs). Byte-exact by construction (nothing reads the dropped objects on the canonical path) — validated by profiling/equiv-check.sh with COLOC_TRIM forwarded: both auto and keep reproduce the pre-change golden byte-identically (aggregated summaries + all 7,382 full_results), i.e. auto ≡ keep ≡ main, plus the full 202 COLOC + 4750 lncRNA oracle on default auto. Result: per-fork peak RSS (COLOC_PROFILE VmHWM, direct 19-chr COLOC run, 4 cores) drops from a tight ~3.2–3.5 GB (keep: min 3182 / median 3332 / max 3467 MB) to ~1.0 GB typical (auto: min 989 / median 1003 MB) — a ~2.3 GB-per-fork reduction, matching the directly-measured 2,493 MB parent dead-weight. This is a RAM/headroom lever (lets the pipeline fit / avoid OOM on a low-RAM box — where RAM, not cores, binds), core-count-independent, zero new dependencies., findings profiling/coloc-parent-ram-trim-findings.md.

Plot-skip toggle (INFERNO_PLOTS=on|off shipped, Phase A). A profiling pass (profiling/plotprof-*, which times the graphics-device lifecycle via R_PROFILE_USER) found that figure generation is ~24–26 % of the pipeline — 874 figures / ~615–685 s wall, the bulk of the previously "unaccounted" phase-breakdown block. Of that, ~315 s is serial, single-threaded plotting in the COLOC (47 figures, 145 s) and lncRNA (357 figures, 170 s) scripts — core-count-independent, so an even larger share on a low-core box. A new INFERNO_PLOTS=on|off env toggle (default on) skips it: every figure goes through a make_graphic(name,…); <plot>; dev.off() span, and a deterministic transform (profiling/gate-plots.awk) wraps each of the 47 COLOC + 24 lncRNA spans in if(MAKE_PLOTS){ … } so off skips the build entirely. Gating only device-drawing spans is byte-safe by construction (no result-table side effect lives inside a span — verified across all 71) and on (the default) runs byte-identically (the braces are transparent). No INFERNO.py change — INFERNO_PLOTS is set via docker -e and inherited by every child Rscript. Validated: COLOC equiv-check.sh reproduces the pre-change golden byte-identically in both on and off (aggregated summaries + all 7,382 full_results); off writes 0 COLOC/lncRNA figures with the 202 COLOC + 4750 lncRNA oracle still exact, while on restores the full 404-figure set — both VALIDATION PASSED. Measured win: full pipeline 2563 s → 2190 s (−372.86 s, both warm). Core-count-independent, zero new dependencies. The already-parallelized annotation plotting (~300–370 s) is Phase B; keeping all figures while speeding them on a multi-core box is the separate SP2 (lncRNA plot parallelization)., findings profiling/plot-skip-toggle-findings.md.

lncRNA correlation (Tier-0 shipped). With COLOC down to ~30 min, the lncRNA step became the new long-pole. Profiling (profiling/run-lncrna-profile.sh) found §4 — the correlation loop — was ~95% of it, and that psych::corr.test spent most of its time computing p-values the result-defining 0.5/0.5 correlation threshold never reads. Tier-0: replace corr.test with cor + a vectorized p-value reconstruction that reproduces corr.test's output byte-for-byte (matched to its exact (r·√(n-2))/√(1-r²) / -2·expm1(pt(…,log.p)) arithmetic, incl. the use="pairwise" cor mode); hoist the lncRNA-independent gene_resids out of the per-lncRNA loop; and parallelize the per-chromosome loop (parallel::mclapply, reuses COLOC_CORES, deterministic chromosome-order combine + fail-loud worker guard). Result on the 60-region IBD run (24 cores): §4 17,303 s → 1,480 s — 11.7× faster (lncRNA step ~5 h → ~40 min), reproducing the same 4750 target genes bit-exact (and the per-lncRNA tables byte-identical), at ~3.3 GB/fork peak. Gated bit-exact by profiling/lncrna-equiv-check.sh (per-lncRNA table diff). See profiling/ for the instrumentation and the remaining roadmap (full-corpus tabix re-encode → Parquet/Tier-2 → S3/Batch region-grain for COLOC; lncRNA §1–3 setup + §5 if they prove material).

Baseline annotation pipeline parallelism (Increment 1 shipped). With COLOC + lncRNA optimized, the pre-COLOC annotation phase became the next serial chunk: src/expand_and_annotate_snps.py (LD expansion + per-source overlap annotation, ~592 s, single-threaded Python 2.7) and analysis_scripts/Rscript_run_full_analysis.R (the "Analyzing X" summaries + plots, ~808 s). Increment 1 is pure parallelization, byte-exact by construction (no algorithm change; each unit writes its own output): a per-source multiprocessing fan-out in the Python annotator (expand_ld_snps stays serial, then ~9 independent sources run concurrently so HOMER + GTEx stop running back-to-back) + an mclapply wave over the independent R summary blocks (1–11 concurrent, block 12 serial), both bounded by COLOC_CORES, fail-loud, zero new dependencies (py2 stdlib multiprocessing + base-R parallel; no pandas/numpy in the py2.7 image). The existing 202/4750 oracle is downstream of and blind to the annotation outputs, so a new gate profiling/annot-equiv-check.sh byte-diffs the annotation text outputs against a golden — with a baseline byte-reproducibility pre-check that (correctly) flagged the PLINK .log / command_run_date timestamp files and excluded them. The parallelized files override the prebuilt base image's serial copies via a Dockerfile.advanced COPY (we don't rebuild the base, so the deployed base image stays unchanged; serial-by-default keeps either image safe). Result (60-region IBD, COLOC_CORES=12): annotation phase 1400.7 s → 1136.0 s (−19%) — the R summary 808 → 629 s (−22%) is the clean win (CPU-bound), while the Python annotator (592 → 507 s, −14%) is I/O-bound: its dominant cost is the plain-text 13 GB HOMER + 13 GB GTEx scans, which contend for disk when fanned out, so parallelism yields diminishing returns. 202 + 4750 reproduced exactly and annotation text outputs byte-identical. The honest read: the I/O-bound Python annotator is precisely the signal that the structural lever is Increment 2 — tabix range queries (read kilobytes, not 13 GB) — not more parallelism.

Annotation eQTL/motif reads → tabix (Increment 2 shipped). Acting on that signal, the two 13 GB long poles — gtex_eqtl_overlap (44 GTEx .snpgenes) and homer_motif_overlap (the HOMER motif BED) — now answer their ~2,753 LD-SNP point/interval queries with bgzip+tabix range reads instead of scanning the whole corpus. A new in-image, idempotent, cmp-gated re-encode tool (scripts/reencode-annot-tabix.sh + -one.sh) builds the siblings: GTEx prepends bare-chr+pos parsed from the snp variant-id and indexes tabix -s1 -b2 -e2 -S 1 (placeholder header skipped), asserting cmp(strip 2 cols) == original; HOMER pastes the BED with its line-paired sequence file into one record and indexes tabix -s1 -b2 -e3 (1-based fully-closed) — deliberately not -p bed, whose 0-based half-open model silently drops every snp_pos == motif_start overlap (a latent, dataset-dependent byte-break a position-sampled gate can't reliably catch; the correct preset was proven equivalent in-image before any 13 GB re-encode) — asserting both the bed columns and the seq column cmp-equal to the originals. The read-swaps preserve each function's exact overlap semantics: GTEx the unchanged chr:pos:ref:alt exact-match (byte-identical after the existing final sort), HOMER the unchanged fully-closed motif_start <= snp_pos <= motif_end test with the sequence read from the co-located combined row (row[7]) so the PWM ref/alt/delta math and per-SNP corpus-order emission stay untouched. Selection is per-file ANNOT_READ=auto|tabix|scan (auto = tabix iff the .tbx.gz sibling exists, else the verbatim scan — so deployed behavior is byte-identical until a corpus is re-encoded); tabix + missing sibling → fail-loud, query error → no fallback. Validated in layers: per-file cmp (44 GTEx + HOMER, 0 fidelity failures); a new profiling/annot-tabix-read-equiv.sh (tabix == scan candidate rows, byte+order identical, GTEx exact + HOMER closed-interval); the authoritative profiling/annot-equiv-check.sh in ANNOT_READ=tabix (fail-loud, both sources forced through tabix) — every annotation text output byte-identical to the all-scan golden (gtex_eqtl_overlap 435 + homer_motif_overlap 144 rows); and the full IBD oracle 202 COLOC + 4750 lncRNA reproduced exactly. Result (COLOC_CORES=12): the Python annotation phase 506.7 s → 414.7 s (−18 %) as the two 13 GB scans collapse to kB queries (the residual is the plink LD-expansion + the still-scanned smaller sources). Honest read: a single run on a shared box, so the banked, proven value is the byte-exactness + the eliminated ~26 GB of sequential reads, not a full-pipeline headline (the COLOC/lncRNA/R deltas in this run are run-to-run noise; those phases don't consume the annotation outputs). Full-corpus re-encode (44 GTEx .snpgenes + the 5.5 GB HOMER sibling) is a one-time documented ops step via scripts/reencode-annot-tabix.sh; the remaining smaller scan-bound sources (Roadmap / TargetScan / factorbook / FANTOM5 / dashr) are deferred to a later increment.

R summary plot-render parallelization (Increment 3 shipped). With the annotator tabix'd, the R summary (Rscript_run_full_analysis.R, ~629 s in Increment 1) became the largest serial chunk. A profile-first pass — a standalone in-image harness (profiling/rsummary-standalone.sh) replaying the R step against a staged 60-region input under serial Rprof — proved the step is plotting-bound, not data-prep-bound (~63 % of self-time is ggplot/cairo heatmap rendering + dev.off() flushes; the candidate data-prep loop was 4.7 s / 0.4 %). That fired the design's honesty escape-hatch: no byte-exact data-prep change moves the wall-clock, so the increment pivoted to rendering the plots concurrently. The single biggest cost — the Roadmap per-state ChromHMM heatmap loop (analyze_roadmap_chromhmm_function.R:186, ~346 s rendering ~15 independent heatmaps serially) — was converted from a for loop to mclapply(mc.cores=COLOC_CORES) with a fail-loud check. Each iteration writes its own PNG from read-only globals (no shared state, no RNG) → embarrassingly parallel; default COLOC_CORES=1 ⇒ mclapply runs serially in-process ⇒ byte-identical to the loop (same default-serial contract as Increment 1; parallel already loaded, zero new deps). Verification unlock: the 15 per-state PNGs are cairo-deterministic and serial/parallel-invariant (no tIME chunk), so a plain cmp is an exact appearance gate — no pixel-diff tooling needed (none exists in-image). A new profiling/roadmap-heatmap-cmp.sh enforces it; annot-equiv-check.sh (which excludes *.png) is unaffected and stays green. Validation: ANNOT-EQUIV OK on the rebuilt image; the baked 60-region run text byte-identical + all 15 per-state heatmaps byte-identical (serial and nested-fork parallel paths); full IBD oracle 202 COLOC + 4750 lncRNA reproduced exactly. Result (COLOC_CORES=12, 3-run medians on a shared box): R summary 579 s → 342 s (−41 %, 1.69×) — the recovered ~237 s is the redundant serial cairo rendering now spread across cores. Honest Amdahl note: Wave-2 (analyze_fantom5_eqtl_chromHMM_overlap, also plotting-bound but PDF output → not cmp-verifiable under the zero-deps constraint) is untouched and now dominates the R step; it and the smaller scan-bound annotation sources are deferred. (PIVOT section),

1000G LD-expansion reads → tabix (Increment 4 shipped). The last untouched scan was Phase 1, LD expansion (expand_and_annotate_snps.py::expand_ld_snps): for each tag SNP it collects the variants within ±LD_AREA (±500 kb) from the 1000G reference by scanning the entire per-chromosome VCF in pure-Python gzip.open(), then hands a small per-SNP VCF to plink for the r² math. A profile-first gate (env-gated KG_PROFILE timers, 60-region, COLOC_CORES=12) proved the lever decisively — and as the inverse of Increment 3: gzip-scan 330.5 s vs plink 4.9 s, i.e. the VCF scan is 98.5 % of LD-expansion compute and plink is negligible. The fix is the same tabix range-query lever as Increment 2, on a new corpus: tabix -h <file> <merged windows> streams the header + only each SNP's ±500 kb window instead of the whole chromosome. Approach-A source-swap — the per-chromosome processing loop body is verbatim; only the line source changes (in both chromosome-flush blocks: the in-loop "previous chr" and the after-loop "final chr"). Byte-exact by construction: the merged-window stream is a superset of every SNP window and nothing outside a window touches state, and tabix -h copies the header verbatim. Zero new deps: tabix is already in-image and the 1000G VCFs are already BGZF, so this is index-only (tabix -p vcf → .tbi, no recompression) via scripts/index-1000g-tabix.sh; no bcftools, base image unchanged. Selection is KG_READ=auto|tabix|scan (mirrors ANNOT_READ/COLOC_READ): auto = tabix iff a .tbi exists, else the verbatim scan — so deployed behavior is byte-identical until a corpus is indexed. Two gotchas pinned: ±500 kb windows of nearby tag SNPs overlap, so windows are merged disjoint (overlapping regions make tabix emit duplicate records); and the VCF contigs are bare (1, not chr1). Validation: annot-equiv-check.sh already byte-gates the whole ld_expansion/ tree (neighbor VCFs + plink .ld + final _ld_cutoff_snps_within_*.txt) — ANNOT-EQUIV OK at 4-region with KG_READ=tabix forced, plus a 60-region ld_expansion/ byte-diff identical under both tabix and scan (excluding plink's timestamped *.log). Because LD expansion is upstream of everything and its output is now byte-identical, the COLOC/lncRNA results are identical by construction, so the 2 h 202+4750 oracle is redundant here and demoted to optional. Result (COLOC_CORES=12): the per-chromosome VCF read collapses from ~330 s (clean baseline; up to ~900 s under box contention) to ~12 s of windowed tabix reads — a ~28–77× cut that no longer scales with corpus size; plink (~5 s) is untouched. Honest single-run-on-shared-box caveat (the scan baseline is I/O-bound and load-sensitive). The other four 1000G populations index via scripts/index-1000g-tabix.sh all (one-time ops); plink (4.9 s) is left as-is.

lncRNA correlation §4 rank-hoist + lm vectorization (shipped). With COLOC and the annotation/LD scans optimized, the lncRNA correlation step (Phase 5) was the clean-cold #2 lever — §4 ("Computing all expression correlations") ≈ 1,437 s, ~25 % of the 5,790 s pipeline. A full-scale per-chromosome op-timer (profiling/run-lncrna-prof2.sh + lncRNA_gtex_correlation.prof2.R, production parallelism, real coloc summaries, COLOC_CORES=12) localized it precisely: spearman ranking is ~77 % of §4 CPU, almost all redundant — R computes cor(x, M, method="spearman") by re-ranking the entire gene matrix M on every one of the 40 colocalized lncRNAs (37× more expensive than the identical-shape pearson call, 6,986 s vs 189 s) — and the tissue-specific per-gene lm is another 13.6 %; read.table is only 3.7 % (so fread is not worth the byte-risk; the escape-hatch fired). Two arithmetic-equivalent rewrites, confined to §4: (1) rank-hoist — rank each gene matrix once per chromosome / tissue-class (gene_rank_mat, gene_resids_rank, tiss_rank_mat) and reuse it across all lncRNAs via a new cor_with_p_cached(x, rankM, use=…) helper, since spearman ≡ pearson-on-ranks through R's own C_cor path (use="all.obs" mirrors the use="pairwise" spearman of cor_with_p; use="everything" the default-use partial cor); (2) lm vectorization — replace the per-gene formula-lm with one lm.fit on the shared model.matrix(~ PCs) design (one QR solves all gene columns; lm is lm.fit on that design). No production toggle (unlike the *_READ increments): old and new compute identical bytes, so there is nothing to switch between — the byte-gates are the safety net. Byte-exact by construction, re-proven empirically through three gates: an in-image mechanism unit test (profiling/lncrna-rank-hoist-unit.R, identical() on real chr21 data — including the production-faithful vector shape of this_lncrna_resids, which residuals(lm(n×1 matrix ~ tissues)) drops to a plain vector); the scoped chr21+22 byte-diff (profiling/lncrna-equiv-check.sh → EQUIVALENCE OK, 31 per-lncRNA tables incl. tissue-specific); and the full IBD oracle 202 COLOC + 4750 lncRNA reproduced exactly (both diff files empty). Result (COLOC_CORES=12): §4 1,437 s → 195 s — 7.4× faster (over-delivering on the ~3.9× design estimate, because the two levers together attack ~90 % of §4 CPU); the correlation loop is no longer the long pole — the lncRNA-expression read (~419 s, untouched) is now the step's dominant cost. Parallelism (mclapply over 24 chromosomes, mc.preschedule=FALSE) is unchanged; the chr1 makespan tail remains (finer-grain load-balance is the next lever, only material at >12 cores). The full cold pipeline wall came in at 4,874 s this run (≈ −15 % vs the 5,790 s clean baseline; pipeline totals carry run-to-run / shared-box noise, so the directly-measured, byte-gated §4 7.4× is the honest headline). Profile gate profiling/lncrna-rank-hoist-findings.md,

lncRNA correlation §3 row-prefilter read (shipped). With §4 optimized, the dominant cost of the lncRNA step shifted to §3 ("pull out colocalized lncRNAs' expression vectors") — previously measured but untouched. A full-scale per-op serial timer (profiling/run-lncrna-prof3.sh + lncRNA_gtex_correlation.prof3.R, real coloc summaries, in-image) confirmed the profiling finding: read.table is 83 % of §3 — 14 full per-chromosome expression matrices (chr1 = 462 MB, 8,555 sample columns) read serially just to extract the ~40 colocalized lncRNA rows. The fix reads only the lncRNA gene rows via an awk row-prefilter piped into the same read.table parser (LNCRNA_READ=auto|prefilter|full, default auto ≡ prefilter): awk keeps the header and any row whose Description field is in the chr's target set, so !duplicated picks the same first occurrence and drop=FALSE on single-column subsets preserves the matrix layout — byte-identical by construction. Byte-exact by construction, re-proven through three gates: an in-image unit test (profiling/lncrna-section3-read-unit.R, identical() on multi-row, dedup, and single-column cases); the scoped chr21+22 byte-diff (profiling/lncrna-equiv-check.sh → EQUIVALENCE OK in BOTH prefilter and full modes, the full path being main + the no-op drop=FALSE); and the full IBD oracle 202 COLOC + 4750 lncRNA reproduced exactly. Result (COLOC_CORES=12): §3 402.9 s → 61.9 s — ~6.5× faster, ~341 s off the pipeline (~6 % of the 5,790 s cold baseline). Core-count independent — it removes work rather than spreading it, so the gain holds on the low-resource deployment box. Honest perf: the realized ~6.5× is below the ~15× design estimate because §3 is now disk-scan-bound, not parse-bound — the awk prefilter still streams the full ~2.5 GB of per-chromosome files sequentially, so what was eliminated is the read.table parse (the bulk of the old ~348 s read), while the sequential scan (~62 s) remains; parallelizing the awk scans over chromosomes, or sharing §4's parallel read, is the natural next lever. Profile gate profiling/lncrna-section3-read-findings.md,

COLOC worker overhead hoist — Levers 1 + 2 (shipped). With lncRNA retired as the #2 lever, attention returned to COLOC (Phase 4), now the dominant pipeline cost. A full-scale per-gene op-timer on the COLOC worker loop (profiling/run-coloc-profile.sh, real 60-region input, COLOC_CORES=12) pinned the hot spots precisely: coloc.abf (the CRAN dependency) is only 12.6 % of worker CPU — the dominant reducible costs are the per-gene gene_name grep (~32 %, a tissue-invariant regex over the 204,940-row Ensembl xref, repeated for every gene × tissue pair) and per-gene match + coerce prep (~21 %, a full re-scan + column coercion of the tag's eQTL table, repeated for every gene). Two arithmetic-equivalent rewrites: (Lever 1) build a gene_name_lookup named-vector once before the parallel workers start (COW-shared via mclapply fork) and look up each gene by ID in O(1), replacing the repeated regex; (Lever 2) coerce the numeric eQTL columns once per tag and split() the coerced table by gene once, so each worker receives its pre-sliced rows directly, replacing the per-gene full-table re-scan and per-gene coercion. No production toggle (no *_READ environment variable): old and new compute identical bytes, so there is nothing to switch between — the byte-gates are the safety net. Byte-exact by construction, re-proven through three gates: an in-image unit test (profiling/coloc-overhead-hoist-unit.R, identical() on Lever-1 lookup and Lever-2 split); the scoped 4-region byte-diff (profiling/equiv-check.sh → EQUIVALENCE OK, all COLOC + lncRNA outputs); and the full IBD oracle 202 COLOC + 4750 lncRNA reproduced exactly (both diff files empty). Result (COLOC_CORES=12): COLOC 2,416 s → 1,415 s — 1.71× faster (the two levers together retire ~53 % of the worker's non-coloc.abf overhead). Zero new dependencies; base image unchanged.

Pathway step (optional, runtime-deferred). --run_pathway_analysis (WebGestaltR GO/KEGG over-representation) is wired and its script is corrected (a 94-line verbatim-duplicate block that broke parsing was removed), but it cannot currently run against the live WebGestalt API: the pinned WebGestaltR 0.4.3 fails gene ID mapping ("Organism is missing") against the 2026 service (listOrganism() works; ORA's idMappingGene does not). Pathway is not part of the canonical IBD run and has no oracle. Making it runnable needs a newer WebGestaltR (risks the R 3.6.3 pin that COLOC/lncRNA reproducibility depends on) or an offline GMT-database mode.

Still deferred

• MetaXcan and LDSC — external Python 2 tools; planned as modern Py3 / uv forks in their own layers.
• Pathway runtime — see the WebGestaltR-vs-live-API note above.
• Large website download artifacts (e.g. scz2.snp.results.1kg_annotations.txt.gz) are served separately, not baked into the image.

License & attribution

This repository is released under the MIT License (see LICENSE). It packages and extends INFERNO, which is itself MIT-licensed (© University of Pennsylvania; upstream notice retained at pipeline/vendor/inferno-code/LICENSE.md). Vendored third-party components — bedtools, plink, and the PHP libraries under website/html/vendor/ — retain their own licenses.

Upstream INFERNO: https://bitbucket.org/wanglab-upenn/INFERNO

Citation

If you use INFERNO, please cite:

Amlie-Wolf A, Tang M, Mlynarski EE, et al. INFERNO: inferring the molecular mechanisms of noncoding genetic variants. Nucleic Acids Research 46(17):8740–8753, 2018.