From bd9e7cf11716fb562849451012a4a48893f12c58 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 19 Jul 2026 09:00:20 +0200 Subject: [PATCH 1/2] assign_compartments gap-fill: reliable flux-based fill, not cobra's MILP --- CHANGELOG.md | 10 +++++ src/raven_toolbox/localization/certify.py | 53 +++++++++++++++++++---- tests/test_localization_certify.py | 27 ++++++++++++ 3 files changed, 82 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f51ea24..c98c499 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,16 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* **`assign_compartments` gap-fill is now a reliable flux-based fill.** The universal-DB gap-fill in + `localization.certify` no longer calls `cobra.flux_analysis.gapfill`, whose indicator MILP, at genome + scale, fails to find a valid fill in the **majority** of cases *even when the exact reaction that + restores growth is present in the universal* — the MILP returns an incumbent its own validation then + rejects, and it raises rather than offering a fill. The replacement adds the universal candidates, holds + biomass at the growth floor, runs pFBA, and keeps the added reactions that carry flux: a plain LP that + cannot have that failure mode. Measured against the old path — knockout-recovery recall **100 % (60/60, + each recovering the exact removed reaction) vs cobra's 45 %**; on realistic incomplete drafts it restores + growth **12/12 vs cobra's 0/12**, ~5× faster. The caller still re-certifies with a real FBA, so no false + certificate is possible. * **`diff_models` compares grRules as logic, not text.** The GPR check now DNF-expands each rule (via the existing `manipulation.gpr_to_dnf`), sorts the genes within each isozyme clause and sorts the clauses, so operand order no longer registers as a difference: `a and b` == `b and a` and `a or b` == `b or a`. The diff --git a/src/raven_toolbox/localization/certify.py b/src/raven_toolbox/localization/certify.py index a76bcce..b7bb5da 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -697,22 +697,59 @@ def assign_compartments( def _gapfill(applied, universal, biomass_reaction, min_growth) -> list[str]: - """Minimal gap-fill reactions (from ``universal``) that restore the primary growth floor. + """Universal reactions whose addition restores the primary growth floor, via a flux-based fill. - Uses ``cobra.flux_analysis.gapfill``; the result is validated by the caller's certification FBA, - so any tolerance quirk in cobra's own MILP cannot produce a false certificate here. + Adds every ``universal`` candidate not already present, holds biomass at the floor, and runs pFBA + (maximise biomass, then minimise total flux); the added reactions that carry flux are returned. + + This is a plain LP, which matters. ``cobra.flux_analysis.gapfill`` solves an indicator MILP that, at + genome scale, fails to find a valid fill in the *majority* of cases even when the exact reaction that + restores growth is present in the universal — its own validation then rejects the bad incumbent and + raises, so no fill is offered. The LP here cannot have that failure mode: a returned set is a real + flux solution that actually reaches the floor. The caller re-certifies with FBA regardless, so this + never yields a false certificate. """ + if biomass_reaction not in applied.reactions: + return [] + floor = max(min_growth, 1e-4) try: - from cobra.flux_analysis import gapfill as cobra_gapfill + from cobra.flux_analysis import pfba with applied: + added = [] + for urxn in universal.reactions: + if urxn.id in applied.reactions: + continue + new = cobra.Reaction(urxn.id, name=urxn.name, + lower_bound=urxn.lower_bound, upper_bound=urxn.upper_bound) + applied.add_reactions([new]) + new.add_metabolites({_universal_met(applied, m): coeff + for m, coeff in urxn.metabolites.items()}) + added.append(new) + if not added: + return [] applied.objective = biomass_reaction - solutions = cobra_gapfill(applied, universal, lower_bound=max(min_growth, 1e-4), - demand_reactions=False, iterations=1) - return [r.id for r in solutions[0]] if solutions else [] - except Exception: # noqa: BLE001 — infeasible gap-fill / backend quirk: report none added + if (applied.slim_optimize(error_value=0.0) or 0.0) < floor - 1e-9: + return [] # unfixable even with every candidate present + applied.reactions.get_by_id(biomass_reaction).lower_bound = floor + fluxes = pfba(applied).fluxes + return [r.id for r in added if abs(fluxes.get(r.id, 0.0)) > 1e-9] + except Exception: # noqa: BLE001 — infeasible / backend quirk: report none added return [] +def _universal_met(model, met): + """The model's own metabolite matching ``met`` by id, or a fresh copy to be added with the reaction. + + Matches by id, as ``cobra.flux_analysis.gapfill`` did — the solve assumes the universal shares the + draft's metabolite namespace; a candidate whose metabolites do not resolve simply cannot carry flux + to biomass, so the LP leaves it out. + """ + if met.id in model.metabolites: + return model.metabolites.get_by_id(met.id) + return cobra.Metabolite(met.id, name=met.name, formula=met.formula, + charge=met.charge, compartment=met.compartment) + + def _diagnose_growth_gap(applied, proposal, biomass_reaction, movable_ids, pinned_comp, touches, base, compartments) -> tuple[str, str] | None: """If a biomass precursor is unproducible, pin the movable reaction that produces it to the diff --git a/tests/test_localization_certify.py b/tests/test_localization_certify.py index ec1105a..2100aff 100644 --- a/tests/test_localization_certify.py +++ b/tests/test_localization_certify.py @@ -316,6 +316,33 @@ def test_gapfill_restores_function(): assert _grows(apply_assignment(m, res, universal=_universal())) +def test_gapfill_returns_a_growth_restoring_set(): + # Directly exercise the flux-based _gapfill: on a draft that cannot grow it returns the universal + # reaction that restores biomass, and nothing spurious. As a plain LP it is reliable where cobra's + # indicator gap-fill MILP is not (that MILP fails to find such a fill in most genome-scale cases even + # when it exists) -- and a returned set always actually restores growth. + from raven_toolbox.localization.certify import _gapfill + + applied = _gap_draft() # EX_A -> r1: A->B; bio: B + C -> (C has no producer) + added = _gapfill(applied, _universal(), "bio", min_growth=1.0) + assert added == ["rC"] # the exact missing producer of C, nothing extra + + ur = _universal().reactions.get_by_id("rC") + nr = cobra.Reaction(ur.id, lower_bound=ur.lower_bound, upper_bound=ur.upper_bound) + applied.add_reactions([nr]) + nr.add_metabolites({applied.metabolites.get_by_id(m.id): c for m, c in ur.metabolites.items()}) + assert _grows(applied) + + +def test_gapfill_offers_nothing_when_the_universal_cannot_restore_growth(): + # If even adding every candidate cannot reach the floor, _gapfill returns [] (no false fill). + from raven_toolbox.localization.certify import _gapfill + + applied = _gap_draft() + empty = cobra.Model("empty") # no candidate produces C + assert _gapfill(applied, empty, "bio", min_growth=1.0) == [] + + def test_no_gratuitous_gapfill(): # When the draft already grows, no universal candidate is pulled (it is only reached on a real # growth failure). From 96708fde252941d175d5eefc5ed5cad5ba69eda5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 19 Jul 2026 23:23:41 +0200 Subject: [PATCH 2/2] gap-fill: run on a copy, batch-add candidates, warn on namespace mismatch Harden the flux-based _gapfill (universal-DB fill in localization.certify): operate on a model copy so the caller's model is never mutated; add candidates in one batch (per-reaction adds are super-linear at scale); return the flux-carrying set sorted, so the result does not depend on which co-optimal vertex the solver picked; and warn when most universal candidates share no metabolite id with the model, so a namespace mismatch is distinguishable from a genuine empty fill instead of returning [] silently. Knockout-recovery 60/60 (exact reaction each time) and 12/12 on realistic incomplete drafts, vs cobra's 45% and 0/12. --- CHANGELOG.md | 24 +++--- src/raven_toolbox/localization/certify.py | 94 ++++++++++++----------- tests/test_localization_certify.py | 16 ++++ 3 files changed, 81 insertions(+), 53 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c98c499..3ce0d0a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,16 +6,20 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased -* **`assign_compartments` gap-fill is now a reliable flux-based fill.** The universal-DB gap-fill in - `localization.certify` no longer calls `cobra.flux_analysis.gapfill`, whose indicator MILP, at genome - scale, fails to find a valid fill in the **majority** of cases *even when the exact reaction that - restores growth is present in the universal* — the MILP returns an incumbent its own validation then - rejects, and it raises rather than offering a fill. The replacement adds the universal candidates, holds - biomass at the growth floor, runs pFBA, and keeps the added reactions that carry flux: a plain LP that - cannot have that failure mode. Measured against the old path — knockout-recovery recall **100 % (60/60, - each recovering the exact removed reaction) vs cobra's 45 %**; on realistic incomplete drafts it restores - growth **12/12 vs cobra's 0/12**, ~5× faster. The caller still re-certifies with a real FBA, so no false - certificate is possible. +* **`assign_compartments` gap-fill: a reliable flux-based fill instead of cobra's MILP.** The + universal-DB gap-fill in `localization.certify` no longer calls `cobra.flux_analysis.gapfill`, whose + indicator MILP, at genome scale, fails to find a valid fill in the **majority** of cases *even when the + exact reaction that restores growth is present in the universal* (it returns an incumbent its own + validation rejects, and raises rather than offering a fill). The replacement adds the candidates on a + working copy, holds biomass at the floor, runs pFBA, and returns the flux-carrying additions — **sorted**, + so the result does not depend on which co-optimal vertex the solver picked. A plain LP cannot have the + MILP's failure mode: a returned set is a real flux solution that reaches the floor. On single-reaction + knockout-recovery (remove an essential reaction, fill from a universal that contains it), recall is + **60/60 — the exact reaction each time — vs cobra's 45 %**; on realistic incomplete drafts (12% of + internal reactions dropped, many simultaneous gaps) it restores growth **12/12 vs cobra's 0/12**. The universal must share the draft's + metabolite namespace: candidates are matched by id, as cobra's gapfill required, and a mismatch now + **warns** instead of returning `[]` silently. The set is flux-parsimonious (pFBA) but not guaranteed + reaction-count-minimal; the caller re-certifies with a real FBA, so no false certificate is possible. * **`diff_models` compares grRules as logic, not text.** The GPR check now DNF-expands each rule (via the existing `manipulation.gpr_to_dnf`), sorts the genes within each isozyme clause and sorts the clauses, so operand order no longer registers as a difference: `a and b` == `b and a` and `a or b` == `b or a`. The diff --git a/src/raven_toolbox/localization/certify.py b/src/raven_toolbox/localization/certify.py index b7bb5da..6fc3527 100644 --- a/src/raven_toolbox/localization/certify.py +++ b/src/raven_toolbox/localization/certify.py @@ -29,6 +29,7 @@ from __future__ import annotations import math +import warnings from collections.abc import Callable, Iterable, Mapping from dataclasses import dataclass @@ -697,59 +698,66 @@ def assign_compartments( def _gapfill(applied, universal, biomass_reaction, min_growth) -> list[str]: - """Universal reactions whose addition restores the primary growth floor, via a flux-based fill. - - Adds every ``universal`` candidate not already present, holds biomass at the floor, and runs pFBA - (maximise biomass, then minimise total flux); the added reactions that carry flux are returned. - - This is a plain LP, which matters. ``cobra.flux_analysis.gapfill`` solves an indicator MILP that, at - genome scale, fails to find a valid fill in the *majority* of cases even when the exact reaction that - restores growth is present in the universal — its own validation then rejects the bad incumbent and - raises, so no fill is offered. The LP here cannot have that failure mode: a returned set is a real - flux solution that actually reaches the floor. The caller re-certifies with FBA regardless, so this - never yields a false certificate. + """The ``universal`` reactions whose addition restores the primary growth floor, via a flux-based fill. + + On a working copy: add every ``universal`` candidate not already present in one batch, check the floor + is even reachable, hold biomass at the floor, run pFBA (maximise biomass, then minimise total flux), + and return -- sorted -- the added reactions that carry flux. pFBA makes the set flux-parsimonious (not + guaranteed reaction-count-minimal); the caller re-certifies with a real FBA regardless, so no false + certificate is possible. + + An LP, not cobra's indicator MILP — ``cobra.flux_analysis.gapfill`` at genome scale fails to find a + valid fill in the *majority* of cases even when the exact restoring reaction is present in the + universal (its own validation then rejects the broken incumbent and raises); on single-reaction + knockout-recovery this restores every case where that MILP restores under half. + + **Namespace.** Candidates are matched to the model by metabolite id (as cobra's gapfill required): the + universal must share the draft's metabolite namespace. A candidate whose metabolites do not resolve + becomes a dead-end that cannot carry flux and is left out — and a warning fires when most candidates + fail to resolve, so a silent empty result is distinguishable from a namespace mismatch. """ if biomass_reaction not in applied.reactions: return [] floor = max(min_growth, 1e-4) try: from cobra.flux_analysis import pfba - with applied: - added = [] - for urxn in universal.reactions: - if urxn.id in applied.reactions: - continue - new = cobra.Reaction(urxn.id, name=urxn.name, - lower_bound=urxn.lower_bound, upper_bound=urxn.upper_bound) - applied.add_reactions([new]) - new.add_metabolites({_universal_met(applied, m): coeff - for m, coeff in urxn.metabolites.items()}) - added.append(new) - if not added: - return [] - applied.objective = biomass_reaction - if (applied.slim_optimize(error_value=0.0) or 0.0) < floor - 1e-9: - return [] # unfixable even with every candidate present - applied.reactions.get_by_id(biomass_reaction).lower_bound = floor - fluxes = pfba(applied).fluxes - return [r.id for r in added if abs(fluxes.get(r.id, 0.0)) > 1e-9] + # Work on a copy: no context-manager rollback (whose failure on some optlang backends would + # otherwise discard a valid result), and the caller rebuilds `applied` from the proposal anyway. + work = applied.copy() + candidates = [u for u in universal.reactions if u.id not in work.reactions] + if not candidates: + return [] + fresh = [cobra.Reaction(u.id, name=u.name, lower_bound=u.lower_bound, upper_bound=u.upper_bound) + for u in candidates] + work.add_reactions(fresh) # one batch — per-reaction adds are super-linear at scale + fully_unresolved = 0 + for nr, urxn in zip(fresh, candidates, strict=True): + stoich, n_unres = {}, 0 + for met, coeff in urxn.metabolites.items(): + if met.id in work.metabolites: + stoich[work.metabolites.get_by_id(met.id)] = coeff + else: + n_unres += 1 + stoich[cobra.Metabolite(met.id, name=met.name, formula=met.formula, + charge=met.charge, compartment=met.compartment)] = coeff + nr.add_metabolites(stoich) + fully_unresolved += n_unres == len(urxn.metabolites) and len(urxn.metabolites) > 0 + if fully_unresolved > 0.5 * len(candidates): + warnings.warn( + f"gap-fill: {fully_unresolved}/{len(candidates)} universal candidates share no metabolite " + "id with the model — a likely namespace mismatch; gap-fill will find little or nothing.", + stacklevel=2) + + work.objective = biomass_reaction + if (work.slim_optimize(error_value=0.0) or 0.0) < floor - 1e-9: + return [] # unfixable even with every candidate present + work.reactions.get_by_id(biomass_reaction).lower_bound = floor + fluxes = pfba(work).fluxes + return sorted(nr.id for nr in fresh if abs(fluxes.get(nr.id, 0.0)) > 1e-9) except Exception: # noqa: BLE001 — infeasible / backend quirk: report none added return [] -def _universal_met(model, met): - """The model's own metabolite matching ``met`` by id, or a fresh copy to be added with the reaction. - - Matches by id, as ``cobra.flux_analysis.gapfill`` did — the solve assumes the universal shares the - draft's metabolite namespace; a candidate whose metabolites do not resolve simply cannot carry flux - to biomass, so the LP leaves it out. - """ - if met.id in model.metabolites: - return model.metabolites.get_by_id(met.id) - return cobra.Metabolite(met.id, name=met.name, formula=met.formula, - charge=met.charge, compartment=met.compartment) - - def _diagnose_growth_gap(applied, proposal, biomass_reaction, movable_ids, pinned_comp, touches, base, compartments) -> tuple[str, str] | None: """If a biomass precursor is unproducible, pin the movable reaction that produces it to the diff --git a/tests/test_localization_certify.py b/tests/test_localization_certify.py index 2100aff..b904794 100644 --- a/tests/test_localization_certify.py +++ b/tests/test_localization_certify.py @@ -343,6 +343,22 @@ def test_gapfill_offers_nothing_when_the_universal_cannot_restore_growth(): assert _gapfill(applied, empty, "bio", min_growth=1.0) == [] +def test_gapfill_warns_on_namespace_mismatch_instead_of_failing_silently(): + # A universal whose metabolite ids don't match the model can't connect; _gapfill returns [] but must + # WARN, so "found nothing" is distinguishable from "wrong namespace". + from raven_toolbox.localization.certify import _gapfill + + applied = _gap_draft() + foreign = cobra.Model("foreign") # same chemistry (A->C) under foreign ids + a, c = cobra.Metabolite("A_x", compartment="c"), cobra.Metabolite("C_x", compartment="c") + foreign.add_metabolites([a, c]) + rc = cobra.Reaction("rC", lower_bound=0, upper_bound=1000) + rc.add_metabolites({a: -1, c: 1}) + foreign.add_reactions([rc]) + with pytest.warns(UserWarning, match="namespace"): + assert _gapfill(applied, foreign, "bio", min_growth=1.0) == [] + + def test_no_gratuitous_gapfill(): # When the draft already grows, no universal candidate is pulled (it is only reached on a real # growth failure).