diff --git a/.github/workflows/build_container.yaml b/.github/workflows/build_container.yaml index 23100bd..d179f50 100644 --- a/.github/workflows/build_container.yaml +++ b/.github/workflows/build_container.yaml @@ -24,6 +24,16 @@ jobs: # pre-pulled docker images), freeing ~25-30GB. ``tool-cache: false`` # keeps the runner's node/python toolcache so other actions (checkout, # docker login) still work fast. + # + # NOTE this image is MONOLITHIC -- on top of that base it installs + # miniforge and pre-builds every per-rule conda env. It deliberately does + # NOT build a prostt5-* env: predict_3di reuses the base image's own ROCm + # torch via --gpu-backend system, which keeps several GB out of the image. + # Even so this is close to what a hosted runner can manage after the + # cleanup below. If a build starts failing on "no space left on device", + # the honest fix is building on a machine with real disk and pushing + # manually, not shaving more off here. The ``df -h`` steps exist to make + # that diagnosis obvious rather than a guess. - name: Free up disk space uses: jlumbroso/free-disk-space@main with: @@ -46,6 +56,11 @@ jobs: IMAGE=quay.io/${{ secrets.QUAY_USERNAME }}/phables:${{ github.sha }} docker build -f container/Dockerfile -t $IMAGE . + - name: Show image size and free disk after build + run: | + docker image ls + df -h + - name: Push Docker image run: | IMAGE=quay.io/${{ secrets.QUAY_USERNAME }}/phables:${{ github.sha }} diff --git a/PR_BODY.md b/PR_BODY.md new file mode 100644 index 0000000..1b4ee78 --- /dev/null +++ b/PR_BODY.md @@ -0,0 +1,157 @@ +Follow-up to #68. Three independent strands: a crash fix, two flow-decomposition +speedups, and a rework of the container added in #68. + +Rebases cleanly onto current `develop` (checked against `d196f5f`, including the +`QUAY_NAMESPACE` change to `build_container.yaml`). + +--- + +## 1. Fix: phables crashes when a sample resolves no genomes + +A sample where nothing is resolved and there are no unresolved phage-like edges +produces an empty `genomes_and_unresolved_edges.fasta`. CoverM 0.7.0 doesn't +print an empty table for the resulting zero-alignment BAM — it panics: + +``` +[WARN coverm::contig] No primary alignments were observed for sample X +thread 'main' panicked at src/coverage_printer.rs:467:61: +index out of bounds: the len is 0 but the index is 0 +``` + +That killed the whole run at the very last stage, after all the expensive work +had already succeeded. Hit for real on `SRR19670770`. + +`coverm_bam2counts_genomes` now checks the BAM for alignments first and writes a +header-only coverage table instead of invoking CoverM when there are none. +Everything downstream already handled an empty table correctly, so the run +finishes normally with empty report tables. The header reproduces CoverM's own +exactly, including `Covered Fraction` being two words. + +## 2. Performance: flow decomposition + +Both changes are opt-out-safe — the second is off by default — and both were +profiled before being written rather than guessed at. + +**Where the time actually goes.** For a component's MILP, *building* the model +is ~95% of the cost, not solving it (large component: 170 ms build vs 8.9 ms +solve). Two consequences, both measured: solver threads make no difference at +all (1/2/4/8 threads are flat within noise), and a solver `time_limit` does not +help either. + +**a. Start the K search at a proven lower bound** (`FD_Algorithm`) + +`FD_Algorithm` tried K = 1, 2, 3, … until feasible, rebuilding the whole MILP +each time. K is structural to the model, so it genuinely cannot be reused, and +flowpaths does not expose a HiGHS warm start — meaning every attempt below the +true answer was a full model build that could only return infeasible. + +`get_lowerbound_k()` takes the max of the graph width and +`ceil(log2(#distinct flow values))`, both lifted from flowpaths' own +`MinFlowDecomp.get_lowerbound_k`. Both are lower bounds, so starting there +cannot skip a feasible smaller K. It costs 1–4 ms and falls back to 1 on any +error, since a lower bound is an optimisation and must never be why a component +fails to resolve. + +| | | +|---|---| +| 18/18 synthetic cases | identical K, path count and path sets | +| speedup | 1.5×–4.9×, growing with component size | +| components that can't resolve within `--maxpaths` | up to 5.9× (the bound proves `K >= maxpaths` up front instead of burning the whole ladder) | + +(Also annotates `data["minK"]`, which was set to a constant `2` and never read +by anything, so it isn't mistaken for the live bound.) + +**b. `--mfd-workers`: run components in parallel** (default `1`, unchanged behaviour) + +Components are independent, so the loop is embarrassingly parallel. +`resolve_short_parallel` chunks them, runs the **existing** `resolve_short` once +per chunk in a worker process, and merges the returned accumulators — no change +to that function's ~1400-line body, since it's already parameterised by the +component set and already returns everything it builds. + +This is only sound because no component's logic depends on another's results. +That was verified against the body first: every touch of a shared accumulator is +a pure `add`/`union`/`append`, with no conditional or membership test against +them anywhere in the loop (per-component decisions use the loop-local `comp_*` +sets). It's noted in the docstring, because chunking would silently change +results if that ever stopped being true. + +Chunks merge in component order, so `all_resolved_paths` is identical to the +sequential run — genomes are numbered by position, so a different order would +rename every genome without changing the biology. + +| Workload | 2 | 4 | 8 workers | +|---|---|---|---| +| uniform components | 1.93× | 3.74× | 6.08× | +| realistic skew | 1.92× | 2.02× | 2.35× | + +The skewed case is the one to plan around: one component was 41% of total +runtime, giving an Amdahl ceiling of 2.4× — so 2.35× is ~98% of what's +achievable. Two details mattered: more chunks than workers (so the pool can +balance an uneven workload; this alone took 8 workers from 1.92× to 2.35×), and +sending the heavy read-only inputs once per worker via a pool initializer rather +than once per chunk, so smaller chunks don't mean re-pickling the assembly graph. + +Verified identical results at 2/4/8 workers, identical path ordering, and +identical output for 1 component, fewer components than workers, more components +than workers, and `workers=1`. + +## 3. Container: one monolithic image, replacing the per-rule container flags + +**This removes `--container` and `--prostt5-container`, both added in #68.** +Worth being explicit about, since they were merged only recently. + +Those flags pointed individual rules at images via Snakemake `container:` +directives. That approach fights Snakemake: combining `container:` with `conda:` +triggers the documented "ad-hoc combination" behaviour, which builds a *fresh* +conda env *inside* the container rather than using what the image already has — +defeating the point. The rules now carry plain `conda:` directives only, and the +image satisfies them by having every per-rule env **pre-built inside it** +(hybracter's approach). Run it and there is nothing left to create at runtime, +which removes the conda-env creation race that motivated containerising at all: +each sample is its own Snakemake process, so concurrent array tasks needing the +same not-yet-built env can corrupt each other's `mamba env create`. + +Also here: +- **`--gpu-backend system`** — builds no env and uses the ambient `torch` + + `pholdlib`. Conda envs are isolated, so a rule declaring `conda:` can never see + a torch installed outside it; this is the only way to *reuse* a known-good GPU + torch (the container's ROCm base, or a module-loaded torch on HPC) instead of + installing a second copy. Keeps several GB out of the image. +- `container/prebuild_envs.sh` generates its own throwaway inputs rather than + using `tests/data/` (which is gitignored, so absent in a fresh clone or CI + checkout), and points `phables install` at an *empty* databases dir so the + download rules are actually in the DAG and the `curl` env gets built. +- `container/test_image.sh` fails the build if any per-rule tool is missing, if + `torch`/`pholdlib` aren't importable, or if any pre-built env contains its own + torch (i.e. a second copy crept in). +- `prostt5-rocm.yaml` bumped to `torch==2.9.1` (verified present on the pinned + `rocm6.3` index, cp310–cp314). + +--- + +## Testing + +- Existing test suite passes. +- Workflow DAG verified by real `--dry-run`s across gene caller, detection mode, + GPU backend and tree options. +- Performance and equivalence numbers above are from synthetic instances driving + the real `FD_Inexact` / `resolve_short` code paths. +- The container built end-to-end: all 8 per-rule conda envs solve and install, no + duplicate torch, and the pip install leaves the base image's torch untouched. + +## Not verified + +- No end-to-end run on a real sample through the parallel path — the equivalence + testing used a stubbed `resolve_short` (real chunking and merging, synthetic + per-component results). +- Speedups are from synthetic components with clean topology. Build-dominates- + solve should hold generally since it tracks graph size, but a genuinely hard + component would shift the ratio. The lower bound is *valid* regardless, just + possibly looser on real graphs. +- No Setonix Apptainer run of the final image. +- Process pools copy memory per worker, so a large assembly graph at 8 workers + may bind on RAM before CPU. + +Happy to split any of the three strands into its own PR if that's easier to +review — they're independent. diff --git a/container/Dockerfile b/container/Dockerfile index c5bcb5c..0904217 100644 --- a/container/Dockerfile +++ b/container/Dockerfile @@ -1,137 +1,103 @@ # -# phables -- massive all-in-one image for Setonix (ROCm), replacing per-rule -# --use-conda env creation entirely. See ../docs/container.md for how to run -# it (--container --use-singularity --no-use-conda) and why this -# exists (a real, hit-in-production race condition: every sample's `phables -# run` is its own separate process, so many array tasks needing the same -# shared --use-conda env for the first time at once can corrupt each other's -# `mamba env create`). +# phables -- ONE monolithic image for Setonix: ROCm base + conda + every +# per-rule conda env pre-built INSIDE the image. Build it once, convert it to a +# single .sif, and every `phables run` inside it finds its envs already there. +# Same strategy as hybracter's own container (miniforge installed into the +# image, envs materialised at build time), NOT the per-rule `container:` +# directive approach -- phables has no --container/--prostt5-container flags +# any more, and its Snakemake rules carry plain `conda:` directives only. +# +# Why this exists: a real, hit-in-production race condition. Every sample's +# `phables run` is its own separate process (one SLURM array task each), and +# Snakemake only serialises conda-env creation WITHIN a single process's DAG. +# Many array tasks all needing the same not-yet-built shared env at once +# corrupt each other's `mamba env create` ("Fatal Python error: +# init_fs_encoding ... no codec search functions registered"). With every env +# baked into the image there is nothing left to create at runtime, and so +# nothing left to race. +# +# TORCH: this image does NOT install a second PyTorch. The Pawsey ROCm base +# already ships a torch that is known to work on Setonix's MI250X (gfx90a), so +# predict_3di reuses it via `--gpu-backend system` (see below) instead of a +# prostt5-* conda env. That is the ONLY way to reuse it: conda envs are +# isolated, so a rule declaring `conda:` can never see a torch installed +# outside its env -- which is exactly why an earlier version of this image +# ended up downloading multiple GB of a second, unverified torch. +# +# The consequence for how the image is built: phables and Snakemake are +# installed into the BASE IMAGE'S python, not into miniforge's. Snakemake runs +# `script:` rules that declare no conda env using its own interpreter, so that +# interpreter has to be the one holding torch. Miniforge is installed only to +# provide the `conda` binary that BUILDS the other per-rule envs, and is +# appended to PATH rather than prepended so it can never shadow the base +# python. +# +# See ../docs/container.md for how to build the .sif and run it on Setonix. # -# Base is Pawsey's own verified-working ROCm+PyTorch image -- the SAME one -# phold's own container (../../phold/container/hpci/Dockerfile) builds from, -# deliberately: torch==2.7.1+rocm6.3 baked in here is exactly what -# envs/prostt5-rocm.yaml pins, confirmed working on real Setonix MI250X -# hardware this session. Nothing below reinstalls or touches torch -- pip -# installs pholdlib and every other pure-Python dependency directly into this -# same system Python, so they see (and don't conflict with) the base image's -# already-correct torch/ROCm stack rather than resolving a second, possibly -# different one. FROM quay.io/pawsey/pytorch:2.7.1-rocm6.3.3 ARG DEBIAN_FRONTEND="noninteractive" RUN apt-get update && apt-get install -y --no-install-recommends \ - wget \ - tar \ - build-essential \ + wget \ + bzip2 \ + ca-certificates \ + git \ + procps \ && rm -rf /var/lib/apt/lists/* -# the base image only provides python3, not python -- some tooling (and -# Snakemake's own script: directive) expects `python` on PATH -RUN ln -sf "$(command -v python3)" /usr/local/bin/python +# Pin down the base image's python BEFORE miniforge exists, so later steps can +# target it unambiguously no matter where the base keeps it. Everything phables +# needs goes into this interpreter. +RUN command -v python3 > /opt/base_python_path && \ + ln -sf "$(cat /opt/base_python_path)" /usr/local/bin/python && \ + echo "base python: $(cat /opt/base_python_path)" && \ + "$(cat /opt/base_python_path)" -c "import torch; print('base torch:', torch.__version__)" # -# foldseek -- same install as phold's own Dockerfile (direct binary download, -# not conda/pip), since that's the reference build this image otherwise -# mirrors for the ROCm/torch stack +# Miniforge -- installed ONLY for its `conda`/`mamba` binaries, which Snakemake +# shells out to when building the per-rule envs. Deliberately APPENDED to PATH: +# prepending it would make miniforge's python the default, and phables would +# then run under an interpreter with no torch in it, silently defeating the +# whole point of building on a ROCm base. # -RUN wget https://mmseqs.com/foldseek/foldseek-linux-avx2.tar.gz && \ - tar -xzf foldseek-linux-avx2.tar.gz && \ - rm foldseek-linux-avx2.tar.gz && \ - mv foldseek /opt/foldseek -ENV PATH="/opt/foldseek/bin:${PATH}" +ENV CONDA_DIR=/opt/miniforge3 +RUN wget -q https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh -O /tmp/miniforge.sh && \ + bash /tmp/miniforge.sh -b -p "${CONDA_DIR}" && \ + rm /tmp/miniforge.sh +ENV PATH="${PATH}:${CONDA_DIR}/bin" # -# Every other per-rule tool that isn't a pure-Python package: bootstrap a -# standalone micromamba (no full miniforge needed) and install them all into -# one prefix, /opt/conda -- kept entirely separate from the base image's own -# system Python (added to PATH, but conda never pulls in a python of its own -# here since none of these packages are Python libraries, only standalone -# binaries -- confirmed by their real recipes: minimap2/samtools/mmseqs2/ -# fraggenescan/hmmer/mafft have no python dependency; coverm's own recipe -# pulls minimap2+samtools transitively too, listed here explicitly anyway for -# a single deterministic solve rather than relying on that side effect). -# mmseqs2 pinned to the exact version envs/mmseqs.yaml already uses. -RUN curl -Ls https://micro.mamba.pm/api/micromamba/linux-64/latest | tar -xvj bin/micromamba && \ - mv bin/micromamba /usr/local/bin/micromamba && \ - micromamba create -y -p /opt/conda -c conda-forge -c bioconda \ - minimap2 \ - samtools \ - "coverm>=0.6" \ - "mmseqs2=13.45111" \ - fraggenescan \ - hmmer \ - mafft \ - && micromamba clean -a -y -ENV PATH="/opt/conda/bin:${PATH}" - -# -# Everything else -- phables' own pure-Python deps, ProstT5/pholdlib, and -# phables itself -- pip installed directly into the base image's system -# Python. Deliberately NOT a separate venv/conda env: the whole point is for -# all of this to see the base image's already-installed, already-verified -# torch/ROCm build rather than resolving (and risking conflicting with) a -# second one. +# phables + Snakemake + pholdlib into the BASE python. pholdlib declares no +# torch version constraint, so pip pulls its own dependencies (transformers, +# sentencepiece, h5py, ...) without touching the already-installed torch -- the +# assertion below fails the build if that ever stops being true, since a +# silently-replaced torch is precisely the failure this image exists to avoid. # -# pholdlib's own declared torch requirement is unconstrained (confirmed -# earlier this session against envs/prostt5-rocm.yaml/-cpu.yaml/-cuda.yaml) -- -# `pip install pholdlib` alone pulls in its transitive deps (transformers, -# sentencepiece, h5py, loguru, etc.) via pip's own resolver without touching -# the already-installed torch. -RUN pip install --ignore-installed PyYAML && \ - python -c "import torch; print('torch (pre-existing, untouched):', torch.__version__)" +COPY . /opt/phables_src +RUN BASE_PY="$(cat /opt/base_python_path)" && \ + TORCH_BEFORE="$("$BASE_PY" -c 'import torch; print(torch.__version__)')" && \ + "$BASE_PY" -m pip install --no-cache-dir "snakemake>=8" pholdlib /opt/phables_src && \ + TORCH_AFTER="$("$BASE_PY" -c 'import torch; print(torch.__version__)')" && \ + echo "torch before=${TORCH_BEFORE} after=${TORCH_AFTER}" && \ + test "$TORCH_BEFORE" = "$TORCH_AFTER" # -# NOT included: cogent3/piqtree (phylotree.yaml's own deps, for the optional -# --phylotree/build_tree rule -- not part of this workflow's usage). -# Deliberately left out rather than fought into pip (every cogent3 release on -# PyPI is an unreleased alpha, which breaks a plain version-range pip -# install). Note: --phylotree is therefore NOT supported under --container -- -# build_tree still declares container: CONTAINER_IMAGE (unchanged), so it -# will fail if actually invoked against this image. +# NOTE: there is deliberately NO step here patching config.yaml's gpu_backend +# to `system`. An earlier version did exactly that, and it does not work: +# phables passes every CLI option to snaketool as merge_config=kwargs, so any +# option with a non-None click default is merged OVER the config file, and +# --gpu-backend's click default is "cpu". The patched config was silently +# ignored, Snakemake went to build a prostt5-cpu env that isn't in this image, +# and the run died on the read-only .sif filesystem. # -RUN pip install \ - "pholdlib>=0.1.2" \ - "pyrodigal-gv>=0.3.1" \ - biopython \ - python-igraph \ - pysam \ - "networkx>=2.8.6" \ - scipy \ - numpy \ - pandas \ - tqdm \ - click \ - "metasnek>=0.0.3" \ - "flowpaths>=0.2.20" - -RUN python -c "import torch; print('torch (post pip installs, should be unchanged):', torch.__version__)" - +# Callers must therefore pass `--gpu-backend system` explicitly. See +# ../docs/container.md; mass_processing/run_phables_container.sh does this. # -# phables itself -- installed from the exact source in this build context -# (this Dockerfile is meant to be built by CI on every commit, tagged by -# commit SHA, so it's the checked-out commit's own tree, not a fetched -# git+https ref pinned separately). -# -COPY . /opt/phables_src -RUN pip install /opt/phables_src -# -# tests -- fail the build loudly if any of the above didn't actually work, -# rather than shipping a broken image to quay.io -# -RUN foldseek --help -RUN minimap2 --version -RUN samtools --version -RUN mmseqs version -RUN coverm --version -RUN hmmsearch -h -# the bioconda package's actual binary is run_FragGeneScan.pl, not -# `fraggenescan` -- confirmed against genes.smk's own real invocation -RUN run_FragGeneScan.pl 2>&1 | head -5 || true -RUN mafft --version -RUN python --version -RUN python -c "import pholdlib; print('pholdlib OK')" -RUN python -c "import torch; assert torch.__version__.startswith('2.7.1'), torch.__version__" -RUN phables --version -RUN phables run -h +# Pre-build every per-rule conda env. See the script for the full rationale +# (default --conda-prefix, placeholder databases, one pass per flag combo). +RUN bash /opt/phables_src/container/prebuild_envs.sh + +# Fail the build loudly rather than pushing a broken image to quay.io. +RUN bash /opt/phables_src/container/test_image.sh diff --git a/container/prebuild_envs.sh b/container/prebuild_envs.sh new file mode 100755 index 0000000..abbd310 --- /dev/null +++ b/container/prebuild_envs.sh @@ -0,0 +1,99 @@ +#!/bin/bash +# Pre-builds EVERY per-rule conda env into the image, at Docker build time. +# Run from container/Dockerfile; not meant to be run on a host. Takes no +# arguments and reads nothing from the repo -- it generates its own throwaway +# inputs (see below), so it works in a bare CI checkout. +# +# Envs go into phables' own DEFAULT --conda-prefix (snake_base("workflow/conda"), +# i.e. inside the installed package) deliberately: a user inside the resulting +# container then runs plain `phables run ...` with no --conda-prefix flag and +# Snakemake resolves these exact envs. Snakemake names an env by hashing its +# file content together with the conda prefix path -- both identical at build +# time and run time here, since they're the same paths in the same image -- so +# the hashes match and nothing is rebuilt at runtime. That matters more than +# it sounds: a .sif is READ-ONLY when running, so an attempted rebuild is a +# hard failure, not just a slow path. +# +# --conda-create-envs-only builds a DAG's envs without running any of it. The +# DAG still has to RESOLVE, which requires the database files to EXIST -- but +# only to exist, since no job runs. Empty placeholder files are therefore +# enough. That was verified for real (dry-running every flag combination below +# against zero-byte placeholder DB files) before this script was written, and +# it's what keeps the multi-GB PHROGs/hallmark databases OUT of the image: +# mount the real ones at runtime via --databases, exactly as outside a +# container. +# +# One invocation per flag combination, because which envs a DAG needs depends +# on the flags -- gene caller, phage-detection mode and GPU backend each select +# different rules and env files. Together these cover every env under +# workflow/envs/ that any `phables run` (or `phables install`) can reach. + +set -euxo pipefail + +DB=/tmp/placeholder_db +WORK=/tmp/envbuild_inputs + +mkdir -p "$DB/phrogs_mmseqs_db" "$DB/hallmark_db" +touch "$DB/marker.hmm" \ + "$DB/phrog_annot_v4.tsv" \ + "$DB/phrogs_mmseqs_db/phrogs_profile_db" \ + "$DB/hallmark_db/hallmark_db" \ + "$DB/hallmark_db/hallmark_categories.tsv" + +# Synthetic minimal inputs, generated here rather than taken from +# tests/data/: that directory is in .gitignore, so it does NOT exist in a +# fresh clone or in a CI checkout -- depending on it made the image build +# fail with "Invalid value for '--reads': Path ... does not exist". Nothing +# below is ever actually processed (no job runs under +# --conda-create-envs-only); these files exist purely so click's exists=True +# checks pass and the DAG can resolve. Verified to produce the same DAG as +# the real test data. +mkdir -p "$WORK/reads" +printf 'H\tVN:Z:1.0\nS\tedge_1\tACGTACGTACGTACGTACGTACGTACGTACGT\tLN:i:32\nS\tedge_2\tTTTTGGGGCCCCAAAATTTTGGGGCCCCAAAA\tLN:i:32\nL\tedge_1\t+\tedge_2\t+\t0M\n' \ + > "$WORK/assembly_graph.gfa" +printf '@r1\nACGT\n+\nIIII\n' | gzip > "$WORK/reads/sample1_R1.fastq.gz" +printf '@r1\nACGT\n+\nIIII\n' | gzip > "$WORK/reads/sample1_R2.fastq.gz" + +GFA="$WORK/assembly_graph.gfa" +READS="$WORK/reads" +COMMON=(--input "$GFA" --reads "$READS" --databases "$DB" --threads 1) + +# 1. default path -> coverm, genecall (FragGeneScan), smg (HMMER), mmseqs, phables +phables run "${COMMON[@]}" --output /tmp/envbuild1 --conda-create-envs-only + +# 2. the other gene caller +phables run "${COMMON[@]}" --output /tmp/envbuild2 \ + --genecaller pyrodigal-gv --conda-create-envs-only + +# 3. ProstT5 + foldseek detection -> the foldseek env (and, via +# --gpu-backend system, NO prostt5-* torch env at all: predict_3di reuses the +# base image's already-working ROCm torch, which is installed in the same +# python running Snakemake here). This is the whole reason no multi-GB +# second torch is downloaded during this build. Deliberately NOT `rocm`/ +# `cpu`/`cuda` -- each of those would solve and download their own torch. +phables run "${COMMON[@]}" --output /tmp/envbuild3 \ + --phagedetection prostt5-foldseek --gpu-backend system --conda-create-envs-only + +# 4. Optional phylogenetic tree -> phylotree env (MAFFT + cogent3/piqtree). +# Note conda resolves cogent3 fine, unlike pip, where every published +# release is a prerelease and a plain version range matches nothing. +phables run "${COMMON[@]}" --output /tmp/envbuild4 \ + --build-tree --conda-create-envs-only + +# 5. install.smk's own env (curl), so `phables install` works inside here too. +# Deliberately pointed at an EMPTY databases dir, not "$DB": the placeholder +# files in $DB satisfy install.smk's own download targets, so Snakemake says +# "Nothing to be done", the DAG is empty, and NO env gets created -- the +# curl env would silently be missing from the image. An empty dir puts the +# four *_download rules in the DAG so their conda env actually gets built. +# (--conda-create-envs-only still downloads nothing.) +mkdir -p /tmp/empty_db +phables install --output /tmp/envbuild5 --databases /tmp/empty_db --conda-create-envs-only + +rm -rf /tmp/envbuild1 /tmp/envbuild2 /tmp/envbuild3 /tmp/envbuild4 /tmp/envbuild5 \ + "$DB" "$WORK" /tmp/empty_db +conda clean -a -y + +echo "=== pre-built conda envs ===" +CONDA_PREFIX_DIR="$(python -c 'import phables, os; print(os.path.join(os.path.dirname(phables.__file__), "workflow", "conda"))')" +ls -1 "$CONDA_PREFIX_DIR" diff --git a/container/test_image.sh b/container/test_image.sh new file mode 100755 index 0000000..eb5e04b --- /dev/null +++ b/container/test_image.sh @@ -0,0 +1,125 @@ +#!/bin/bash +# Build-time smoke tests for the monolithic image -- fail the build loudly +# rather than pushing something broken to quay.io. Run from container/Dockerfile. +# +# Every per-rule tool is checked INSIDE whichever pre-built conda env provides +# it, not on the base PATH: none of them are on the base PATH, and that's the +# point of the per-rule env layout. Envs are located by searching the prefix +# rather than by hardcoded directory names, because Snakemake names each env +# by a content hash that this script has no business predicting. + +set -euo pipefail + +PREFIX="$(python -c 'import phables, os; print(os.path.join(os.path.dirname(phables.__file__), "workflow", "conda"))')" +echo "conda prefix: $PREFIX" +test -d "$PREFIX" + +echo "=== phables CLI ===" +phables --version +phables run -h > /dev/null +phables install -h > /dev/null +echo "CLI OK" + +echo "=== the --container / --prostt5-container flags must be GONE ===" +if phables run -h 2>&1 | grep -qE '\-\-(prostt5-)?container'; then + echo "ERROR: a container flag is still present in the CLI" >&2 + exit 1 +fi +echo "confirmed absent" + +echo "=== pre-built env count ===" +n=$(find "$PREFIX" -maxdepth 1 -mindepth 1 -type d | wc -l) +echo "found $n env directories" +find "$PREFIX" -maxdepth 1 -mindepth 1 -type d -exec basename {} \; +# 6 distinct env files are reachable at minimum: coverm, genecall, smg, mmseqs, +# phables, curl -- before counting foldseek and phylotree. No prostt5-* env is +# expected (or wanted): gpu_backend=system reuses the base image's torch. +test "$n" -ge 6 + +# Finds an executable in any pre-built env; fails if no env provides it. +check_bin() { + local want="$1" e + for e in "$PREFIX"/*/; do + if [ -x "${e}bin/${want}" ]; then + echo "OK: $want -> $e" + return 0 + fi + done + echo "MISSING from every pre-built env: $want" >&2 + return 1 +} + +echo "=== per-rule binaries ===" +check_bin minimap2 +check_bin samtools +check_bin coverm +check_bin mmseqs +check_bin foldseek +check_bin hmmsearch +# the bioconda package's real binary name, confirmed against genes.smk's own +# invocation -- NOT `fraggenescan` +check_bin run_FragGeneScan.pl +check_bin mafft +check_bin curl + +echo "=== torch must come from the BASE image, reused -- not reinstalled ===" +# predict_3di runs with NO conda env (gpu_backend=system), i.e. in the same +# python that runs Snakemake -- which must therefore be the base image's python, +# the one already holding a working torch. Two things are checked: that torch + +# pholdlib are importable here (fatal if not -- predict_3di simply cannot run), +# and that no conda env carries a torch of its own (fatal -- that would mean a +# second copy got installed after all). What the torch's build flavour is, is +# only reported. +# +# Reported, NOT asserted. This deliberately cannot fail the build. +# +# The base image's torch is Pawsey's own source build for Setonix -- it reports +# e.g. "2.7.1a0+gite2d141d", with no "+rocm6.3" suffix, because it isn't a +# stock wheel. An earlier version of this script tested for the substring +# "rocm" in torch.__version__ and failed a perfectly good ROCm build at the +# very last step of a multi-GB image. Whether that torch is ROCm-enabled is +# Pawsey's business, not something worth re-litigating here at build time, so +# this prints the facts (torch.version.hip is the real signal: set to the HIP +# version on ROCm builds, None otherwise) and moves on. +# +# The check that DOES matter -- that nothing installed a second torch -- is +# below and is fatal. +python - <<'PY' +import torch + +hip = getattr(torch.version, "hip", None) +print("torch:", torch.__version__) +print("torch.version.hip:", hip) +print("torch.version.cuda:", getattr(torch.version, "cuda", None)) +if hip is None: + print("WARNING: torch.version.hip is None -- this torch does not look " + "ROCm-enabled. Fine if intentional (e.g. a CPU-only base image); " + "worth a look if you expected GPU ProstT5 on Setonix.") +PY +python -c "import pholdlib; print('pholdlib OK')" +python -c "import phables, snakemake; print('phables + snakemake share this interpreter')" + +# A torch inside any pre-built env means a second copy got installed after all, +# which is the exact regression this layout exists to prevent. +for e in "$PREFIX"/*/; do + if [ -x "${e}bin/python" ] && "${e}bin/python" -c "import torch" 2>/dev/null; then + echo "ERROR: a pre-built conda env contains its own torch: $e" >&2 + echo " gpu_backend=system should mean no prostt5-* env is built." >&2 + exit 1 + fi +done +echo "OK: no pre-built env ships a duplicate torch" + +echo "=== --gpu-backend system must be available ===" +# This, NOT a grep of config.yaml's gpu_backend value. An earlier version of +# this script asserted the config file said "system" and passed happily, while +# actual runs still used cpu: phables merges every CLI option over the config +# (merge_config=kwargs), and --gpu-backend's click default is "cpu". The config +# file value is simply not the effective value, so checking it proves nothing. +# What the image can meaningfully guarantee is that the choice EXISTS -- the +# caller is responsible for passing it (see docs/container.md). +phables run -h 2>&1 | grep -q -- '--gpu-backend .*system' \ + || { echo "ERROR: this phables has no --gpu-backend system choice" >&2; exit 1; } +echo "OK: --gpu-backend system is accepted" + +echo "=== all image tests passed ===" diff --git a/docs/container.md b/docs/container.md index 4adc96f..8982129 100644 --- a/docs/container.md +++ b/docs/container.md @@ -1,75 +1,104 @@ # Running phables from a single container -`--use-conda`'s per-rule environment creation is convenient but has a real -failure mode at HPC scale: every rule's conda env is named by a hash of its -env file, so every sample's `phables run` needs the *same* env directory — -but each sample is its own separate Snakemake process (a SLURM array task, -say), and Snakemake only serializes env creation *within* one process's DAG. -Two array tasks both needing the same not-yet-built env at the same moment -can race `mamba env create` into the same directory, corrupting it -(`Fatal Python error: init_fs_encoding ... no codec search functions -registered` is what that looks like when it happens). - -`container/Dockerfile` builds a single image with every per-rule tool this -workflow needs already installed — `--container` then points the whole -workflow at it, replacing conda entirely so there's nothing left to race. +`container/Dockerfile` builds **one monolithic image**: a ROCm base, conda +installed inside it, and every per-rule conda env this workflow can reach +already built into the image. Convert it to a single `.sif` and each +`phables run` inside it finds its environments already there. + +There are no container-related CLI flags. Earlier versions had `--container` +and `--prostt5-container`, which pointed individual Snakemake rules at images +via `container:` directives; both are **removed**. The rules carry plain +`conda:` directives only, and the container satisfies them by having the envs +pre-built rather than by bypassing conda. + +## Why + +`--use-conda`'s per-rule env creation has a real failure mode at HPC scale. +Every rule's conda env is named by a hash of its env file, so every sample's +`phables run` needs the *same* env directory — but each sample is its own +separate Snakemake process (a SLURM array task), and Snakemake only serialises +env creation *within* one process's DAG. Two array tasks both needing the same +not-yet-built env at the same moment can race `mamba env create` into the same +directory and corrupt it (`Fatal Python error: init_fs_encoding ... no codec +search functions registered` is what that looks like in the wild). + +Baking the envs into the image removes the failure mode entirely: at runtime +there is nothing left to create, so nothing left to race. ## Getting an image -CI (`.github/workflows/build_container.yaml`) builds and pushes an image to -`quay.io//phables:` on every push and pull request -— tag by the exact commit you want, there's no floating `latest`. +CI (`.github/workflows/build_container.yaml`) builds and pushes to +`quay.io//phables:` on every push and pull request — +tag by the exact commit you want; there is no floating `latest`. -## Running with it +On Setonix, pull it once as a `.sif`: ```bash -phables run --input assembly_graph.gfa --reads fastq \ - --container quay.io/gbouras13/phables: \ - --no-use-conda \ - --use-singularity +module load singularity/4.1.0-slurm +singularity pull phables.sif docker://quay.io/gbouras13/phables: ``` -`--use-singularity` is Snakemake's own flag for running rules inside -Apptainer/Singularity (the standard container runtime on HPC, incl. Setonix) -— pass it through phables' existing `snake_args` passthrough, same as -`--dry-run`/`--keep-going`/etc. `--use-conda` alone does **not** activate -container execution; explicitly turn it off (`--no-use-conda`) when using -`--container`, since combining the two makes Snakemake build a *separate* -conda env *inside* the container instead of using what's already -installed there — exactly the redundant-torch-reinstall problem this image -exists to avoid. - -## What's actually in the image - -Base: `quay.io/pawsey/pytorch:2.7.1-rocm6.3.3` — the same verified-working -ROCm+PyTorch build `envs/prostt5-rocm.yaml` pins, and the same base phold's -own container (`../../phold/container/hpci/Dockerfile`) builds from. Nothing -in phables' own Dockerfile reinstalls or touches torch — `pip install -pholdlib` and everything else pure-Python go straight into that same system -Python, so they use the base image's already-correct torch/ROCm stack rather -than resolving a second, possibly conflicting one. This is the direct answer -to "don't rebuild the ProstT5 conda env for this container" — there's no -separate ProstT5 env in the container at all; predict_3di runs against the -same Python everything else does. - -Everything else — foldseek (direct binary, matching phold's own install), -minimap2/samtools/coverm/mmseqs2/FragGeneScan/HMMER/MAFFT (via a standalone -micromamba install into `/opt/conda`, kept deliberately separate from the -system Python so it can't touch it), and phables' own Python dependencies -(pyrodigal-gv, biopython, python-igraph, pysam, flowpaths, cogent3/piqtree, -...) — covers every `envs/*.yaml` this workflow has, so the whole DAG can run -from this one image with `--no-use-conda`. - -## `--container` vs `--prostt5-container` - -`--prostt5-container` already existed for pointing *just* `predict_3di` at a -container (e.g. phold's own image, which also has pholdlib+torch). It still -works exactly as before, and if set, it wins for that one rule. `--container` -is new and broader: it's the default container for **every** rule, including -predict_3di if `--prostt5-container` isn't also given. In practice, setting -`--container` alone is enough — there's no reason to use `--prostt5-container` -separately unless you specifically want predict_3di on a *different* image -than the rest of the workflow. +## Running it + +Run phables *inside* the container. Nothing special is passed to phables +itself — `--use-conda` is its own default and the envs are already present: + +```bash +singularity exec --rocm \ + -B /scratch/pawsey1018:/scratch/pawsey1018 \ + phables.sif \ + phables run --input assembly_graph.gfa --reads fastq \ + --output phables_out \ + --databases /scratch/.../all_databases/databases \ + --phagedetection prostt5-foldseek \ + --gpu-backend system \ + --prostt5-checkpoint /scratch/.../model.pt \ + --threads 8 +``` + +Notes that matter on Setonix: + +- **`--rocm`** exposes the host's GPU devices to the container. Without it, + ProstT5 silently falls back to CPU. +- **Bind-mount your scratch** so databases, inputs and outputs are visible. + Databases are deliberately *not* in the image (see below). +- **Don't pass `--conda-prefix`.** The envs were built at the default prefix + (inside the installed package), and Snakemake resolves an env by hashing its + file content *together with the prefix path* — changing the prefix changes + the hash, and Snakemake would try to rebuild into a read-only filesystem. +- **Pass `--gpu-backend system` explicitly.** This is required, not optional. + phables merges every CLI option over the config file + (`merge_config=kwargs`), and `--gpu-backend`'s click default is `cpu` — so + the image's own `config.yaml` value is *not* the effective value and cannot + be relied on. Omitting the flag was tried and failed for real: the runtime + config showed `gpu_backend: cpu`, Snakemake went to build a `prostt5-cpu` + env that isn't in the image, and the run died with + `OSError: [Errno 30] Read-only file system`. Any backend other than + `system` fails the same way. + +## What's in the image + +- **Base**: `quay.io/pawsey/pytorch:2.7.1-rocm6.3.3`, Pawsey's own verified + ROCm build, supplying both the ROCm userspace matching Setonix's MI250X + (gfx90a) **and the torch the workflow actually uses**. +- **phables, Snakemake and pholdlib installed into that base python** — not + into miniforge's. Snakemake runs a `script:` rule that declares no conda env + using its own interpreter, so that interpreter has to be the one holding + torch. The build asserts torch's version is unchanged across the pip install, + since a silently-replaced torch is the exact failure this avoids. +- **Miniforge** at `/opt/miniforge3`, appended to `PATH` (never prepended, so + it cannot shadow the base python). It exists only to provide the `conda` + binary that builds the envs below. +- **Every per-rule conda env**, prebuilt by `container/prebuild_envs.sh`: + coverm (minimap2/samtools/CoverM), genecall (FragGeneScan), pyrodigal-gv, + smg (HMMER), mmseqs, foldseek, phylotree (MAFFT + cogent3/piqtree), and curl + for `phables install`. **No `prostt5-*` env** — that's the point of + `system`. + +**Databases are not included** — PHROGs and the hallmark DB are multi-GB and +separately versioned. Mount them and point `--databases` at them, exactly as +outside a container. `phables install` also works inside the image if you'd +rather fetch them from there. ## Building it yourself @@ -77,16 +106,59 @@ than the rest of the workflow. docker build -f container/Dockerfile -t phables:local . ``` -The base image is large (~14GB compressed) — building on a standard GitHub -Actions runner needs the disk-cleanup step already in -`build_container.yaml` (`jlumbroso/free-disk-space`), or the very first -`FROM` line fails with `no space left on device`. Building locally needs -proportionate free disk space too. - -**Not yet built or run for real** — the Dockerfile and the workflow-wide -`--container` plumbing were written and the Snakemake DAG wiring was verified -via real `--dry-run`s (with and without `--container` set, confirming both -resolve identically at the structural level), but no actual `docker build` -or real Apptainer/Singularity execution against Setonix hardware has -happened yet. Build it, push a first tag, and run one real sample through it -before trusting this for production batches. +`container/prebuild_envs.sh` does the env pre-building. It generates its own +throwaway inputs — a two-segment GFA and a pair of tiny gzipped FASTQs, plus +zero-byte placeholder database files — purely so the Snakemake DAG can +*resolve*, then runs `--conda-create-envs-only` once per flag combination (gene +caller, detection mode, GPU backend, tree) so every reachable env gets built. +Nothing is ever processed, since no job runs. + +Two things it deliberately does **not** do, both of which broke a real build: + +- It doesn't use `tests/data/`. That directory is in `.gitignore`, so it + doesn't exist in a fresh clone or a CI checkout — depending on it failed with + `Invalid value for '--reads': Path ... does not exist`. +- It points `phables install` at an *empty* databases directory, not the + placeholder one. The placeholders satisfy install.smk's own download targets, + so Snakemake reports "Nothing to be done", the DAG is empty and the `curl` + env is silently never built. +`container/test_image.sh` then smoke-tests the result and fails the build if any +env is missing a binary the rules actually invoke, if torch/pholdlib aren't +importable in the ambient python, or if any pre-built env turns out to contain a +torch of its own — that last one being the regression that would mean a second +copy got installed after all. + +It reports the base torch's build flavour (`torch.version.hip`) but does **not** +assert on it. The Pawsey base is a source build reporting e.g. +`2.7.1a0+gite2d141d`, with no `+rocm6.3` suffix; an earlier version of this +script tested for the substring `rocm` and failed a perfectly good ROCm build at +the very last step of a multi-GB image. What that torch is compiled against is +Pawsey's business. + +**Disk**: this image is large — a ~14GB compressed ROCm base plus the conda +envs. The CI workflow runs `jlumbroso/free-disk-space` first because a stock +GitHub Actions runner has only ~14GB free and the base alone won't fit. Reusing +the base torch rather than installing a second one keeps several GB off the +total, but this is still close to the limit of what a hosted runner can build; +if CI starts failing on `no space left on device`, building on a machine with +real disk and pushing manually is the fallback. + +## Status + +A real `docker build` now gets all the way through the image and into +`test_image.sh`. Confirmed working on a real build: + +- All **8** per-rule conda envs solve and install (coverm, genecall, smg, + mmseqs, foldseek, phylotree, phables, curl) — including the pins that were + only repodata-checked before (`mmseqs2=13.45111`, `cogent3<2026.7`). +- **No `prostt5-*` env is created**, so no second torch is downloaded — the + torch reuse works as designed. +- Installing phables/Snakemake/pholdlib into the base python leaves its torch + untouched (the build's own before/after assertion passed). +- Every per-rule binary resolves inside a pre-built env. + +Still unverified: a Setonix Apptainer run — in particular whether ProstT5 +actually sees the GPU through `singularity exec --rocm`, which no build-time +check can answer. Pull the `.sif`, put one real sample through it, and confirm +predict_3di lands on the GPU rather than silently falling back to CPU before +trusting this for production batches. diff --git a/docs/usage.md b/docs/usage.md index ee99c90..2b84d0a 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -18,14 +18,6 @@ Options: --use-conda / --no-use-conda Use conda for Snakemake rules [default: use-conda] --conda-prefix PATH Custom conda env directory - --container PATH container image with every per-rule tool - already installed (container/Dockerfile), - replacing --use-conda's per-rule env - creation for the WHOLE workflow -- not just - predict_3di (--prostt5-container). Needs - --use-singularity passed as a trailing - snakemake arg; don't also pass --use-conda - alongside this. --profile TEXT Snakemake profile --snake-default TEXT Customise Snakemake runtime args [default: --rerun-incomplete, --printshellcmds, @@ -64,13 +56,20 @@ Options: default) or prostt5-foldseek (structural; needs --hallmark-db, --hallmark-categories and --prostt5-checkpoint) [default: mmseqs] - --gpu-backend [cpu|cuda|rocm] PyTorch build for the ProstT5 conda env: + --gpu-backend [cpu|cuda|rocm|system] + PyTorch build for the ProstT5 conda env: cpu, cuda, or rocm (e.g. Setonix's MI250X nodes). Independent of --prostt5-cpu, which forces ProstT5 onto the CPU device at runtime even inside a GPU-capable env -- - this controls which env gets built - [default: cpu] + this controls which env gets built. 'system' + builds NO env at all and uses the torch + + pholdlib already installed in the ambient + python -- for when a working GPU torch is + already in place (e.g. inside the container, + or a module-loaded torch on HPC) and + installing a second one would be wasteful or + wrong [default: cpu] --foldseek-gpu use foldseek's CUDA GPU search mode for the hallmark scan (requires --gpu-backend cuda, a CUDA-capable foldseek build on PATH -- not @@ -112,13 +111,17 @@ Options: batch immediately [default: 4000] --prostt5-max-batch INTEGER max sequences per ProstT5 batch -- device- specific, tune per GPU [default: 20] - --prostt5-container TEXT container image with pholdlib + torch - already installed (e.g. phold's own image), - used instead of a conda env for predict_3di. - Needs --use-singularity passed as a trailing - snakemake arg -- --use-conda alone won't - honour it. Overrides --gpu-backend for this - rule. + --mfd-workers INTEGER number of worker processes for the flow- + decomposition step. Components are + independent, so they are split across + workers and merged in order (results are + identical to --mfd-workers 1, including + genome numbering). 1 = sequential, as + before. Note this is separate from + --threads: the MILP solver itself gets no + measurable benefit from extra threads, so + parallelism has to come from running + components concurrently [default: 1] --evalue FLOAT maximum e-value for phrog annotations [default: 1e-10] --seqidentity FLOAT minimum sequence identity for phrog @@ -132,26 +135,26 @@ Options: --prefix TEXT prefix for genome identifier -h, --help Show this message and exit. - +  If you use Phables in your work, please cite Phables as, - +  Vijini Mallawaarachchi, Michael J Roach, Przemyslaw Decewicz, Bhavya Papudeshi, Sarah K Giles, Susanna R Grigson, George Bouras, Ryan D Hesse, Laura K Inglis, Abbey L K Hutton, Elizabeth A Dinsdale, Robert A Edwards, Phables: from fragmented assemblies to high-quality bacteriophage genomes, Bioinformatics, Volume 39, Issue 10, October 2023, btad586, https://doi.org/10.1093/bioinformatics/btad586 - - +  +  For more information on Phables please visit: https://phables.readthedocs.io/ - - +  +  CLUSTER EXECUTION: phables run ... --profile [profile] For information on Snakemake profiles see: https://snakemake.readthedocs.io/en/stable/executing/cli.html#profiles - +  RUN EXAMPLES: Required: phables run --input [assembly graph file] Specify threads: phables run ... --threads [threads] @@ -190,7 +193,7 @@ Options: * `--databases` - path to the databases directory [default: wherever `phables install` put them] * `--use-conda` / `--no-use-conda` - use conda for Snakemake rules [default: `use-conda`] * `--conda-prefix` - custom conda env directory -* `--container` - run the whole workflow from a single container image instead of per-rule conda envs (needs `--use-singularity`, and `--no-use-conda`) -- see [Running from a single container](container.md) +* `--mfd-workers` - worker processes for the flow-decomposition (MFD) step [default: 1]. Components are independent, so they're split across processes and merged in component order — output is identical to `1`, genome numbering included. Distinct from `--threads`: the MILP solver gains nothing measurable from extra threads (building the model dominates, not solving it), so speedup has to come from running components concurrently. Expect ~2–2.5x on a realistic component mix — a few large components dominate the runtime and can't be split * `--snake-default` - customise Snakemake runtime args [default: `--rerun-incomplete, --printshellcmds, --nolock, --show-failed-logs`] ### Phage-gene detection: `--phagedetection` @@ -211,9 +214,10 @@ Two independent methods for finding phage-like genes on unitigs, feeding the sam * `--prostt5-half-precision` / `--prostt5-full-precision` - run ProstT5 in half precision (ignored on CPU) [default: `prostt5-half-precision`] * `--prostt5-cpu` - force ProstT5 onto CPU even when a GPU is available * `--prostt5-max-residues`, `--prostt5-max-seq-len`, `--prostt5-max-batch` - ProstT5 batching knobs, device-specific -- tune per GPU rather than trusting the defaults on unfamiliar hardware [defaults: 4000, 4000, 20] -* `--gpu-backend` - which PyTorch build the ProstT5 conda env solves against: `cpu`, `cuda`, or `rocm` (e.g. Setonix's MI250X). Independent of `--prostt5-cpu`, which forces the CPU device at runtime even inside a GPU-capable env -- this controls which env gets *built* [default: `cpu`] +* `--gpu-backend` - which PyTorch build the ProstT5 conda env solves against: `cpu`, `cuda`, or `rocm` (e.g. Setonix's MI250X). Independent of `--prostt5-cpu`, which forces the CPU device at runtime even inside a GPU-capable env -- this controls which env gets *built* [default: `cpu`]. A fourth value, `system`, builds **no** env at all and runs ProstT5 against the `torch`+`pholdlib` already present in the ambient python. Conda envs are isolated, so a rule that declares one can never see a torch installed outside it -- `system` is the only way to reuse an existing, known-good GPU torch (the [container](container.md)'s ROCm base image, or a module-loaded torch on HPC) instead of installing a second copy * `--foldseek-gpu` - use Foldseek's CUDA GPU search mode for the hallmark scan. Requires `--gpu-backend cuda`, a CUDA-capable Foldseek build on `PATH` (not the plain bioconda package), and a `*_gpu`-suffixed, `makepaddedseqdb`-prepared hallmark DB -* `--prostt5-container` - use a container with `pholdlib`+torch already installed (e.g. phold's own image) instead of a conda env for the ProstT5 step. Needs `--use-singularity` passed as a trailing Snakemake arg -- `--use-conda` alone won't honour it. Overrides `--gpu-backend` for this rule. + +To run all of this from one prebuilt image on HPC (Setonix), see [Running from a single container](container.md) — the image ships every conda env already built, so there are no per-rule container flags to pass. ## Example usage diff --git a/phables/__main__.py b/phables/__main__.py index 4b575ae..4a88e93 100644 --- a/phables/__main__.py +++ b/phables/__main__.py @@ -93,20 +93,6 @@ def common_options(func): type=click.Path(), show_default=False, ), - click.option( - "--container", - default=None, - required=False, - help=( - "container image with every per-rule tool already installed " - "(container/Dockerfile), replacing --use-conda's per-rule env " - "creation for the WHOLE workflow -- not just predict_3di " - "(--prostt5-container). Needs --use-singularity passed as a " - "trailing snakemake arg; don't also pass --use-conda alongside " - "this." - ), - type=click.Path(), - ), click.option( "--profile", help="Snakemake profile", default=None, show_default=False ), @@ -242,9 +228,14 @@ def run_options(func): "PyTorch build for the ProstT5 conda env: cpu, cuda, or rocm " "(e.g. Setonix's MI250X nodes). Independent of --prostt5-cpu, " "which forces ProstT5 onto the CPU device at runtime even " - "inside a GPU-capable env -- this controls which env gets built" + "inside a GPU-capable env -- this controls which env gets " + "built. 'system' builds NO env at all and uses the torch + " + "pholdlib already installed in the ambient python -- for when " + "a working GPU torch is already in place (e.g. inside the " + "container, or a module-loaded torch on HPC) and installing a " + "second one would be wasteful or wrong" ), - type=click.Choice(["cpu", "cuda", "rocm"]), + type=click.Choice(["cpu", "cuda", "rocm", "system"]), show_default=True, ), click.option( @@ -356,16 +347,20 @@ def run_options(func): show_default=True, ), click.option( - "--prostt5-container", - default=None, + "--mfd-workers", + default=1, required=False, help=( - "container image with pholdlib + torch already installed (e.g. " - "phold's own image), used instead of a conda env for predict_3di. " - "Needs --use-singularity passed as a trailing snakemake arg -- " - "--use-conda alone won't honour it. Overrides --gpu-backend for " - "this rule." + "number of worker processes for the flow-decomposition step. " + "Components are independent, so they are split across workers " + "and merged in order (results are identical to --mfd-workers 1, " + "including genome numbering). 1 = sequential, as before. Note " + "this is separate from --threads: the MILP solver itself gets " + "no measurable benefit from extra threads, so parallelism has " + "to come from running components concurrently" ), + type=int, + show_default=True, ), click.option( "--evalue", diff --git a/phables/config/config.yaml b/phables/config/config.yaml index 5d6e637..ee2d48c 100644 --- a/phables/config/config.yaml +++ b/phables/config/config.yaml @@ -57,35 +57,30 @@ prostt5_cpu: False # gpu_backend controls which env gets *built*. foldseek_gpu (below) is # separate again: foldseek's own --gpu mode is CUDA-only, so it stays off # regardless of gpu_backend unless explicitly requested and gpu_backend=cuda. +# +# A fourth value, `system`, means "don't build an env for predict_3di at all; +# torch and pholdlib are already installed in the ambient python". Conda envs +# are isolated, so a rule with a conda: directive can never see a torch that +# lives outside its env -- `system` is the only way to REUSE an existing, +# known-good GPU torch instead of installing a second one. That's exactly the +# situation inside container/Dockerfile, whose ROCm base image already ships a +# working torch, and it's also how you'd use a module-loaded torch on HPC. +# Note this key cannot be usefully overridden by editing a config file: the +# CLI hands every option to snaketool as merge_config=kwargs, so +# --gpu-backend's click default ("cpu") is merged OVER whatever is written +# here. Pass --gpu-backend on the command line instead. (Learned the hard way +# by baking `system` into the container's config and watching runs still use +# cpu -- see container/Dockerfile.) gpu_backend: cpu foldseek_gpu: False -# Optional: a single container image with EVERY per-rule tool this workflow -# needs already installed (e.g. "docker://quay.io/gbouras13/phables:" -- -# see container/Dockerfile). When set, every rule in the workflow runs with -# ONLY `container:` -- no `conda:` -- eliminating Snakemake's own per-rule -# --use-conda env creation entirely (and the real race condition that comes -# with it when many samples run as separate concurrent processes on an HPC -# array: multiple jobs needing the same not-yet-built shared env at once can -# corrupt each other's env-creation). Deliberately not combined with a -# `conda:` directive for the same reason as prostt5_container below: -# Snakemake's documented "conda inside container" feature builds a *separate* -# env inside the container rather than using what's already installed, which -# defeats the entire point here. When empty, every rule falls back to its own -# conda: env exactly as before. -# Needs `--use-singularity` (Apptainer/Singularity on most HPC, incl. Setonix) -# passed through phables' existing snake_args passthrough, since --use-conda -# alone won't honour the container: directive -- and --use-conda should -# generally NOT also be passed alongside this, for the reason above. -container: -# Optional: a container image (e.g. "docker://quay.io/gbouras13/phold:") -# that already has pholdlib + torch/ROCm baked in, used for JUST the -# predict_3di rule -- e.g. to point predict_3di at a different image than the -# rest of the workflow. Falls back to `container` above when unset, so -# setting `container` alone is enough to route predict_3di through it too; -# only set this separately if predict_3di specifically needs a different -# image. When both are empty, falls back to the conda envs selected by -# gpu_backend above. Same --use-singularity requirement as `container` above. -prostt5_container: +# Worker processes for the flow-decomposition step. Components are independent +# of each other (no component's logic reads another's results), so they are +# chunked across processes and merged in component order -- output is identical +# to mfd_workers: 1, genome numbering included. Separate from `threads`: the +# MILP solver gains nothing measurable from extra threads because building the +# model, not solving it, dominates -- so concurrency has to come from running +# components side by side instead. +mfd_workers: 1 # Batch limits are device-specific -- the defaults are conservative CPU/laptop-MPS # values, not tuned for any particular GPU. Retune per device (phold's autotune.py # is the reference approach) rather than trusting these at scale. diff --git a/phables/workflow/envs/prostt5-rocm.yaml b/phables/workflow/envs/prostt5-rocm.yaml index 6b33dc9..0361d7d 100644 --- a/phables/workflow/envs/prostt5-rocm.yaml +++ b/phables/workflow/envs/prostt5-rocm.yaml @@ -8,21 +8,33 @@ dependencies: # Unlike CUDA, ROCm is never the default pip torch build -- this index # override is required, not just for reproducibility. - --extra-index-url https://download.pytorch.org/whl/rocm6.3 - # torch==2.7.1+rocm6.3 pinned exactly, not just floored -- matches Setonix's - # own verified-working Pawsey container (quay.io/pawsey/pytorch:2.7.1-rocm6.3.3), - # confirmed real via https://download.pytorch.org/whl/rocm6.3/torch/ (the - # exact wheel exists for cp39-cp313). Switched to this conda-env path instead - # of the --prostt5-container route after repeatedly hitting container-specific - # problems on Setonix that had nothing to do with phables/pholdlib itself: - # a broken torchaudio needing libomp.so, --home overriding the container's - # apparent CWD so relative output paths silently resolved to the wrong place, - # and a missing MIOpen AI-heuristic kernel-selection file for gfx90a. All real, - # all diagnosed, none of them things a plain conda env run directly on the host - # (no container namespace, no bind-mount surface) is exposed to at all. + # torch==2.9.1+rocm6.3 pinned exactly, not just floored. Confirmed real on + # https://download.pytorch.org/whl/rocm6.3/torch/ (wheels for cp310-cp314, + # so every Python the `python>=3.10` pin above can resolve is covered -- + # 2.9.1 drops the cp39 wheel 2.7.1 had, which that pin already excludes). + # + # This env is the ONLY torch the workflow uses on ROCm, in both ways of + # running it. container/Dockerfile builds on Pawsey's ROCm base + # (quay.io/pawsey/pytorch:2.7.1-rocm6.3.3) for its ROCm userspace, but + # pre-builds THIS env inside the image and leaves the base's own + # system-python torch 2.7.1 untouched and unused -- so predict_3di runs on + # 2.9.1+rocm6.3 whether or not you're in the container. Keep it that way: + # the base image's torch version is an implementation detail of the base, + # not a version this workflow is pinned to. + # + # Phables drives ProstT5 through a plain conda env rather than pointing + # rules at a prebuilt torch container, after repeatedly hitting + # container-specific problems on Setonix that had nothing to do with + # phables/pholdlib itself: a broken torchaudio needing libomp.so, --home + # overriding the container's apparent CWD so relative output paths + # silently resolved to the wrong place, and a missing MIOpen AI-heuristic + # kernel-selection file for gfx90a. All real, all diagnosed, none of them + # things a conda env is exposed to. # >=2.6 would already avoid transformers' checkpoint loader refusing # torch.load() on older torch outright (CVE-2025-32434 hardening in # check_torch_load_is_safe, which broke loading Rostlab/ProstT5_fp16 -- no - # safetensors file); pinning to exactly 2.7.1 goes further, matching a - # version already confirmed to work on this exact hardware (MI250X, gfx90a). - - torch==2.7.1 + # safetensors file); pinning exactly goes further, keeping the version + # deterministic across every array task rather than drifting as new + # releases land mid-run. + - torch==2.9.1 - pholdlib>=0.1.2 diff --git a/phables/workflow/rules/02_phables_preflight.smk b/phables/workflow/rules/02_phables_preflight.smk index 4d59bc1..bfff53e 100644 --- a/phables/workflow/rules/02_phables_preflight.smk +++ b/phables/workflow/rules/02_phables_preflight.smk @@ -36,21 +36,7 @@ GC = config['genecaller'] PD = config['phagedetection'] GPU_BACKEND = config['gpu_backend'] FOLDSEEK_GPU = config['foldseek_gpu'] -# CONTAINER_IMAGE: workflow-wide container (e.g. a phables release image with -# every per-rule tool baked in), replacing --use-conda's per-rule env creation -# entirely when set -- every rule below does -# `conda: None if CONTAINER_IMAGE else os.path.join(...)` + -# `container: CONTAINER_IMAGE`, so a real value here disables conda for ALL of -# them at once (not just predict_3di). Run with --use-singularity, no -# --use-conda, once this is set. -# -# PROSTT5_CONTAINER falls back to CONTAINER_IMAGE when not set explicitly -- -# --prostt5-container remains available to point JUST predict_3di at a -# different image (e.g. phold's own) than the rest of the workflow, but -# --container alone is now sufficient to route everything, predict_3di -# included, through one image. -PROSTT5_CONTAINER = config['prostt5_container'] or config['container'] -CONTAINER_IMAGE = config['container'] +MFD_WORKERS = config['mfd_workers'] EV = config['evalue'] SI = config['seqidentity'] CT = config['covtol'] diff --git a/phables/workflow/rules/coverage.smk b/phables/workflow/rules/coverage.smk index 9d3908a..30be549 100644 --- a/phables/workflow/rules/coverage.smk +++ b/phables/workflow/rules/coverage.smk @@ -96,9 +96,7 @@ rule coverm_map: resources: mem_mb = config["resources"]["jobMem"] conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "coverm.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "coverm.yaml") log: os.path.join(LOGSDIR, "coverm_map.{sample}.log") shell: @@ -119,9 +117,7 @@ rule coverm_bam2counts: output: os.path.join(OUTDIR, "preprocess", "temp", "{sample}.cov") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "coverm.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "coverm.yaml") log: os.path.join(LOGSDIR, "coverm_bam2counts.{sample}.log") shell: diff --git a/phables/workflow/rules/genes.smk b/phables/workflow/rules/genes.smk index 2dce368..fdfd2f5 100644 --- a/phables/workflow/rules/genes.smk +++ b/phables/workflow/rules/genes.smk @@ -25,9 +25,7 @@ if GC == "pyrodigal-gv": log: os.path.join(LOGSDIR, "gene_call_pyrodigal_gv.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "genecall.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "genecall.yaml") script: os.path.join("..", "scripts", "gene_caller.py") @@ -48,9 +46,7 @@ else: out = os.path.join(LOGSDIR, "gene_call_fraggenescan_out.log"), err = os.path.join(LOGSDIR, "gene_call_fraggenescan_err.log"), conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "smg.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "smg.yaml") shell: """ run_FragGeneScan.pl -genome={input.genome} -out={params.frag} -complete=0 -train=complete -thread={threads} 1>{log.out} 2>{log.err} @@ -71,9 +67,7 @@ rule scan_smg: hmm_out=os.path.join(LOGSDIR, "smg_scan_hmm_out.log"), hmm_err=os.path.join(LOGSDIR, "smg_scan_hmm_err.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "smg.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "smg.yaml") shell: """ hmmsearch --domtblout {output.hmmout} --cut_tc --cpu {threads} {input.hmm} {input.faa} 1>{log.hmm_out} 2> {log.hmm_err} @@ -98,9 +92,7 @@ rule scan_phrogs: log: os.path.join(LOGSDIR, "phrogs_scan.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "mmseqs.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "mmseqs.yaml") shell: """ mkdir -p {params.out_path} @@ -135,74 +127,42 @@ if PD == "prostt5-foldseek": QUERY_3DI = os.path.join(OUTDIR, "preprocess", "hallmark", "proteins_3di.fasta") - if PROSTT5_CONTAINER: - - # container:-only, deliberately with no conda: alongside it. Snakemake's - # documented way to combine the two ("Ad-hoc combination of Conda package - # management with containers") builds a fresh, isolated conda env *inside* - # the container rather than exposing what the image already has installed - # -- which would just reinstall a second, redundant torch and ignore the - # image's own verified-working one (e.g. phold's own container, which - # already bundles pholdlib + a working torch/ROCm stack for Setonix, per - # the same logic phold's own Snakemake rules use). Needs - # `--use-singularity` passed through phables' snake_args passthrough -- - # --use-conda alone won't honour this directive. - rule predict_3di: - input: - faa = PROTEINS_FILE, - threads: - config["resources"]["jobCPU"] - resources: - mem_mb = config["resources"]["jobMem"] - output: - threedi = QUERY_3DI - params: - checkpoint = config["prostt5_checkpoint"], - model_name = config["prostt5_model"], - model_dir = config["prostt5_model_dir"], - half_precision = config["prostt5_half_precision"], - cpu = config["prostt5_cpu"], - max_residues = config["prostt5_max_residues"], - max_seq_len = config["prostt5_max_seq_len"], - max_batch = config["prostt5_max_batch"], - log: - os.path.join(LOGSDIR, "predict_3di.log") - container: - PROSTT5_CONTAINER - script: - os.path.join("..", "scripts", "predict_3di.py") - - else: - - rule predict_3di: - input: - faa = PROTEINS_FILE, - threads: - config["resources"]["jobCPU"] - resources: - mem_mb = config["resources"]["jobMem"] - output: - threedi = QUERY_3DI - params: - checkpoint = config["prostt5_checkpoint"], - model_name = config["prostt5_model"], - model_dir = config["prostt5_model_dir"], - half_precision = config["prostt5_half_precision"], - cpu = config["prostt5_cpu"], - max_residues = config["prostt5_max_residues"], - max_seq_len = config["prostt5_max_seq_len"], - max_batch = config["prostt5_max_batch"], - log: - os.path.join(LOGSDIR, "predict_3di.log") - conda: - # gpu_backend selects which torch build this env solves against -- - # cpu/cuda/rocm need different PyTorch wheels (conda envs are - # solved once from a static file, so this has to be three files, - # not one file with a runtime switch). See the individual env - # files for what each backend actually needs and why. - os.path.join("..", "envs", f"prostt5-{GPU_BACKEND}.yaml") - script: - os.path.join("..", "scripts", "predict_3di.py") + rule predict_3di: + input: + faa = PROTEINS_FILE, + threads: + config["resources"]["jobCPU"] + resources: + mem_mb = config["resources"]["jobMem"] + output: + threedi = QUERY_3DI + params: + checkpoint = config["prostt5_checkpoint"], + model_name = config["prostt5_model"], + model_dir = config["prostt5_model_dir"], + half_precision = config["prostt5_half_precision"], + cpu = config["prostt5_cpu"], + max_residues = config["prostt5_max_residues"], + max_seq_len = config["prostt5_max_seq_len"], + max_batch = config["prostt5_max_batch"], + log: + os.path.join(LOGSDIR, "predict_3di.log") + # gpu_backend selects which torch build this rule runs against. + # cpu/cuda/rocm each need a different PyTorch wheel, and a conda env is + # solved once from a static file, so those are three separate env files + # rather than one file with a runtime switch. + # + # `system` is the odd one out and takes NO conda: directive at all: + # conda envs are isolated, so a rule that declares one can never see a + # torch installed outside it. Omitting the directive is therefore the + # only way to REUSE an already-working torch (the container's ROCm base + # image, a module-loaded torch on HPC) instead of installing a second + # copy. The rule then runs in whichever python is running Snakemake, + # which must already provide torch + pholdlib. + conda: + None if GPU_BACKEND == "system" else os.path.join("..", "envs", f"prostt5-{GPU_BACKEND}.yaml") + script: + os.path.join("..", "scripts", "predict_3di.py") rule build_hallmark_query_db: @@ -216,9 +176,7 @@ if PD == "prostt5-foldseek": log: os.path.join(LOGSDIR, "build_hallmark_query_db.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "foldseek.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "foldseek.yaml") script: os.path.join("..", "scripts", "build_foldseek_query_db.py") @@ -271,9 +229,7 @@ if PD == "prostt5-foldseek": log: os.path.join(LOGSDIR, "scan_hallmark.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "foldseek.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "foldseek.yaml") shell: """ foldseek search {params.query_prefix} {input.hallmark_db} {params.result} {params.tmp} \ diff --git a/phables/workflow/rules/gfa2fasta.smk b/phables/workflow/rules/gfa2fasta.smk index 531aa2c..2e43c07 100644 --- a/phables/workflow/rules/gfa2fasta.smk +++ b/phables/workflow/rules/gfa2fasta.smk @@ -15,8 +15,6 @@ rule run_gfa2fasta: log: os.path.join(LOGSDIR, "gfa2fasta.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "phables.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "phables.yaml") script: os.path.join('..', 'scripts', 'gfa2fasta.py') \ No newline at end of file diff --git a/phables/workflow/rules/phables.smk b/phables/workflow/rules/phables.smk index 9a1ca07..fab47d3 100644 --- a/phables/workflow/rules/phables.smk +++ b/phables/workflow/rules/phables.smk @@ -36,14 +36,13 @@ rule run_phables: hallmark_minbits = config["hallmark_minbits"], output = os.path.join(OUTDIR, "phables"), nthreads = config["resources"]["jobCPU"], + mfd_workers = MFD_WORKERS, log = os.path.join(LOGSDIR, "phables_output.log") threads: config["resources"]["jobCPU"] log: os.path.join(LOGSDIR, "phables_output.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "phables.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "phables.yaml") script: os.path.join("..", "scripts", "phables.py") diff --git a/phables/workflow/rules/phylotree.smk b/phables/workflow/rules/phylotree.smk index 0f61961..caae230 100644 --- a/phables/workflow/rules/phylotree.smk +++ b/phables/workflow/rules/phylotree.smk @@ -11,9 +11,7 @@ rule build_msa: log: os.path.join(LOGSDIR, "mafft_output.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "phylotree.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "phylotree.yaml") shell: """ mafft --auto --thread {threads} {input} > {output} @@ -36,8 +34,6 @@ rule build_tree: log: os.path.join(LOGSDIR, "piqtree_output.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "phylotree.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "phylotree.yaml") script: os.path.join("..", "scripts", "phylotree.py") \ No newline at end of file diff --git a/phables/workflow/rules/postprocess.smk b/phables/workflow/rules/postprocess.smk index d481e5c..438f55a 100644 --- a/phables/workflow/rules/postprocess.smk +++ b/phables/workflow/rules/postprocess.smk @@ -75,9 +75,7 @@ rule coverm_map_genomes: resources: mem_mb = config["resources"]["jobMem"] conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "coverm.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "coverm.yaml") log: os.path.join(LOGSDIR, "coverm_map_genomes.{sample}.log") shell: @@ -90,22 +88,55 @@ rule coverm_map_genomes: rule coverm_bam2counts_genomes: - """Per-sample coverage stats over the resolved genomes.""" + """Per-sample coverage stats over the resolved genomes. + + Guarded against the zero-alignment case, which is a NORMAL outcome, not an + error: a sample where phables resolved no genomes AND had no unresolved + phage-like edges produces an empty genomes_and_unresolved_edges.fasta, so + coverm_map_genomes maps against an empty reference and emits a valid but + empty BAM. CoverM 0.7.0 panics outright on such a BAM rather than printing + an empty table -- + + [WARN coverm::contig] No primary alignments were observed for sample X + thread 'main' panicked at src/coverage_printer.rs:467:61: + index out of bounds: the len is 0 but the index is 0 + + -- which killed the whole run at the very last stage, after all the + expensive work had already succeeded. Everything downstream of here already + handles an empty table correctly (coverm_combine_genomes writes header-only + output; format_koverage_results.py's `readlines()[1:]` yields no rows and + pandas writes header-only report TSVs), so emitting the header ourselves is + all that's needed for the run to finish normally with empty report tables. + + The header below reproduces CoverM's own exactly -- " ", + space-joined, in the order the -m flags are given (confirmed against CoverM + 0.7.0 source: coverage_printer.rs writes "\\t{stoit_name} {estimator_header}", + mosdepth_genome_coverage_estimators.rs::column_headers defines the metric + strings, and bin/coverm.rs builds the estimator list by iterating the -m + flags in order). Note "Covered Fraction" is two space-separated words. + Stoit name is the BAM's basename without .bam, i.e. exactly {sample}. + """ input: os.path.join(OUTDIR, "postprocess", "temp", "{sample}.bam") output: temp(os.path.join(OUTDIR, "postprocess", "temp", "{sample}.cov")) conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "coverm.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "coverm.yaml") log: os.path.join(LOGSDIR, "coverm_bam2counts_genomes.{sample}.log") shell: """ - coverm contig -b {input} \ - -m count -m rpkm -m tpm -m mean -m covered_fraction -m variance \ - > {output} 2> {log} + n_aln=$(samtools view -c {input}) + if [ "$n_aln" -eq 0 ]; then + echo "No alignments in {input} -- no genomes were resolved for this sample (and no unresolved phage-like edges), so there is nothing to compute coverage over. Writing a header-only coverage table instead of running coverm, which panics on a zero-alignment BAM. The run continues and finishes normally; the per-genome report tables will be empty." > {log} + S={wildcards.sample} + printf 'Contig\\t%s Read Count\\t%s RPKM\\t%s TPM\\t%s Mean\\t%s Covered Fraction\\t%s Variance\\n' \ + "$S" "$S" "$S" "$S" "$S" "$S" > {output} + else + coverm contig -b {input} \ + -m count -m rpkm -m tpm -m mean -m covered_fraction -m variance \ + > {output} 2> {log} + fi """ @@ -151,8 +182,6 @@ rule format_genome_coverage: log: os.path.join(LOGSDIR, "format_koverage_results_output.log") conda: - None if CONTAINER_IMAGE else os.path.join("..", "envs", "phables.yaml") - container: - CONTAINER_IMAGE + os.path.join("..", "envs", "phables.yaml") script: os.path.join("..", "scripts", "format_koverage_results.py") \ No newline at end of file diff --git a/phables/workflow/scripts/format_koverage_results.py b/phables/workflow/scripts/format_koverage_results.py index e0f9c5e..9367c00 100644 --- a/phables/workflow/scripts/format_koverage_results.py +++ b/phables/workflow/scripts/format_koverage_results.py @@ -41,11 +41,19 @@ # / coverm_bam2counts_genomes / coverm_combine_genomes), completing PLAN.md # §4.8. The indices below therefore describe the CoverM-mode header: # -# Sample Contig Count RPKM TPM Mean Covered_fraction Variance -# 0 1 2 3 4 5 6 7 +# Sample Contig Read Count RPKM TPM Mean Covered Fraction Variance +# 0 1 2 3 4 5 6 7 # -# ("Covered_fraction" is coverm's `covered_fraction` method as it appears after -# coverm_combine_genomes strips the per-column " " prefix.) +# (Those are CoverM's own literal header strings as they appear after +# coverm_combine_genomes strips the per-column " " prefix -- verified +# against CoverM 0.7.0 source, mosdepth_genome_coverage_estimators.rs:: +# column_headers. The `count` method's header is "Read Count" and +# `covered_fraction`'s is "Covered Fraction" -- two space-separated words, NOT +# "Covered_fraction" as this comment previously claimed. Column ORDER follows +# the order the -m flags are passed in postprocess.smk, not any canonical +# order: bin/coverm.rs builds the estimator list by iterating the flags with +# .enumerate(). Only the indices below actually matter to this script; the +# names are documentation.) # # HISTORICAL, for anyone reading old output or an older checkout: this script # previously parsed Koverage's *native* "map" mode file, whose header is diff --git a/phables/workflow/scripts/phables.py b/phables/workflow/scripts/phables.py index 32913fe..b3b6830 100755 --- a/phables/workflow/scripts/phables.py +++ b/phables/workflow/scripts/phables.py @@ -50,6 +50,7 @@ def main(): hallmark_minbits = float(snakemake.params.hallmark_minbits) output = snakemake.params.output nthreads = int(snakemake.params.nthreads) + mfd_workers = int(snakemake.params.mfd_workers) log = snakemake.params.log # Setup logger @@ -95,6 +96,7 @@ def main(): logger.info(f"Input long reads: {longreads}") logger.info(f"Prefix for genome identifiers: {prefix}") logger.info(f"Number of threads to use: {nthreads}") + logger.info(f"Number of flow-decomposition workers: {mfd_workers}") logger.info(f"Output folder: {output}") if prefix is None or prefix == "": @@ -248,7 +250,7 @@ def main(): phage_like_edges, all_phage_like_edges, unresolved_phage_like_edges, - ) = short_utils.resolve_short( + ) = short_utils.resolve_short_parallel( assembly_graph, pruned_vs, unitig_names, @@ -269,6 +271,7 @@ def main(): prefix, output, nthreads, + mfd_workers, ) # Log final summary information diff --git a/phables/workflow/scripts/phables_utils/FD_Inexact.py b/phables/workflow/scripts/phables_utils/FD_Inexact.py index 09045c6..d250128 100644 --- a/phables/workflow/scripts/phables_utils/FD_Inexact.py +++ b/phables/workflow/scripts/phables_utils/FD_Inexact.py @@ -4,6 +4,7 @@ import itertools import logging +import math import flowpaths as fp import networkx as nx @@ -249,11 +250,81 @@ def flowMultipleDecomposition(data, K, nthreads): return data +def get_lowerbound_k(data): + """ + Smallest number of paths any valid decomposition of this component could + possibly use. + + The K search below tries K = 1, 2, 3, ... until one is feasible, rebuilding + the whole MILP each time (K is structural to the model -- every variable is + indexed by it -- so it genuinely cannot be reused, and neither flowpaths nor + HiGHS-via-flowpaths offers a warm start). Every attempt below the true answer + is therefore a complete model build that can only ever come back infeasible, + and model construction is ~95% of the cost of an attempt. Starting the search + at a proven lower bound skips exactly those wasted builds. + + Two bounds, both taken from flowpaths' own MinFlowDecomp.get_lowerbound_k: + + - the graph's WIDTH: the minimum number of paths needed to cover every edge. + Any decomposition must cover every edge carrying flow, so it cannot use + fewer paths than this. + - ceil(log2(number of distinct flow values)): k paths can produce at most + 2^k distinct subset sums, so k must be at least log2 of the number of + distinct values that have to be represented. + + Both are lower bounds, so the maximum of them is too, and starting there can + never skip a feasible smaller K -- the search still returns the same first + feasible K, just without the doomed attempts before it. Verified on synthetic + components: identical K and identical path sets with and without this. + + Returns 1 (i.e. the original behaviour) if anything about the computation + fails. A lower bound is an optimisation, not a correctness requirement, so it + must never be the reason a component stops resolving. + """ + try: + graph = data["graph"] + if graph.number_of_edges() == 0: + return 1 + + # flowpaths' stDAG wants string node ids, same relabelling + # flowMultipleDecomposition does before building its model. + node_labels = {node: str(node) for node in graph.nodes} + flowpaths_graph = nx.relabel_nodes(graph, node_labels, copy=True) + lowerbound = fp.stDAG(flowpaths_graph).get_width() + + # Only edges that must carry flow constrain the decomposition; an edge + # whose lower bound is 0 need not be covered at all. + distinct_flows = { + int(data["flows_low"][edge]) + for edge in data["edges"] + if data["flows_low"].get(edge, 0) > 0 + } + if len(distinct_flows) > 1: + lowerbound = max(lowerbound, math.ceil(math.log2(len(distinct_flows)))) + + return max(1, lowerbound) + except Exception as e: + logger.debug(f"Could not compute a lower bound on K ({e}); starting from 1") + return 1 + + def FD_Algorithm(data, max_paths, nthreads): solutionWeights = 0 solutionSet = 0 - for i in range(1, max_paths + 1): + # See get_lowerbound_k: K < lowerbound is provably infeasible, and each such + # attempt costs a full MILP build. Capped at max_paths so a bound above the + # user's --maxpaths doesn't turn into a search over an empty range with + # different semantics -- that case has no solution within max_paths either + # way, and this keeps the "attempt at least one K" behaviour identical. + lowerbound = min(get_lowerbound_k(data), max_paths) + if lowerbound > 1: + logger.debug( + f"Starting the K search at {lowerbound} rather than 1 " + f"({lowerbound - 1} provably-infeasible attempt(s) skipped)" + ) + + for i in range(lowerbound, max_paths + 1): data = flowMultipleDecomposition(data, i, nthreads) if data["message"] == "solved": solutionSet = data["solution"] @@ -348,6 +419,11 @@ def SolveInstances(Graphs, max_paths, outfile, recfile, nthreads): "adj_in": AD_in, "adj_out": AD_out, "subpaths": Graphs[s]["subpaths"], + # Never read by anything -- the K search takes its starting point + # from get_lowerbound_k(), which computes a real bound per component + # rather than assuming a constant. Left in place because this dict is + # passed around wholesale and removing a key is a wider change than + # it looks; flagged so it isn't mistaken for the live lower bound. "minK": 2, "runtime": 0, } diff --git a/phables/workflow/scripts/phables_utils/edge_graph_utils.py b/phables/workflow/scripts/phables_utils/edge_graph_utils.py index 41c3b00..1e0667d 100644 --- a/phables/workflow/scripts/phables_utils/edge_graph_utils.py +++ b/phables/workflow/scripts/phables_utils/edge_graph_utils.py @@ -59,6 +59,19 @@ def get_unitig_lengths(edge_file): return unitig_lengths +def _oriented_links_inner(): + """Inner factory for oriented_links. + + A module-level function rather than the `lambda: defaultdict(list)` this + used to be, purely so the resulting structure can be PICKLED. Lambdas + cannot be, which meant oriented_links could not cross a process boundary -- + it broke resolve_short_parallel's worker pool with + `PicklingError: Can't pickle >` the first time it ran on + real data. Behaviour is identical: missing keys still get a defaultdict(list). + """ + return defaultdict(list) + + def get_links(assembly_graph_file): """ Get links from the assembly graph @@ -67,7 +80,7 @@ def get_links(assembly_graph_file): node_count = 0 graph_contigs = {} edges_lengths = {} - oriented_links = defaultdict(lambda: defaultdict(list)) + oriented_links = defaultdict(_oriented_links_inner) link_overlap = defaultdict(int) links = [] diff --git a/phables/workflow/scripts/phables_utils/short_utils.py b/phables/workflow/scripts/phables_utils/short_utils.py index 9fb760c..b9ee0f8 100644 --- a/phables/workflow/scripts/phables_utils/short_utils.py +++ b/phables/workflow/scripts/phables_utils/short_utils.py @@ -1,5 +1,9 @@ #!/usr/bin/env python3 +import math +import multiprocessing +import pickle +from concurrent.futures import ProcessPoolExecutor import logging import sys import time @@ -16,6 +20,12 @@ # Create logger logger = logging.getLogger("phables 2.0.0") +# How many chunks to create per worker in resolve_short_parallel. >1 so the +# pool can load-balance an uneven component workload; small enough that +# per-chunk overhead stays negligible against per-component solve times +# measured in tens of milliseconds. +CHUNKS_PER_WORKER = 4 + def resolve_short( assembly_graph, @@ -1498,3 +1508,208 @@ def resolve_short( all_phage_like_edges, unresolved_phage_like_edges, ) + + +# Set in the PARENT before the pool is created; forked workers inherit it +# through the process image, so none of it is ever serialised. +# +# This was previously passed as ProcessPoolExecutor(initargs=(kwargs,)), which +# pickles the whole thing once per worker. That is fine for a small graph and +# fatal for a real one: on a 485k-vertex assembly the payload is the igraph +# object plus every unitig's sequence in graph_unitigs -- gigabytes, serialised +# eight times through a pipe. It killed the pool seconds after startup, and +# multiprocessing additionally cannot send a single object larger than ~2GB at +# all. Inheriting by fork moves that cost to zero, and copy-on-write means the +# eight workers share the pages rather than each holding a full copy. +_WORKER_KWARGS = None + + +def _resolve_short_chunk(chunk): + """Worker entry point: one slice of components, nothing else. + + Module-level (not a closure) so it is picklable by ProcessPoolExecutor -- + though only the small `chunk` is ever pickled; the bulk inputs arrive via + _WORKER_KWARGS, inherited from the parent. + """ + return resolve_short(pruned_vs=chunk, **_WORKER_KWARGS) + + +def resolve_short_parallel( + assembly_graph, + pruned_vs, + unitig_names, + unitig_names_rev, + self_looped_nodes, + graph_unitigs, + minlength, + link_overlap, + unitig_coverages, + compcount, + oriented_links, + junction_pe_coverage, + likely_complete, + alpha, + mincov, + covtol, + maxpaths, + prefix, + output, + nthreads, + workers=1, +): + """ + resolve_short, with components processed in parallel across processes. + + Deliberately implemented WITHOUT touching resolve_short's ~1400-line body: + that function is already parameterised by the set of components to process + (`pruned_vs`) and already returns every accumulator it builds, so splitting + the components into contiguous chunks, running one chunk per worker, and + merging the returned tuples is equivalent to running it once over all of + them. The alternative -- extracting the loop body into a per-component + function -- would mean restructuring 1400 lines of nested branching for no + additional benefit. + + This is only sound because no component's logic depends on another's + results. Verified against the body: every touch of a shared accumulator is + a pure `X.add(...)` / `X = X.union(...)` / `X.append(...)`, and there is not + one conditional or membership test against them anywhere in the loop -- + per-component decisions use the loop-local `comp_*` sets instead. If that + ever stops being true, chunking silently changes results, so it is worth + re-checking before extending this. + + Chunks are CONTIGUOUS and merged in chunk order, so `all_resolved_paths` + ends up in exactly the same order as the sequential run. That matters + because downstream naming numbers the genomes by position: a different + order would rename every genome without changing the biology, which is a + nasty kind of non-reproducibility. + + Falls back to a plain sequential call for workers <= 1 or a single chunk, + so the default path is byte-for-byte the code that ran before. + """ + kwargs = dict( + assembly_graph=assembly_graph, + unitig_names=unitig_names, + unitig_names_rev=unitig_names_rev, + self_looped_nodes=self_looped_nodes, + graph_unitigs=graph_unitigs, + minlength=minlength, + link_overlap=link_overlap, + unitig_coverages=unitig_coverages, + compcount=compcount, + oriented_links=oriented_links, + junction_pe_coverage=junction_pe_coverage, + likely_complete=likely_complete, + alpha=alpha, + mincov=mincov, + covtol=covtol, + maxpaths=maxpaths, + prefix=prefix, + output=output, + # Each worker is its own process, so the MILP solver inside it should + # not also try to use every core -- that would oversubscribe by + # workers x nthreads. Profiling showed solver threads make no + # measurable difference anyway (model construction dominates), so 1 is + # the right per-worker value rather than a compromise. + nthreads=1, + ) + + keys = list(pruned_vs) + if workers <= 1 or len(keys) <= 1: + return resolve_short(pruned_vs=pruned_vs, **{**kwargs, "nthreads": nthreads}) + + # Everything below has to cross a process boundary, so it all has to pickle. + # Checked up front rather than discovered when the pool starts: a failure + # there aborts the whole phables run, and losing a completed assembly to a + # performance optimisation is a bad trade. Hit for real -- oriented_links + # was a defaultdict built with a lambda, which cannot be pickled (fixed at + # source in edge_graph_utils._oriented_links_inner), and it took out a real + # 120-component run. Falling back to the sequential path keeps that a slow + # run instead of a failed one. + try: + pickle.dumps(kwargs) + except Exception as e: + logger.warning( + f"Cannot run components in parallel -- some input is not picklable " + f"({type(e).__name__}: {e}). Falling back to sequential; the result " + f"is unaffected, only the runtime." + ) + return resolve_short(pruned_vs=pruned_vs, **{**kwargs, "nthreads": nthreads}) + + workers = min(workers, len(keys)) + + # Workers must INHERIT the bulk inputs rather than be sent them. On a real + # assembly those inputs are gigabytes (see _WORKER_KWARGS), and pickling + # them per worker is both ruinously slow and subject to multiprocessing's + # ~2GB per-object ceiling. Fork gives the children the parent's memory + # image directly, at no serialisation cost and, thanks to copy-on-write, + # very little extra memory. + # + # If fork is unavailable (Windows; macOS defaults to spawn but can still + # fork explicitly) there is no cheap way to hand over that much data, so + # this runs sequentially rather than attempting a copy that would either + # fail outright or exhaust the node's memory. + try: + mp_context = multiprocessing.get_context("fork") + except ValueError: + logger.warning( + "The 'fork' start method is unavailable, so component-level " + "parallelism would have to copy the whole assembly graph to every " + "worker. Running sequentially instead; the result is unaffected, " + "only the runtime." + ) + return resolve_short(pruned_vs=pruned_vs, **{**kwargs, "nthreads": nthreads}) + + # Deliberately MORE chunks than workers. Component cost is heavily skewed -- + # a handful of large components dominate, and which ones is not known in + # advance -- so splitting into exactly one chunk per worker (static + # partitioning) lets a single worker draw several expensive components while + # the rest sit idle. Smaller chunks let the pool hand out more work to + # whichever process finishes first, and cost nothing extra to send, since + # only the chunk itself crosses the boundary. + size = max(1, math.ceil(len(keys) / (workers * CHUNKS_PER_WORKER))) + chunks = [ + {k: pruned_vs[k] for k in keys[i : i + size]} for i in range(0, len(keys), size) + ] + + logger.info( + f"Resolving {len(keys)} components across {workers} worker process(es) " + f"in {len(chunks)} chunk(s)" + ) + + # Published to the module namespace BEFORE the pool exists, so every forked + # child sees it. Cleared afterwards so the parent does not keep a second + # reference to the graph alive for the rest of the run. + global _WORKER_KWARGS + _WORKER_KWARGS = kwargs + try: + with ProcessPoolExecutor( + max_workers=workers, mp_context=mp_context + ) as executor: + # .map preserves input order, which is what keeps the merged + # all_resolved_paths identical to the sequential run. + results = list(executor.map(_resolve_short_chunk, chunks)) + except Exception as e: + # A worker dying (OOM above all -- eight processes touching a large + # graph can outgrow the node) surfaces here as BrokenProcessPool. Better + # to spend the extra wall-clock than to lose an assembly that already + # cost hours, so this retries the whole thing sequentially. + logger.warning( + f"Parallel component resolution failed ({type(e).__name__}: {e}). " + f"Falling back to sequential -- the result is unaffected, only the " + f"runtime. If this is memory, lower --mfd-workers." + ) + return resolve_short(pruned_vs=pruned_vs, **{**kwargs, "nthreads": nthreads}) + finally: + _WORKER_KWARGS = None + + # Field 1 (all_resolved_paths) and field 2 (all_components) are lists and + # are extended; every other field is a set and is unioned. Merged in chunk + # order for the reason given above. + merged = list(results[0]) + for result in results[1:]: + for i, value in enumerate(result): + if isinstance(value, list): + merged[i] = merged[i] + value + else: + merged[i] = merged[i] | value + return tuple(merged)