Skip to content

Reproduction

Jump to bottom
Wenyu (Eddy) Huang edited this page Jun 30, 2026 · 3 revisions

Reproduction

There are three tiers of reproduction depending on what you want to do:

Tier What you do What you need Where to look
1 - Use the shipped tables Render any figure (Figs 2-8 + Tables S1, S2) from the analysis tables shipped with the repo; sweep thresholds, change the dataset, or write your own notebooks against the preprocessed CSVs. No profiler runs, no databases, no reads. Conda env only Quick reproduction
2 - Re-process from raw reads Trust the shipped reference databases, but re-run all six profilers - to validate the published results end-to-end, test custom parameters, or extend the analysis to new samples. Env + reference DBs + reads (Data Fetch) Full reproduction
3 - Rebuild / extend the reference Update the RefSeq snapshot, add new genomes, or audit the build end-to-end. Env + raw RefSeq genomes + taxdump + everything in Tier 2 Database Build

Most readers only need Tier 1. Tier 2 is for users who want to re-process the pipeline (validate, or modify parameters); Tier 3 is for researchers extending or rebuilding the unified reference.

Run everything from the bakeoff/ repository root. Notebook import paths and default flags assume this working directory.

What the repository provides

To make the quick path possible, the repository ships the pre-computed analysis tables - the output of the heavy pipeline stages - so you can jump straight to the figure notebooks:

Provided in the repo Path Feeds
Preprocessed per-(sample, tool, db) tables results/metadata/<ts>/analysis-prep/preprocessed/{detection,abundance}/ Figs 2-6
Ground-truth tables + read totals results/metadata/<ts>/analysis-prep/ground_truth/, .../preprocessed/detection/totals.csv Figs 2-6
DYN cohort cache (derived, de-identified profiles) results/metadata/<ts>/dyn-prep/tables/cohorts/ - see Datasets Figs 7-8

<ts> is the YYYYMMDD_HHMMSS of the analysis_prep.py / dyn_prep.py run that wrote those tables. Multiple runs accumulate side-by-side under results/metadata/; figure notebooks call find_latest_analysis_prep_ts() / find_latest_dyn_prep_ts() from scripts/analysis/utils/results.py to auto-pick the most recent. Pin a specific timestamp by setting TS = "<YYYYMMDD_HHMMSS>" near the top of the notebook config cell.

These tables are what analysis_prep.py (mock + simulated datasets) and dyn_prep.py (DYN cohort) produce; shipping them lets you skip stages 1–3 of the pipeline.

Not metadata-reproducible (require the full pipeline / raw data, or are hand-drawn): Fig 1 (read-length, needs the reads), Fig 9 (resource benchmarking, needs per-run timing/memory logs), Fig 10 (recommended-workflow diagram, hand-drawn), Fig S1 (error rates, needs alignments). Use full reproduction for these.

Quick reproduction

Reproduce Figs 2-8 + Tables S1, S2 from the shipped tables.

Step 1 - Environment

conda activate bakeoff      # see Installation

The quick path only exercises the analysis stack (Python, Jupyter, pandas, scipy, seaborn). The profiler binaries are not needed here - if a profiler failed to install, quick reproduction still works.

Step 2 - Confirm the tables are present

# the most-recent run dir is what notebooks auto-pick
ls -dt results/metadata/2*/ | head -1

# inside that run dir, expect:
TS=$(ls -dt results/metadata/2*/ | head -1)
ls $TS/analysis-prep/preprocessed/detection/    # per-tool detection CSVs
ls $TS/analysis-prep/preprocessed/abundance/    # per-tool abundance CSVs
ls $TS/analysis-prep/ground_truth/              # truth tables
ls $TS/dyn-prep/tables/cohorts/                 # DYN cohort cache (Figs 7-8)

Step 3 - Run the figure notebook

Open the notebook for the panel you want, set the config cell at the top (DATASET, and threshold where applicable), and run all cells. Headless example:

jupyter nbconvert --to notebook --execute scripts/analysis/analysis_detection.ipynb
Fig Content Notebook Config cell
2 Simulated PacBio + ONT detection (P/R/F1) analysis_detection.ipynb DATASET = a sim key
3 Mock detection (P/R/F1) - PacBio standard + ONT Zymo as representatives; supps S2 cover the other two mocks analysis_detection.ipynb DATASET = a mock key
4 Detection operating characteristic across the threshold sweep (all four D6331 mocks pooled) auxiliary/sensitivity_floor.ipynb - (pools all four mocks)
5 Mock relative-abundance stacked bars, PacBio standard + ONT Zymo analysis_abundance.ipynb DATASET = a mock key
6 Mock abundance concordance (ρ, 1-BC, 1-JSD), PacBio standard + ONT Zymo analysis_abundance.ipynb DATASET = a mock key
7 DYN per-sample heatmaps (species + genus) dyn_heatmap.ipynb DATASET, MIN_ABUNDANCE
8 DYN alpha-diversity (Shannon, Simpson) dyn_alpha_div.ipynb cohorts, MIN_ABUNDANCE

Outputs land under results/<ts>/{detection,abundance,dyn-alpha-div,dyn-heatmap}/..., where <ts> is the same timestamp the notebook reused from analysis_prep.py / dyn_prep.py. Each save site is gated by the SAVE = True/False knob at the top of every notebook - set SAVE = False to render inline only.

Table S1

jupyter nbconvert --to notebook --execute scripts/analysis/database_comparison.ipynb

Full reproduction

Regenerate everything from raw reads. Requires the reference databases and read datasets (Data Fetch) and all six profiler binaries (Installation). See Pipeline Overview for the mechanics of each stage.

Step 1 - Activate the environment

conda activate bakeoff

Step 2 - Run the profilers

process.sh accepts any mix of datasets in a single --dataset-list; the output <data_group>/<technology> subpath is derived from each fastq's own path, so mock + simulated + DYN entries can live in one manifest. The canonical "everything" manifest is data/datasets_all.txt - every long-read dataset, with section headers (mock / simulated / DYN) so subsets are an easy comment-out:

# Run every long-read dataset against both DB families
bash scripts/process/process.sh --db unified --tools all --dataset-list data/datasets_all.txt
bash scripts/process/process.sh --db default --tools all --dataset-list data/datasets_all.txt

To run only a subset, either edit datasets_all.txt (comment out the sections you don't need) or use one of the focused per-group manifests already in data/:

Manifest Contents
data/datasets_mock_pacbio.txt PacBio HiFi D6331 mock (Figs 3, 5-6 + supp panels in Figs S2-S5)
data/datasets_mock_ont.txt ONT D6331 mock (Figs 3, 5-6 + supp panels in Figs S2-S5)
data/datasets_simulated.txt Simulated PacBio + ONT (Fig 2)
data/datasets_dyn_pacbio.txt, ..._dyn_ont.txt DYN clinical cohort (Figs 7-8 + supp panels in Figs S6-S8)

Common optional flags:

  • --db {unified|default} - which reference family to profile against.
  • --tools <list> - comma-separated subset (e.g. kraken2,sylph), or all.
  • --dataset <fq> / --dataset-list <file> - one or many input FASTQs. Repeatable; combinable.
  • --db-dir <path> - parent of the per-DB indexes (default data/ref_db); override when the DBs live elsewhere.
  • --threads N - per-tool threads (default 50).

For paired-end Illumina baselines (Fig 8), use scripts/process/process_sr.sh with the same flag conventions and data/datasets_dyn_illumina.txt.

Step 3 - Post-process raw outputs

bash scripts/process/postprocess.sh --db unified --jobs 4
bash scripts/process/postprocess.sh --db default --jobs 4

Converts profiler-native outputs into the kreport / .tre / .sylphmpa files the notebooks read. Common optional flags:

  • --data-groups <csv> - default ZymoMockD6331,simulated (DYN is opt-in; add it to include the clinical reports).
  • --db-dir <path> - same semantics as Step 2.
  • --jobs N - per-sample parallelism inside each sub-script.

Step 4 - Ensemble into the analysis-prep tables

python scripts/analysis/analysis_prep.py --jobs 4

Writes per-(sample, tool, db) CSVs + ground-truth tables to results/metadata/<YYYYMMDD_HHMMSS>/analysis-prep/.... Common optional flags:

  • --mode {detection|abundance|both} - default both.
  • --data-groups <csv> - same default as Step 3.
  • --db-dir <path> - ETE3 sqlite snapshots live under here; override if your DB tree differs.
  • --dry-run - prints the planned (mode, tool, db, sample) rows with [ok] / [MISSING] markers and writes a manifest CSV under results/metadata/dry-run/ without parsing or producing real outputs. Use it before a long run to confirm the report tree matches the registry.

Step 4b - (DYN only) build the cohort cache

python scripts/analysis/dyn_prep.py --threads 8

Required only if you want to reproduce Figs 7-8. Writes to results/metadata/<YYYYMMDD_HHMMSS>/dyn-prep/tables/cohorts/. Supports the same --dry-run semantics as Step 4.

Step 5 - Run the figure notebooks

Same as Quick reproduction Step 3. The notebooks auto-discover the most recent timestamp directory written by Steps 4 / 4b, so no extra wiring is needed.

Figure 1, Figure 9, Figure 10, Figure S1

These sit outside both the quick and standard figure notebooks:

Fig Content Source
1 Read-length distributions read-analysis/read_info.ipynb (reads a FASTQ list)
9 Resource benchmarking (runtime, RAM, F1 vs CPU) resource_benchmark.ipynb (per-tool logs from process_*.sh + read-counts from data/read_stats/read_stats.csv)
10 Recommended-workflow diagram hand-drawn - no script
S1 Per-read error-rate distributions read-analysis/error_analysis.ipynb (upstream: error_rates.py)

Tables

Table Content Source
1 Dataset overview read-analysis/read_stats.py → manual tabulation
2 Methods + default DBs hand-tabulated; no script
S1 Taxonomic breadth across DBs and ranks analysis/database_comparison.ipynb

Customizing (optional)

The figure notebooks are parameterized by their config cell, so you can explore beyond the published settings without touching any other code:

  • DATASET - switch among the mock / simulated keys in the DATASETS dict to render the same analysis for a different sample.
  • THRESHOLD_PERCENT (detection) - change the support threshold (published value: 0.001 %); re-run to get precision/recall/F1 at a different cutoff. Result tables are written under results/<ts>/detection/threshold_<x>/.
  • MIN_ABUNDANCE (abundance, DYN heatmap, alpha-div) - change the abundance floor (published value: 0 for mock abundance, 0.1 % for DYN). Outputs are namespaced by threshold so runs don't overwrite each other.
  • SAVE - at the top of every figure notebook; flip to False to render plots inline without touching disk (useful when iterating).

Each re-run writes a new table/figure namespaced by the threshold value, so exploratory runs never clobber the published outputs. This is purely optional - the default config values reproduce the manuscript.

Sweeping a parameter

scripts/analysis/sweep_notebook.sh sweeps one config variable across a list of values without editing the notebook itself, and optionally repeats the whole sweep across multiple datasets (cross-product, datasets × values) in a single command:

# Detection - 6 datasets × 5 thresholds = 30 runs
scripts/analysis/sweep_notebook.sh \
    --notebook scripts/analysis/analysis_detection.ipynb \
    --var THRESHOLD_PERCENT \
    --datasets "pacbio_low_input pacbio_standard_input ont_zymo_kit ont_qiagen_kit sim_pacbio sim_ont"

# Abundance - 4 datasets × 5 thresholds = 20 runs (no simulated; abundance is mock-only)
scripts/analysis/sweep_notebook.sh \
    --notebook scripts/analysis/analysis_abundance.ipynb \
    --var MIN_ABUNDANCE \
    --datasets "pacbio_low_input pacbio_standard_input ont_zymo_kit ont_qiagen_kit"

Both calls inherit the default --values "0.0 0.0001 0.001 0.01 0.1" (the published sensitivity panel). Override with --values "..." to use a different range.

Use --set 'TS="..."' to pin a specific analysis-prep timestamp instead of letting the notebook auto-pick the latest one; --set is repeatable for additional fixed overrides. The script patches the notebook in-memory per iteration, runs it with jupyter nbconvert, and discards the executed .ipynb copies - the per-iteration result tables and figures land in the notebook's normal results/<TS>/.../<dataset>/ subdirs. --help lists all options.

Execution notes

  • NCBI sqlite: notebooks that resolve taxonomy build NCBITaxa(dbfile=...) against the prebuilt taxa032025.sqlite. Do not pass taxdump_file= - it triggers a minutes-long rebuild that overwrites the pinned snapshot. (The quick-path figure notebooks read pre-resolved tables and don't need this.)
  • Working directory: run from bakeoff/. jupyter nbconvert sets the kernel CWD to the notebook's own folder, so execute in place with the repo as CWD (see FAQ and Troubleshooting).
  • Hardcoded paths: some notebooks carry project-specific path constants in the first cell; update them on a fresh checkout.
  • DB drift: if you do a full rebuild of the unified reference, rebuild every tool index from the same taxdump snapshot - see Database Build.

Getting Started

Reproducing the Paper

Reference

Contributing

Clone this wiki locally