Feature/prostt5 foldseek detection by gbouras13 · Pull Request #69 · Vini2/phables

gbouras13 · 2026-08-20T23:56:28Z

fixes 3 issues:

crashing when no genomes are resolved that my previous PR introduced (instead it has a gracious exit like previous)
Multi-threading and huge efficiency gains for the flow decomposition step
Removal of all previous container stuff and replacement of one monolithic Setonix container ala Hybracter (way easier to deal with, you just need to download the sif file )

I have tested widely on setonix.

Claude in more detail below

George

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-runs 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.

… panic) A sample where phables resolves no genomes and has no unresolved phage-like edges produces an empty genomes_and_unresolved_edges.fasta, so coverm_map_genomes emits a valid but empty BAM. CoverM 0.7.0 panics on such a BAM ('index out of bounds: the len is 0 but the index is 0', coverage_printer.rs:467) instead of printing an empty table, killing the whole run at the very last stage after all the expensive work had already succeeded. Guard the rule: if the BAM has no alignments, write coverm's header with no data rows rather than invoking coverm. Everything downstream already handles an empty table (verified end-to-end: 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 the run now completes normally with empty report tables. The fallback header reproduces CoverM's own exactly -- verified against CoverM 0.7.0 source rather than guessed: coverage_printer.rs joins '<stoit> <metric>' with a space, 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. That also disproved an existing comment in format_koverage_results.py claiming the covered_fraction column reads 'Covered_fraction' -- it is actually 'Covered Fraction', two space-separated words, and count is 'Read Count'. Corrected; the IDX_* indices it documents were already right.

torch==2.9.1+rocm6.3 confirmed to exist on the pinned ROCm 6.3 index (download.pytorch.org/whl/rocm6.3/torch/), with cp310-cp314 wheels -- so every Python the env's own python>=3.10 pin can resolve is covered (2.9.1 drops the cp39 wheel 2.7.1 shipped, which that pin already excludes anyway). The --extra-index-url is unchanged; rocm6.3 still carries this version. This intentionally diverges the conda-env ROCm path from the container path, which still builds on quay.io/pawsey/pytorch:2.7.1-rocm6.3.3 and asserts torch 2.7.1 -- noted inline in the yaml so the mismatch reads as deliberate rather than an oversight.

Removes --container and --prostt5-container entirely: both CLI options, both config keys, the CONTAINER_IMAGE/PROSTT5_CONTAINER preflight vars, all 15 per-rule 'container:' directives, and genes.smk's PROSTT5_CONTAINER if/else split around predict_3di. Every rule is back to a plain 'conda:' directive. In their place, container/Dockerfile is now one monolithic image built the way hybracter's is: Pawsey's ROCm base for the gfx90a userspace, miniforge installed inside, and EVERY per-rule conda env pre-built into the image by container/prebuild_envs.sh. Convert it to a single .sif and 'phables run' inside it finds its envs already there -- which is what actually kills the conda-env race condition (many array tasks each being their own Snakemake process, all creating the same shared env at once), since at runtime there is nothing left to create. Envs are built at phables' DEFAULT --conda-prefix on purpose: Snakemake hashes an env file's content together with the prefix path, so building and running at the same prefix in the same image makes the hashes match and nothing is rebuilt -- important because a .sif is read-only, so a rebuild attempt is a hard failure rather than a slow path. Databases stay OUT of the image. prebuild_envs.sh creates zero-byte placeholder DB files purely so the DAG can resolve, then runs --conda-create-envs-only once per flag combination (gene caller, detection mode, GPU backend, tree, install). Verified for real: the DAG resolves against zero-byte placeholder databases, and all six flag combinations dry-run clean after the container removal. container/test_image.sh smoke-tests the built image and fails the build if any pre-built env is missing a binary the rules actually invoke, or if torch didn't come from the prostt5-rocm env rather than the base image's own unused 2.7.1. Not yet built: no docker build or Setonix Apptainer run has happened. The pinned conda packages were confirmed present for linux-64/noarch (mmseqs2=13.45111, foldseek, coverm, fraggenescan, cogent3<2026.7 on bioconda/noarch; piqtree is pip-only, which phylotree.yaml already reflects).

1. tests/data/ is gitignored, so it doesn't exist in a fresh clone or CI checkout -- the build died at the pre-build step with 'Invalid value for --reads: Path ... does not exist'. prebuild_envs.sh now generates its own throwaway inputs (a two-segment GFA and a pair of tiny gzipped FASTQs) instead of depending on untracked data. Verified these produce an identical DAG to the real test data. 2. 'phables install' was pointed at the placeholder databases dir, whose files satisfy install.smk's own download targets -- so Snakemake said 'Nothing to be done', the DAG was empty, and the curl env would have been silently missing from the image (with test_image.sh's own 'check_bin curl' the only thing that would have caught it). It now uses an empty databases dir so the four *_download rules are in the DAG and their env actually gets built. prebuild_envs.sh has now been run end-to-end with its --conda-create-envs-only calls swapped for dry-runs: all six invocations succeed, the synthetic inputs resolve, and cleanup fires. The conda solves themselves and test_image.sh against a real image remain unverified.

Adds --gpu-backend system: predict_3di declares NO conda env and runs against the torch + pholdlib already in the ambient python. Conda envs are isolated, so a rule declaring conda: can never see a torch installed outside its env -- omitting the directive is the only way to reuse an existing, known-good GPU torch. Verified via --list-conda-envs: system declares 5 envs, rocm declares 6 (the extra being prostt5-rocm). The container now uses it, which required inverting how the image is built: phables, Snakemake and pholdlib install into the BASE IMAGE'S python (the one holding Pawsey's verified ROCm torch), not miniforge's. Snakemake runs a script: rule with no conda env using its own interpreter, so that interpreter has to be the one with torch. Miniforge is now APPENDED to PATH, never prepended, so it cannot shadow the base python -- prepending it is exactly why the previous build ran phables under a torch-less python and needed its own torch env. Consequences: no prostt5-* env is built at all (the prostt5-cpu build is gone too), the image drops several GB and one whole conda solve, and it runs the torch already verified on gfx90a rather than an unverified 2.9.1. The image sets its own config default to gpu_backend: system so a plain works with no extra flag -- without that, the config default (cpu) would look for an env the image deliberately lacks, which is fatal on a read-only .sif. The build asserts torch's version is unchanged across the pip install, and test_image.sh now fails if any pre-built env contains a torch of its own -- i.e. if a second copy got installed after all. prostt5-rocm.yaml is untouched and still pins 2.9.1 for non-container use.

The first full image build succeeded through every conda solve and then failed at the very last step on my own over-strict check: it tested for the substring 'rocm' in torch.__version__, but Pawsey's base torch is a source build reporting '2.7.1a0+gite2d141d' with no +rocm6.3 suffix, so a perfectly good ROCm torch was rejected. Now prints torch.__version__, torch.version.hip and torch.version.cuda and continues, warning only if hip is None. What that torch is compiled against is Pawsey's business, and failing a multi-GB build at the last step to re-litigate it is a bad trade. The checks that actually matter stay fatal: torch and pholdlib must be importable in the ambient python (predict_3di cannot run otherwise), and no pre-built conda env may contain a torch of its own (which would mean a second copy was installed after all). That build otherwise confirmed the design end to end: all 8 per-rule envs solve and install, no prostt5-* env is created, and the pip install into the base python leaves its torch untouched.

The first real container run on Setonix failed: Creating conda environment .../envs/prostt5-cpu.yaml... OSError: [Errno 30] Read-only file system: '/usr/local/.../phables/workflow/conda/4a0d2f6693d2fae2ac9de8108a44f6cc_.yaml' The image's config.yaml did say 'gpu_backend: system', but that is not the effective value: 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 default is 'cpu'. So predict_3di asked for a prostt5-cpu env the image deliberately does not contain, on a read-only filesystem. Reproduced locally before fixing: with config.yaml patched to 'system' and no flag passed, the runtime config still reports gpu_backend: cpu; passing --gpu-backend system reports system. Removes the Dockerfile's sed step (it never worked) and replaces test_image.sh's config-file grep, which asserted the wrong thing and passed while real runs used cpu. It now checks what the image can actually guarantee -- that the 'system' choice exists -- and the docs state that callers must pass the flag.

FD_Algorithm tried K = 1, 2, 3, ... until a feasible decomposition was found, rebuilding the whole MILP for each K. K is structural to the model (every variable is indexed by it), so the model genuinely cannot be reused, and flowpaths does not expose a HiGHS warm start -- so every attempt below the true answer was a full model build that could only ever return infeasible. Profiling showed why that dominates: model CONSTRUCTION is ~95% of the cost of an attempt, not solving (a large component: 170ms build vs 8.9ms solve). That also explains why --threads never helped -- threads only affect the 5% -- and why a solver time_limit does nothing here, both of which were measured and rejected before landing on this. get_lowerbound_k() takes the maximum of two bounds, both from flowpaths' own MinFlowDecomp.get_lowerbound_k: the graph width (minimum paths needed to cover every flow-carrying edge) and ceil(log2(#distinct flow values)). Both are lower bounds, so starting there cannot skip a feasible smaller K. It costs 1-4ms 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. Measured on synthetic components, against the previous behaviour through the real SolveInstances entry point: - 18/18 cases returned IDENTICAL K, path count and path sets - 1.5x-4.9x faster, growing with component size - components that cannot resolve within --maxpaths: up to 5.9x, since the bound proves K >= maxpaths up front instead of burning the whole ladder to learn nothing - edge cases verified identical: single edge, all-zero lower bounds, one distinct flow value, bound above --maxpaths, unresolvable Also annotates data["minK"], which was set to a constant 2 and never read by anything, so it isn't mistaken for the live lower bound.

Components are independent, so the flow-decomposition loop is embarrassingly parallel. resolve_short_parallel chunks the components, runs the EXISTING resolve_short once per chunk in a worker process, and merges the returned accumulators. Deliberately no change to resolve_short's ~1400-line body: it is already parameterised by the component set (pruned_vs) and already returns every accumulator it builds, so chunk-and-merge is equivalent to one call over all components. Extracting the loop body into a per-component function would have meant restructuring 1400 lines of nested branching for the same result. Soundness rests on no component's logic depending on another's results, which was verified against the body before writing any of this: every touch of a shared accumulator is a pure add/union/append, and there is not one conditional or membership test against them in the loop -- per-component decisions use the loop-local comp_* sets. Noted in the docstring, since chunking would silently change results if that stops being true. Chunks are merged in component order, so all_resolved_paths ends up identical to the sequential run. That matters: genomes are numbered by position, so a different order would rename every genome without changing the biology. Two implementation details that turned out to matter: - MORE chunks than workers (CHUNKS_PER_WORKER=4). Component cost is heavily skewed and which components are expensive is not known in advance, so one chunk per worker leaves workers idle. On a realistic skewed workload this took 8 workers from 1.92x to 2.35x -- 98% of the 2.4x Amdahl ceiling imposed by the single largest component. - Heavy read-only inputs (the assembly graph above all) are sent once per worker via a pool initializer, not once per chunk, so smaller chunks do not mean repeatedly re-pickling the graph. Per-worker solver threads are pinned to 1: workers x nthreads would oversubscribe, and profiling showed solver threads make no measurable difference anyway. Measured against sequential, through the real wrapper: - identical 18-tuple results at 2, 4 and 8 workers - all_resolved_paths order identical (genome numbering preserved) - uniform workload: 6.1x at 8 workers - realistic skewed workload: 2.35x at 8 workers (ceiling 2.4x) - edge cases identical: 1 component, fewer components than workers, more components than workers, and workers=1 Default is 1, i.e. the sequential path is unchanged.

The comment still described the image as pre-building two torch envs (rocm + cpu). That stopped being true when predict_3di switched to reusing the base image's ROCm torch via --gpu-backend system -- no prostt5-* env is built at all now.

The first real run with --mfd-workers 8 (120 components, SRR27716023) died immediately after the pool started. Cause: edge_graph_utils built oriented_links as defaultdict(lambda: defaultdict(list)), and lambdas cannot be pickled -- so the structure could not cross a process boundary and every worker failed at startup. The sequential path never pickles anything, which is why this only appeared with workers > 1. Fixed at source: the inner factory is now a module-level function, which pickles. Behaviour is identical -- missing keys still get a defaultdict(list) at both levels, verified across a pickle round-trip. Also adds a preflight in resolve_short_parallel: the kwargs are test-pickled before the pool is created, and if anything is unpicklable it logs a clear warning and runs sequentially instead. An optimisation should never be able to destroy a completed assembly's run, and this failure mode is invisible until it happens on real data -- same reasoning as get_lowerbound_k falling back to 1 on error. Verified: with the fix the pool runs and merges all components; with the old lambda form the preflight catches it, warns, falls back, and returns results identical to the parallel path.

…pickle --mfd-workers 8 still died on a real sample (SRR12983552: 485,547 vertices, 2733 components), seconds after the pool started, even with the oriented_links pickling fix. Cause: ProcessPoolExecutor(initargs=(kwargs,)) serialises the ENTIRE input set once per worker -- the igraph object plus every unitig sequence in graph_unitigs. On a real assembly that is gigabytes, sent eight times through a pipe, and multiprocessing cannot transfer a single object larger than ~2GB at all. It worked on the 22k-vertex sample and fell off that cliff on the 485k-vertex one. Workers now INHERIT the inputs: _WORKER_KWARGS is set in the parent before the pool is created, and the pool uses an explicit 'fork' context, so children get the parent's memory image with no serialisation and, via copy-on-write, without eight full copies. Measured with a 200MB stand-in payload: four workers receive it intact in 0.08s, where pickling it per worker costs ~0.5s and ~1.6GB copied for eight -- and the real payload is far larger. Two safety nets, because this failure destroyed a run that had already cost ~5 hours of GPU time: - if 'fork' is unavailable, run sequentially rather than attempt a copy that would fail or exhaust the node - if the pool breaks anyway (a worker OOM-killed surfaces as BrokenProcessPool), warn and retry sequentially instead of failing the run; verified by SIGKILLing a worker mid-flight _WORKER_KWARGS is cleared afterwards so the parent does not hold a second reference to the graph for the rest of the run.

codacy-production · 2026-08-20T23:57:35Z

Not up to standards ⛔

🔴 Issues 1 critical · 1 high · 1 medium

Alerts:
⚠ 3 issues (≤ 0 issues of at least minor severity)

Results:
3 new issues

Category Results

ErrorProne 1 high

Security 1 critical
1 medium

View in Codacy

🟢 Metrics 20 complexity · 0 duplication

Metric Results

Complexity 20

Duplication 0

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer}
_{TIP This summary will be updated as you push new changes.}

gbouras13 added 12 commits August 13, 2026 23:02

Navigation Menu

Feature/prostt5 foldseek detection - #69

Feature/prostt5 foldseek detection#69
gbouras13 wants to merge 12 commits into
Vini2:developfrom
gbouras13:feature/prostt5-foldseek-detection

gbouras13 commented Aug 20, 2026

Uh oh!

codacy-production Bot commented Aug 20, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

gbouras13 commented Aug 20, 2026

1. Fix: phables crashes when a sample resolves no genomes

2. Performance: flow decomposition

3. Container: one monolithic image, replacing the per-rule container flags

Testing

Not verified

Uh oh!

codacy-production Bot commented Aug 20, 2026

Not up to standards ⛔

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant