From 329c57e8bc5c0c6650030a8f612459e731fac896 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 18:15:34 +0100 Subject: [PATCH 01/27] Rank get_enzyme_bottlenecks by the re-solved change of the objective The order no longer depends on the shadow price, which depends on the basis at a degenerate optimum. Enzymes at their bound are solved again with the bound raised by fold_change; multiple_active marks rows where another capacity constraint is active. --- CHANGELOG.md | 1 + docs/migrating_from_gecko_matlab.md | 4 +- src/geckopy/utilities/bottlenecks.py | 98 ++++++++++++------ tests/test_bottlenecks.py | 149 ++++++++++++++++++++++++++- tests/test_scipy_solver_support.py | 8 +- 5 files changed, 217 insertions(+), 43 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e50d57..2a01a3d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `get_enzyme_bottlenecks` ranks the enzymes by `objective_gain`, the change of the objective after the upper bound of one enzyme's usage reaction is raised by `fold_change` (default 2.0) and the model is solved again, instead of by the absolute shadow price of the `prot_` constraint. At a degenerate optimum several capacity constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis. Only the enzymes at their bound (`upper_bound > 0` and `|flux| / upper_bound >= 1 - 1e-6`) are solved again; every other enzyme has `objective_gain` 0. Ties are broken by the position in `model.enzymes`, so the order does not depend on the solver. The columns are `gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active` (`shadow_price` is a reference column and moved from position 2 to 3), so callers that read columns by position must use the column names. `multiple_active` is True when another capacity constraint (an enzyme at its bound, or `prot_pool_exchange` at its upper bound) is active at the same optimum. `df.attrs["n_resolved"]` and `df.attrs["active_constraints"]` hold the number of solves and of active constraints. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index 6c3b465..5ea6eba 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/docs/migrating_from_gecko_matlab.md @@ -324,7 +324,7 @@ Each row maps a MATLAB function to its geckopy equivalent. The | `startGECKOproject` | `geckopy init` *(CLI)* | **`[Py]`** Project skeleton generator runs from the `geckopy` CLI (`geckopy init my_project`). | | `updateGECKOdoc` | *(not applicable)* | Python docs are docstring-driven; no MATLAB Toolbox metadata to regenerate. | | — | `pfba_enzymes` *(new)* | **`[Py]`** pFBA variant that minimises total `usage_prot_*` flux instead of all fluxes. Raises `NotImplementedError` on gecko-light. | -| — | `get_enzyme_bottlenecks` *(new)* | **`[Py]`** Per-enzyme dual-based bottleneck ranker; returns a DataFrame. | +| — | `get_enzyme_bottlenecks` *(new)* | **`[Py]`** Per-enzyme bottleneck ranker; returns a DataFrame sorted by the re-solved change of the objective when the bound of one enzyme at its bound is raised by `fold_change`. The shadow price is a reference column. | ### 5.7 Model adapter (`src/geckomat/model_adapter`) @@ -456,7 +456,7 @@ for the MATLAB-side items the port surfaced. ### geckopy-only (no MATLAB counterpart) - `pfba_enzymes` — enzyme-pool-minimising pFBA. -- `get_enzyme_bottlenecks` — dual-based bottleneck ranking. +- `get_enzyme_bottlenecks` — bottleneck ranking by the re-solved change of the objective. - `relax_proteomics_greedy` — greedy alternative to `flexibilize_enz_concs`. - `merge_kcats` — generalised n-source kcat merge (replaces the diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index ffdfdef..0c31ead 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -1,43 +1,42 @@ -"""Find the enzymes that are limiting the current objective. - -When you solve an ecModel, some enzymes are saturated and some -have spare capacity. The saturated ones are the ones whose -removal (or whose increased availability) would shift the -objective — they're the *bottlenecks*. The LP's shadow price for -each ``prot_`` mass-balance constraint quantifies that -sensitivity directly. - -``get_enzyme_bottlenecks`` runs the solve, ranks every enzyme by -absolute shadow price, and returns the top N as a DataFrame. The -top row is the enzyme whose extra availability would help the -objective the most. - -Ported from the legacy geckopy package described in Carrasco et al. -(2023, https://doi.org/10.1128/spectrum.01705-23), file -geckopy/flux_analysis.py:275-281 (get_protein_bottlenecks). +"""Find the enzymes that limit the current objective. + +When an ecModel is solved, some enzymes are used up to their upper bound and some have spare capacity. An enzyme that is not used up cannot limit the objective: raising its bound does not change the optimum. An enzyme that is used up may limit the objective, but it does not have to, because another capacity constraint can be active at the same optimum. + +``get_enzyme_bottlenecks`` ranks the enzymes by the re-solved change of the objective when the upper bound of one enzyme is raised by a factor ``fold_change``. Only the enzymes at their bound are re-solved; every other enzyme has a change of 0 without a solve. The LP shadow price of the ``prot_`` constraint is reported for reference and is not used for the order: at a degenerate optimum, which enzyme-constrained models routinely have, several constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis and not on the model. + +Ported from the legacy geckopy package described in Carrasco et al. (2023, https://doi.org/10.1128/spectrum.01705-23), file geckopy/flux_analysis.py:275-281 (get_protein_bottlenecks). """ from __future__ import annotations +import math from typing import TYPE_CHECKING import pandas as pd +from ..ec_model.constants import POOL_EXCHANGE_ID from ..solver_support import constraint_shadow_prices, optimize_model if TYPE_CHECKING: from ..ec_model.ec_model import EcModel +# An enzyme is at its bound when |flux| / upper_bound is at least 1 - _AT_CAP_TOL. +_AT_CAP_TOL = 1e-6 +# Changes of the objective below this value are set to 0. +_GAIN_TOL = 1e-10 + + +def get_enzyme_bottlenecks( + model: "EcModel", *, top: int = 10, fold_change: float = 2.0 +) -> pd.DataFrame: + """Top-N enzymes by the change of the objective when their bound is raised. -def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: - """Top-N enzymes by absolute shadow price. + Solves the model. An enzyme is at its bound when its upper bound is positive and ``|flux| / upper_bound`` is at least ``1 - 1e-6``. For each such enzyme the upper bound of its usage reaction is raised to ``upper_bound * fold_change``, the model is solved again and ``objective_gain`` is the new objective minus the objective of the first solve (0 when the new solve is not optimal or the change is below 1e-10). The model is restored after each solve. Every other enzyme has ``objective_gain`` 0 without a solve, because raising a bound that does not bind cannot change the optimum. The change is signed in the direction of the objective as the solver returns it, so it is positive for a maximised objective that increases. - Solves the model, ranks enzymes by ``|shadow_price|``, returns a - DataFrame indexed by uniprot id with columns - ``[gene, shadow_price, flux, cap_usage, upper_bound]``. + The rows are sorted by ``objective_gain`` descending, compared at 9 decimals, and ties keep the order of ``model.enzymes``. ``shadow_price`` (the dual of the ``prot_`` constraint) is reported for reference and does not enter the order, because at a degenerate optimum it depends on the basis the solver ended in. - Ported from the legacy geckopy package (Carrasco et al., 2023, - https://doi.org/10.1128/spectrum.01705-23), - geckopy/flux_analysis.py:275-281. + ``multiple_active`` is True for an enzyme when a capacity constraint other than its own is active at the same optimum, which is where the dual says little about the enzyme. The capacity constraints counted are the enzymes at their bound and the protein pool exchange (``prot_pool_exchange``) when its flux is within 1e-6 relative of a finite, positive upper bound. + + Ported from the legacy geckopy package (Carrasco et al., 2023, https://doi.org/10.1128/spectrum.01705-23), geckopy/flux_analysis.py:275-281. Parameters ---------- @@ -45,14 +44,13 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: EcModel to solve. top Number of top-ranked enzymes to return. Default 10. + fold_change + Factor by which the upper bound of an enzyme at its bound is raised. Default 2.0. Returns ------- pandas.DataFrame - Indexed by uniprot id, columns - ``[gene, shadow_price, flux, cap_usage, upper_bound]``, - sorted by descending ``|shadow_price|`` and truncated to - ``top`` rows. + Indexed by uniprot id, columns ``[gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active]``, truncated to ``top`` rows. ``df.attrs["n_resolved"]`` is the number of enzymes that were solved again and ``df.attrs["active_constraints"]`` the number of active capacity constraints. Raises ------ @@ -64,6 +62,7 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: raise RuntimeError( f"Solver status {sol.status!r}; cannot rank bottlenecks" ) + base = sol.objective_value enzymes = list(model.enzymes) prices = constraint_shadow_prices(model, [enz.prot_metabolite_id for enz in enzymes]) @@ -72,13 +71,46 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: rows.append({ "uniprot": enz.id, "gene": enz.gene, + "objective_gain": 0.0, "shadow_price": float(price), "flux": enz.flux, "cap_usage": enz.cap_usage, "upper_bound": enz.upper_bound, }) + + at_cap = [ + i for i, row in enumerate(rows) + if row["upper_bound"] > 0 and row["cap_usage"] >= 1 - _AT_CAP_TOL + ] + active = len(at_cap) + if model.reactions.has_id(POOL_EXCHANGE_ID): + pool = model.reactions.get_by_id(POOL_EXCHANGE_ID) + pool_ub = float(pool.upper_bound) + if 0 < pool_ub < math.inf and abs(float(pool.flux)) >= pool_ub * (1 - _AT_CAP_TOL): + active += 1 + + for i in at_cap: + usage = enzymes[i].usage_reaction + with model: + usage.upper_bound = rows[i]["upper_bound"] * fold_change + new = optimize_model(model) + if new.status == "optimal": + gain = new.objective_value - base + rows[i]["objective_gain"] = gain if gain >= _GAIN_TOL else 0.0 + if at_cap: + optimize_model(model) + + at_cap_set = set(at_cap) + for i, row in enumerate(rows): + row["multiple_active"] = (active - (1 if i in at_cap_set else 0)) >= 1 + + rows.sort(key=lambda row: -round(row["objective_gain"], 9)) df = pd.DataFrame(rows).set_index("uniprot") - df = df.reindex( - df["shadow_price"].abs().sort_values(ascending=False).index - ) - return df.head(top) + df = df[[ + "gene", "objective_gain", "shadow_price", "flux", "cap_usage", + "upper_bound", "multiple_active", + ]] + df = df.head(top) + df.attrs["n_resolved"] = len(at_cap) + df.attrs["active_constraints"] = active + return df diff --git a/tests/test_bottlenecks.py b/tests/test_bottlenecks.py index 8283845..d8a2246 100644 --- a/tests/test_bottlenecks.py +++ b/tests/test_bottlenecks.py @@ -2,9 +2,13 @@ from pathlib import Path import cobra +import numpy as np import pytest +from scipy import sparse from geckopy import ModelAdapter, make_ec_model +from geckopy.ec_model import EcModel +from geckopy.ec_model.ec_data import EcData from geckopy.ec_model.pipeline import apply_kcat_constraints from geckopy.utilities import get_enzyme_bottlenecks @@ -43,20 +47,155 @@ def test_returns_at_most_top_when_fewer_enzymes(): assert len(df) == len(model.ec.enzymes) -def test_sorted_by_absolute_shadow_price(): +def test_sorted_by_objective_gain(): model = _ectestgem_solvable() - df = get_enzyme_bottlenecks(model, top=10) - abs_sp = df["shadow_price"].abs().to_list() - assert abs_sp == sorted(abs_sp, reverse=True) + df = get_enzyme_bottlenecks(model, top=999) + rounded = df["objective_gain"].round(9).to_list() + assert rounded == sorted(rounded, reverse=True) + assert (df["objective_gain"] >= 0).all() def test_columns(): model = _ectestgem_solvable() df = get_enzyme_bottlenecks(model, top=3) assert list(df.columns) == [ - "gene", "shadow_price", "flux", "cap_usage", "upper_bound", + "gene", "objective_gain", "shadow_price", "flux", "cap_usage", + "upper_bound", "multiple_active", ] assert df.index.name == "uniprot" + assert "n_resolved" in df.attrs + assert "active_constraints" in df.attrs + + +def _toy(solver="glpk", *, ubs=(4.0, 1.0, 5.0, 3.0), pool_ub=1e6): + """Four enzymes. E1 and E2 are parallel steps A -> B, E3 turns B into C, E4 turns A into D. + The objective is SK_C + SK_D, the substrate is not limiting and R4 is capped at 1000. + + With the default bounds E1, E2 and E3 are at their bound (4 + 1 = 5) and raising any one of + them alone leaves the optimum unchanged, while E4 is at its bound and raising it adds flux. + """ + model = EcModel("toy") + mets = {n: cobra.Metabolite(n, compartment="c") for n in + ["A_e", "A", "B", "C", "D", "pool", "prot_E1", "prot_E2", "prot_E3", "prot_E4"]} + model.add_metabolites(list(mets.values())) + + def rxn(rid, stoich, lb, ub): + r = cobra.Reaction(rid) + r.add_metabolites({mets[k]: v for k, v in stoich.items()}) + r.bounds = (lb, ub) + return r + + reactions = [ + rxn("EX_A", {"A_e": -1.0}, -1e6, 0.0), + rxn("TR_A", {"A_e": -1.0, "A": 1.0}, 0.0, 1e6), + rxn("R1", {"A": -1.0, "B": 1.0, "prot_E1": -1.0}, 0.0, 1000.0), + rxn("R2", {"A": -1.0, "B": 1.0, "prot_E2": -1.0}, 0.0, 1000.0), + rxn("R3", {"B": -1.0, "C": 1.0, "prot_E3": -1.0}, 0.0, 1000.0), + rxn("R4", {"A": -1.0, "D": 1.0, "prot_E4": -1.0}, 0.0, 1000.0), + rxn("SK_C", {"C": -1.0}, 0.0, 1000.0), + rxn("SK_D", {"D": -1.0}, 0.0, 1000.0), + rxn("prot_pool_exchange", {"pool": 1.0}, 0.0, pool_ub), + ] + for i, ub in enumerate(ubs, start=1): + reactions.append(rxn(f"usage_prot_E{i}", {"pool": -1.0, f"prot_E{i}": 1.0}, 0.0, ub)) + model.add_reactions(reactions) + model.objective = {model.reactions.SK_C: 1.0, model.reactions.SK_D: 1.0} + model.ec = EcData( + rxns=["R1", "R2", "R3", "R4"], + kcat=np.ones(4), + source=[""] * 4, + notes=[""] * 4, + eccodes=[""] * 4, + genes=["g1", "g2", "g3", "g4"], + enzymes=["E1", "E2", "E3", "E4"], + mw=np.full(4, 100.0), + sequence=[""] * 4, + concs=np.array(ubs, dtype=float), + rxn_enz_mat=sparse.identity(4, format="csr"), + ) + try: + model.solver = solver + except Exception as exc: + pytest.skip(f"solver {solver} unavailable: {exc}") + return model + + +def test_toy_ranks_by_resolved_gain_not_by_dual(): + model = _toy() + df = get_enzyme_bottlenecks(model, top=10) + assert list(df.index) == ["E4", "E1", "E2", "E3"] + assert df.loc["E4", "objective_gain"] == pytest.approx(3.0) + assert df.loc[["E1", "E2", "E3"], "objective_gain"].to_list() == [0.0, 0.0, 0.0] + assert df.attrs["n_resolved"] == 4 + assert df.attrs["active_constraints"] == 4 + assert df["multiple_active"].all() + assert ((df["cap_usage"] - 1.0).abs() < 1e-6).all() + # A zero gain with a nonzero dual: the dual is reported but does not enter the order. + assert (df.loc[["E1", "E2", "E3"], "shadow_price"].abs() > 1e-9).any() + # The model is restored. + assert model.reactions.usage_prot_E4.upper_bound == 3.0 + assert model.enzymes.get_by_id("E4").flux == pytest.approx(3.0) + + +def test_non_binding_enzymes_are_not_resolved(): + model = _toy(ubs=(4.0, 3.0, 20.0, 5000.0)) + df = get_enzyme_bottlenecks(model, top=10) + # E1 and E2 carry 7 through E3 (ub 20) and are at their bound; E3 is not, and E4 is limited by R4. + assert df.attrs["n_resolved"] == 2 + assert df.attrs["active_constraints"] == 2 + assert df.loc[["E3", "E4"], "objective_gain"].to_list() == [0.0, 0.0] + assert df.loc["E1", "objective_gain"] > 0 + assert df.loc["E1", "multiple_active"] + assert list(df.index[:2]) == ["E1", "E2"] + + +def test_single_active_constraint_is_not_multiple_active(): + model = _toy(ubs=(10.0, 10.0, 2.0, 5000.0)) + df = get_enzyme_bottlenecks(model, top=10) + assert df.attrs["n_resolved"] == 1 + assert df.attrs["active_constraints"] == 1 + assert df.loc["E3", "objective_gain"] == pytest.approx(2.0) + assert not df.loc["E3", "multiple_active"] + assert df.loc[["E1", "E2", "E4"], "multiple_active"].all() + assert list(df.index)[0] == "E3" + + +def test_pool_exchange_at_its_bound_is_an_active_constraint(): + model = _toy(ubs=(10.0, 10.0, 2.0, 5000.0), pool_ub=2.0) + df = get_enzyme_bottlenecks(model, top=10) + assert df.attrs["active_constraints"] == 1 + assert df.attrs["n_resolved"] == 0 + assert df["multiple_active"].all() + assert (df["objective_gain"] == 0).all() + + +def test_pool_exchange_and_an_enzyme_at_their_bounds(): + model = _toy(ubs=(10.0, 10.0, 2.0, 5000.0), pool_ub=1004.0) + df = get_enzyme_bottlenecks(model, top=10) + assert df.attrs["active_constraints"] == 2 + assert df.attrs["n_resolved"] == 1 + assert df.loc["E3", "multiple_active"] + assert df.loc["E3", "objective_gain"] == 0.0 + + +def test_fold_change_is_respected(): + twice = get_enzyme_bottlenecks(_toy(), top=1, fold_change=2.0) + thrice = get_enzyme_bottlenecks(_toy(), top=1, fold_change=3.0) + assert twice.index[0] == thrice.index[0] == "E4" + assert twice["objective_gain"].iloc[0] == pytest.approx(3.0) + assert thrice["objective_gain"].iloc[0] == pytest.approx(6.0) + + +@pytest.mark.parametrize("solver", ["glpk", "scipy", "gurobi"]) +def test_order_does_not_depend_on_the_solver(solver): + reference = get_enzyme_bottlenecks(_toy("glpk"), top=10) + df = get_enzyme_bottlenecks(_toy(solver), top=10) + assert list(df.index) == list(reference.index) + assert df["objective_gain"].to_numpy() == pytest.approx( + reference["objective_gain"].to_numpy(), abs=1e-6 + ) + assert df["multiple_active"].to_list() == reference["multiple_active"].to_list() + assert df.attrs["n_resolved"] == reference.attrs["n_resolved"] def test_raises_on_infeasible(): diff --git a/tests/test_scipy_solver_support.py b/tests/test_scipy_solver_support.py index b0a155f..8e4c102 100644 --- a/tests/test_scipy_solver_support.py +++ b/tests/test_scipy_solver_support.py @@ -161,13 +161,15 @@ def test_get_enzyme_bottlenecks_on_scipy_equals_the_reference(reference, scipy_d scipy_model = _on("scipy", tight=True)[0] actual = get_enzyme_bottlenecks(scipy_model, top=999) assert (expected["shadow_price"].abs() > 1e-9).any() - assert list(actual.index[:1]) == list(expected.index[:1]) - assert set(actual.index) == set(expected.index) + assert list(actual.index) == list(expected.index) + assert actual["objective_gain"].to_numpy() == pytest.approx( + expected["objective_gain"].to_numpy(), abs=1e-6 + ) aligned = actual.loc[expected.index] assert aligned["shadow_price"].to_numpy() == pytest.approx( expected["shadow_price"].to_numpy(), abs=1e-6 ) - assert actual["shadow_price"].abs().is_monotonic_decreasing + assert actual["objective_gain"].round(9).is_monotonic_decreasing for uniprot, row in actual.iterrows(): assert row["flux"] == pytest.approx(scipy_model.enzymes.get_by_id(uniprot).flux) From dfa2172a2b2bfd24ab110fd6a531329239772bbe Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 18:17:11 +0100 Subject: [PATCH 02/27] get_enzyme_bottlenecks reports the gain in the objective's direction and requires a fold_change above 1 --- src/geckopy/utilities/bottlenecks.py | 10 +++++++--- tests/test_bottlenecks.py | 16 ++++++++++++++++ 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index 0c31ead..6695438 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -32,7 +32,7 @@ def get_enzyme_bottlenecks( Solves the model. An enzyme is at its bound when its upper bound is positive and ``|flux| / upper_bound`` is at least ``1 - 1e-6``. For each such enzyme the upper bound of its usage reaction is raised to ``upper_bound * fold_change``, the model is solved again and ``objective_gain`` is the new objective minus the objective of the first solve (0 when the new solve is not optimal or the change is below 1e-10). The model is restored after each solve. Every other enzyme has ``objective_gain`` 0 without a solve, because raising a bound that does not bind cannot change the optimum. The change is signed in the direction of the objective as the solver returns it, so it is positive for a maximised objective that increases. - The rows are sorted by ``objective_gain`` descending, compared at 9 decimals, and ties keep the order of ``model.enzymes``. ``shadow_price`` (the dual of the ``prot_`` constraint) is reported for reference and does not enter the order, because at a degenerate optimum it depends on the basis the solver ended in. + ``objective_gain`` is the improvement of the objective in its own direction: a rise for a maximised objective and a fall for a minimised one. The rows are sorted by ``objective_gain`` descending, compared at 9 decimals, and ties keep the order of ``model.enzymes``. ``shadow_price`` (the dual of the ``prot_`` constraint) is reported for reference and does not enter the order, because at a degenerate optimum it depends on the basis the solver ended in. ``multiple_active`` is True for an enzyme when a capacity constraint other than its own is active at the same optimum, which is where the dual says little about the enzyme. The capacity constraints counted are the enzymes at their bound and the protein pool exchange (``prot_pool_exchange``) when its flux is within 1e-6 relative of a finite, positive upper bound. @@ -45,7 +45,7 @@ def get_enzyme_bottlenecks( top Number of top-ranked enzymes to return. Default 10. fold_change - Factor by which the upper bound of an enzyme at its bound is raised. Default 2.0. + Factor, greater than 1, by which the upper bound of an enzyme at its bound is raised. Default 2.0. Returns ------- @@ -54,9 +54,13 @@ def get_enzyme_bottlenecks( Raises ------ + ValueError + If ``fold_change`` is not greater than 1. RuntimeError If the solve status is not optimal. """ + if not fold_change > 1: + raise ValueError(f"fold_change must be greater than 1, got {fold_change!r}") sol = optimize_model(model) if sol.status != "optimal": raise RuntimeError( @@ -95,7 +99,7 @@ def get_enzyme_bottlenecks( usage.upper_bound = rows[i]["upper_bound"] * fold_change new = optimize_model(model) if new.status == "optimal": - gain = new.objective_value - base + gain = (new.objective_value - base) * (1.0 if model.objective_direction == "max" else -1.0) rows[i]["objective_gain"] = gain if gain >= _GAIN_TOL else 0.0 if at_cap: optimize_model(model) diff --git a/tests/test_bottlenecks.py b/tests/test_bottlenecks.py index d8a2246..59761c9 100644 --- a/tests/test_bottlenecks.py +++ b/tests/test_bottlenecks.py @@ -186,6 +186,22 @@ def test_fold_change_is_respected(): assert thrice["objective_gain"].iloc[0] == pytest.approx(6.0) +@pytest.mark.parametrize("fold_change", [1.0, 0.5, float("nan")]) +def test_fold_change_must_exceed_one(fold_change): + with pytest.raises(ValueError, match="fold_change"): + get_enzyme_bottlenecks(_toy(), fold_change=fold_change) + + +def test_gain_is_the_improvement_of_a_minimised_objective(): + maximised = get_enzyme_bottlenecks(_toy(), top=4) + model = _toy() + model.objective = {model.reactions.get_by_id("SK_C"): -1.0, model.reactions.get_by_id("SK_D"): -1.0} + model.objective_direction = "min" + minimised = get_enzyme_bottlenecks(model, top=4) + assert list(minimised.index) == list(maximised.index) + assert minimised["objective_gain"].to_numpy() == pytest.approx(maximised["objective_gain"].to_numpy()) + + @pytest.mark.parametrize("solver", ["glpk", "scipy", "gurobi"]) def test_order_does_not_depend_on_the_solver(solver): reference = get_enzyme_bottlenecks(_toy("glpk"), top=10) From a2d533795b20da0b40a5d3f6583ca10477d9049c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:13:12 +0100 Subject: [PATCH 03/27] get_enzyme_bottlenecks logs a warning when no enzyme is at its usage bound --- CHANGELOG.md | 2 +- src/geckopy/utilities/bottlenecks.py | 15 +++++++++++++++ tests/test_bottlenecks.py | 25 +++++++++++++++++++++++++ 3 files changed, 41 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a01a3d..6abf96d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,7 @@ ## Unreleased -* `get_enzyme_bottlenecks` ranks the enzymes by `objective_gain`, the change of the objective after the upper bound of one enzyme's usage reaction is raised by `fold_change` (default 2.0) and the model is solved again, instead of by the absolute shadow price of the `prot_` constraint. At a degenerate optimum several capacity constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis. Only the enzymes at their bound (`upper_bound > 0` and `|flux| / upper_bound >= 1 - 1e-6`) are solved again; every other enzyme has `objective_gain` 0. Ties are broken by the position in `model.enzymes`, so the order does not depend on the solver. The columns are `gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active` (`shadow_price` is a reference column and moved from position 2 to 3), so callers that read columns by position must use the column names. `multiple_active` is True when another capacity constraint (an enzyme at its bound, or `prot_pool_exchange` at its upper bound) is active at the same optimum. `df.attrs["n_resolved"]` and `df.attrs["active_constraints"]` hold the number of solves and of active constraints. +* `get_enzyme_bottlenecks` ranks the enzymes by `objective_gain`, the change of the objective after the upper bound of one enzyme's usage reaction is raised by `fold_change` (default 2.0) and the model is solved again, instead of by the absolute shadow price of the `prot_` constraint. At a degenerate optimum several capacity constraints are active and the solver assigns the dual to one of them, so the dual depends on the basis. Only the enzymes at their bound (`upper_bound > 0` and `|flux| / upper_bound >= 1 - 1e-6`) are solved again; every other enzyme has `objective_gain` 0. Ties are broken by the position in `model.enzymes`, so the order does not depend on the solver. The columns are `gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active` (`shadow_price` is a reference column and moved from position 2 to 3), so callers that read columns by position must use the column names. `multiple_active` is True when another capacity constraint (an enzyme at its bound, or `prot_pool_exchange` at its upper bound) is active at the same optimum. `df.attrs["n_resolved"]` and `df.attrs["active_constraints"]` hold the number of solves and of active constraints. When no enzyme is at its usage bound, as on ecYeastGEM where the protein pool is the active constraint, a warning says so and the rows are in model order. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index 6695438..d7d96da 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -8,6 +8,7 @@ """ from __future__ import annotations +import logging import math from typing import TYPE_CHECKING @@ -19,6 +20,8 @@ if TYPE_CHECKING: from ..ec_model.ec_model import EcModel +logger = logging.getLogger(__name__) + # An enzyme is at its bound when |flux| / upper_bound is at least 1 - _AT_CAP_TOL. _AT_CAP_TOL = 1e-6 # Changes of the objective below this value are set to 0. @@ -52,6 +55,10 @@ def get_enzyme_bottlenecks( pandas.DataFrame Indexed by uniprot id, columns ``[gene, objective_gain, shadow_price, flux, cap_usage, upper_bound, multiple_active]``, truncated to ``top`` rows. ``df.attrs["n_resolved"]`` is the number of enzymes that were solved again and ``df.attrs["active_constraints"]`` the number of active capacity constraints. + Warns + ----- + When no enzyme is at its usage bound a warning is logged: no single enzyme limits the objective, every ``objective_gain`` is 0 and the rows are in model order. + Raises ------ ValueError @@ -87,11 +94,19 @@ def get_enzyme_bottlenecks( if row["upper_bound"] > 0 and row["cap_usage"] >= 1 - _AT_CAP_TOL ] active = len(at_cap) + pool_active = False if model.reactions.has_id(POOL_EXCHANGE_ID): pool = model.reactions.get_by_id(POOL_EXCHANGE_ID) pool_ub = float(pool.upper_bound) if 0 < pool_ub < math.inf and abs(float(pool.flux)) >= pool_ub * (1 - _AT_CAP_TOL): active += 1 + pool_active = True + if not at_cap: + logger.warning( + "get_enzyme_bottlenecks: no enzyme is at its usage bound, so no single enzyme limits " + "the objective; %s. Every objective_gain is 0 and the rows are in model order.", + "the protein pool exchange is at its bound" if pool_active else "no capacity constraint is active", + ) for i in at_cap: usage = enzymes[i].usage_reaction diff --git a/tests/test_bottlenecks.py b/tests/test_bottlenecks.py index 59761c9..6a4e989 100644 --- a/tests/test_bottlenecks.py +++ b/tests/test_bottlenecks.py @@ -221,3 +221,28 @@ def test_raises_on_infeasible(): r3.lower_bound = 1e9 with pytest.raises(RuntimeError, match="Solver status"): get_enzyme_bottlenecks(model, top=3) + + +def test_warns_when_no_enzyme_is_at_its_bound(caplog): + model = _toy(ubs=(1e5, 1e5, 1e5, 1e5)) + with caplog.at_level("WARNING", logger="geckopy.utilities.bottlenecks"): + df = get_enzyme_bottlenecks(model) + assert df.attrs["n_resolved"] == 0 + assert (df["objective_gain"] == 0).all() + assert "no enzyme is at its usage bound" in caplog.text + assert "no capacity constraint is active" in caplog.text + + +def test_warning_names_the_protein_pool_when_it_is_the_active_constraint(caplog): + model = _toy(ubs=(1e5, 1e5, 1e5, 1e5), pool_ub=6.0) + with caplog.at_level("WARNING", logger="geckopy.utilities.bottlenecks"): + df = get_enzyme_bottlenecks(model) + assert df.attrs["n_resolved"] == 0 + assert df.attrs["active_constraints"] == 1 + assert "the protein pool exchange is at its bound" in caplog.text + + +def test_no_warning_when_an_enzyme_is_at_its_bound(caplog): + with caplog.at_level("WARNING", logger="geckopy.utilities.bottlenecks"): + get_enzyme_bottlenecks(_toy()) + assert "no enzyme is at its usage bound" not in caplog.text From d7bfb7f709c76c5b067fabb167ae2d92d6ae5106 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:13:34 +0100 Subject: [PATCH 04/27] Migration notes: no dash in the get_enzyme_bottlenecks bullet --- docs/migrating_from_gecko_matlab.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/migrating_from_gecko_matlab.md b/docs/migrating_from_gecko_matlab.md index 5ea6eba..ae8128d 100644 --- a/docs/migrating_from_gecko_matlab.md +++ b/docs/migrating_from_gecko_matlab.md @@ -456,7 +456,7 @@ for the MATLAB-side items the port surfaced. ### geckopy-only (no MATLAB counterpart) - `pfba_enzymes` — enzyme-pool-minimising pFBA. -- `get_enzyme_bottlenecks` — bottleneck ranking by the re-solved change of the objective. +- `get_enzyme_bottlenecks`: bottleneck ranking by the re-solved change of the objective. - `relax_proteomics_greedy` — greedy alternative to `flexibilize_enz_concs`. - `merge_kcats` — generalised n-source kcat merge (replaces the From 89fbddad88c735ebdb7a96411b0b0a4e980acd2e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:32:16 +0100 Subject: [PATCH 05/27] Evotune: skip NaN rmse, stable screen order, full generations on a flat objective best_parsimonious ignores points with a NaN rmse, screen_kcat_leverage sorts stably, cmaes_kcat_tuning does not stop on a flat generation, and compute_excarbon warns for a model without metabolite formulas. --- CHANGELOG.md | 3 ++ src/geckopy/kcat_tuning/evotune/distance.py | 11 +++++- src/geckopy/kcat_tuning/evotune/parsimony.py | 9 +++-- src/geckopy/kcat_tuning/evotune/tuning.py | 21 +++++++++-- tests/test_evotune.py | 38 ++++++++++++++++++++ tests/test_evotune_parsimony.py | 14 ++++++++ 6 files changed, 90 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e50d57..f816a95 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,9 @@ ## Unreleased +* `best_parsimonious` skips points with a NaN `rmse`, as MATLAB's `bestParsimonious` does with `min(rmses)`: a NaN at the first point of the frontier made every comparison false and raised `ValueError`. It still raises when every `rmse` is NaN. +* `screen_kcat_leverage` keeps rows of equal rank in the order of the tie groups (a stable sort), as MATLAB's `screenKcatLeverage` does, and `compute_excarbon` warns when no metabolite of the model has a formula (every exchange then gets the carbon weight 1). +* `cmaes_kcat_tuning` runs `max_generations` generations on a flat objective instead of stopping after the first (the flat-fitness criteria of `cma` do not end the search), as MATLAB's `cmaesKcatTuning` does. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/kcat_tuning/evotune/distance.py b/src/geckopy/kcat_tuning/evotune/distance.py index f6cc757..ee84195 100644 --- a/src/geckopy/kcat_tuning/evotune/distance.py +++ b/src/geckopy/kcat_tuning/evotune/distance.py @@ -16,6 +16,7 @@ """ from __future__ import annotations +import warnings from typing import TYPE_CHECKING, Iterable import numpy as np @@ -49,7 +50,9 @@ def compute_excarbon( hand-rolled formula parser). A missing/unparseable formula, or a computed carbon count of exactly 0 (e.g. O2, water), both fall back to 1 -- matching MATLAB's ``Ematrix(isnan(Ematrix)) = 1`` and - ``excarbon(excarbon == 0) = 1`` respectively. + ``excarbon(excarbon == 0) = 1`` respectively. A model in which no + metabolite has a formula gives every exchange the weight 1 and + emits a warning. Parameters ---------- @@ -67,6 +70,12 @@ def compute_excarbon( dict[str, float] """ excarbon: dict[str, float] = {} + if not any(m.formula for m in model.metabolites): + warnings.warn( + "The model has no metabolite formulas; every exchange " + "reaction gets the carbon weight 1.", + UserWarning, stacklevel=2, + ) for rxn_id in rxn_ids: if rxn_id == bio_rxn_id: excarbon[rxn_id] = BIOMASS_CARBON_EQUIV diff --git a/src/geckopy/kcat_tuning/evotune/parsimony.py b/src/geckopy/kcat_tuning/evotune/parsimony.py index 7f3ffd3..3da3b74 100644 --- a/src/geckopy/kcat_tuning/evotune/parsimony.py +++ b/src/geckopy/kcat_tuning/evotune/parsimony.py @@ -239,10 +239,15 @@ def best_parsimonious( Defaults to 2%, which is the scale of the alternate-optimum noise measured on this model -- below that, two vectors are not - meaningfully different in fit, so the sparser one wins. + meaningfully different in fit, so the sparser one wins. A point with + a NaN RMSE is never chosen and does not set the best RMSE; when every + RMSE is NaN no point qualifies and ``ValueError`` is raised. """ if not points: raise ValueError("points is empty.") - best_rmse = min(p.rmse for p in points) + rmses = [p.rmse for p in points if not np.isnan(p.rmse)] + if not rmses: + raise ValueError("every point has a NaN rmse.") + best_rmse = min(rmses) within = [p for p in points if p.rmse <= best_rmse * (1.0 + tolerance)] return min(within, key=lambda p: p.n_changed) diff --git a/src/geckopy/kcat_tuning/evotune/tuning.py b/src/geckopy/kcat_tuning/evotune/tuning.py index 6c448f2..9585e68 100644 --- a/src/geckopy/kcat_tuning/evotune/tuning.py +++ b/src/geckopy/kcat_tuning/evotune/tuning.py @@ -185,6 +185,10 @@ def _resolve_context(model, adapter, params, evotune_data, bio_rxn, okp_method): return params, evotune_data, bio_rxn, okp_method +#: Termination criteria of ``cma`` that a flat objective triggers. +_FLAT_STOPS = frozenset({"tolfun", "tolfunhist", "tolflatfitness"}) + + def _tunable_context( model, params, evotune_data, bio_rxn, okp_method, tunable_mask, ): @@ -351,7 +355,8 @@ def screen_kcat_leverage( df = pd.DataFrame(rows) df["rank"] = df["leverage"] * df["sigma0_log"] - df = df.sort_values("rank", ascending=False, ignore_index=True) + df = df.sort_values( + "rank", ascending=False, ignore_index=True, kind="stable") total = float(df["rank"].sum()) df["cum_leverage_share"] = ( (df["rank"].cumsum() / total) if total > 0 else 0.0 @@ -412,7 +417,9 @@ def cmaes_kcat_tuning( Configuration comes entirely from ``EvotuneParams``: trust tiers, ``max_growth_weight``, ``prior_penalty_weight`` and ``tie_isozymes`` control the search, and ``max_generations``/``rmse_threshold`` are - its stopping conditions. + its stopping conditions. A generation whose candidates all score + the same does not end the search early; ``cma``'s other convergence + criteria (such as ``tolx``) still can. Parameters ---------- @@ -556,12 +563,20 @@ def _score_batch(kcat_matrix, pool): rmse_trace: list[float] = [] objective_trace: list[float] = [] generation = 0 - while not es.stop(): + flat = False + while True: + stop = es.stop() + # cma ends a search whose candidates all score the same + # (tolfun, tolfunhist, tolflatfitness); that is not a stopping + # condition here, so a flat generation does not end the run. + if stop and not (flat and set(stop) <= _FLAT_STOPS): + break asks = es.ask() xm = np.array(asks).T scores = _score_batch(expand(xm), pool) es.tell(asks, list(scores[:, 0])) generation += 1 + flat = bool(np.ptp(scores[:, 0]) == 0.0) j = int(np.argmin(scores[:, 0])) if scores[j, 0] < best_obj: diff --git a/tests/test_evotune.py b/tests/test_evotune.py index b926671..fc73710 100644 --- a/tests/test_evotune.py +++ b/tests/test_evotune.py @@ -443,3 +443,41 @@ def test_result_reports_its_own_corrections(): assert rows[0].source == "brenda" assert rows[0].fold_change == pytest.approx(10.0) assert rows[0].cumulative_share == pytest.approx(3.0 / 4.01) + + +def test_screen_orders_rows_of_equal_rank_as_the_tie_groups(tmp_path, monkeypatch): + from geckopy.kcat_tuning.evotune import tuning + + monkeypatch.setattr(tuning, "_score_kcat_vector", lambda *a, **k: (1.0, 1.0)) + adapter = _adapter(tmp_path) + model = _build_toy(adapter) + screen = screen_kcat_leverage( + model, adapter=adapter, evotune_data=_evotune_data(), n_proc=1, + ) + assert list(screen["leverage"]) == [0.0, 0.0] + assert list(screen["rxn_id"]) == list(model.ec.rxns) + + +def test_cmaes_tuning_runs_all_generations_on_a_flat_objective(tmp_path, monkeypatch): + from geckopy.kcat_tuning.evotune import tuning + + monkeypatch.setattr(tuning, "_score_kcat_vector", lambda *a, **k: (1.0, 1.0)) + adapter = _adapter(tmp_path) + model = _build_toy(adapter) + params = EvotuneParams(max_generations=6, rmse_threshold=-1.0) + result = cmaes_kcat_tuning( + model, adapter=adapter, params=params, evotune_data=_evotune_data(), + n_proc=1, seed=0, verbose=False, + ) + assert result.n_generations == 6 + + +def test_screen_warns_when_the_model_has_no_formulas(tmp_path): + adapter = _adapter(tmp_path) + model = _build_toy(adapter) + for met in model.metabolites: + met.formula = None + with pytest.warns(UserWarning, match="no metabolite formulas"): + screen_kcat_leverage( + model, adapter=adapter, evotune_data=_evotune_data(), n_proc=1, + ) diff --git a/tests/test_evotune_parsimony.py b/tests/test_evotune_parsimony.py index ddb3e5a..3d360da 100644 --- a/tests/test_evotune_parsimony.py +++ b/tests/test_evotune_parsimony.py @@ -72,6 +72,20 @@ def test_best_parsimonious_prefers_fewer_changes_within_tolerance(): assert best_parsimonious(pts, tolerance=0.0).n_changed == 4000 +def test_best_parsimonious_skips_nan_rmse(): + from geckopy.kcat_tuning.evotune.parsimony import FrontierPoint + + pts = [ + FrontierPoint(0.0, np.array([1.0]), 4000, 7.0, float("nan")), + FrontierPoint(3.0, np.array([1.0]), 3000, 6.0, 0.905), + FrontierPoint(6.0, np.array([1.0]), 1000, 3.0, 1.400), + ] + assert best_parsimonious(pts).n_changed == 3000 + nan_only = [FrontierPoint(0.0, np.array([1.0]), 1, 1.0, float("nan"))] + with pytest.raises(ValueError): + best_parsimonious(nan_only) + + def test_fold_change_is_symmetric_and_at_least_one(): from geckopy.kcat_tuning.evotune.parsimony import fold_change From 5fe37b5f06265a5c068180a4e1bef9987939ce19 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:34:38 +0100 Subject: [PATCH 06/27] Export review_assignment, EcStats and ec_stats_from_brenda --- CHANGELOG.md | 1 + src/geckopy/__init__.py | 6 ++++ src/geckopy/gather_kcats/__init__.py | 4 +++ tests/test_review_public_api.py | 42 ++++++++++++++++++++++++++++ 4 files changed, 53 insertions(+) create mode 100644 tests/test_review_public_api.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e50d57..72fecf4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `review_assignment`, `EcStats` and `ec_stats_from_brenda` are exported from `geckopy.gather_kcats` and `geckopy`, and listed in `__all__`. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/__init__.py b/src/geckopy/__init__.py index c5b5711..a23a0ff 100644 --- a/src/geckopy/__init__.py +++ b/src/geckopy/__init__.py @@ -73,10 +73,12 @@ # Gather-kcats from geckopy.gather_kcats import ( + EcStats, apply_kcat_list, assign_kcats, assign_standard_kcat, collect_kcats, + ec_stats_from_brenda, fetch_open_kinetics_predictor, fuzzy_kcat_matching, get_standard_kcat, # deprecated; alias of assign_standard_kcat @@ -91,6 +93,7 @@ read_dlkcat_output, read_kcat_source, remove_standard_kcat, + review_assignment, run_dlkcat, save_kcat_source, select_kcat_value, # deprecated; alias of apply_kcat_list @@ -175,6 +178,7 @@ # Core data + builder + per-enzyme accessor "EcData", "EcModel", + "EcStats", "Enzyme", "make_ec_model", # Databases (loaders + dataclasses) @@ -209,6 +213,7 @@ "assign_kcats", "assign_standard_kcat", "collect_kcats", + "ec_stats_from_brenda", "fetch_open_kinetics_predictor", "fuzzy_kcat_matching", "kcats_from_brenda", @@ -222,6 +227,7 @@ "read_dlkcat_output", "read_kcat_source", "remove_standard_kcat", + "review_assignment", "run_dlkcat", "save_kcat_source", "submit_open_kinetics_predictor", diff --git a/src/geckopy/gather_kcats/__init__.py b/src/geckopy/gather_kcats/__init__.py index 795816b..f5688ab 100644 --- a/src/geckopy/gather_kcats/__init__.py +++ b/src/geckopy/gather_kcats/__init__.py @@ -18,6 +18,7 @@ ) from .read_dlkcat_output import read_dlkcat_output from .remove_standard_kcat import remove_standard_kcat +from .review import EcStats, ec_stats_from_brenda, review_assignment from .run_dlkcat import run_dlkcat from .select_kcat_value import ( apply_kcat_list, @@ -27,6 +28,7 @@ from .write_dlkcat_input import extract_enzyme_substrate_pairs, write_dlkcat_input __all__ = [ + "EcStats", "OKPClient", "OKPError", "apply_kcat_list", @@ -34,6 +36,7 @@ "assign_standard_kcat", "collect_kcats", "ec_rxn_keys", + "ec_stats_from_brenda", "extract_enzyme_substrate_pairs", "fetch_open_kinetics_predictor", "format_kcat_source", @@ -51,6 +54,7 @@ "read_kcat_source", "read_dlkcat_output", "remove_standard_kcat", + "review_assignment", "run_dlkcat", "save_kcat_source", "select_kcat_value", # deprecated; alias of apply_kcat_list diff --git a/tests/test_review_public_api.py b/tests/test_review_public_api.py new file mode 100644 index 0000000..fe704a4 --- /dev/null +++ b/tests/test_review_public_api.py @@ -0,0 +1,42 @@ +"""The kcat review functions are public API of geckopy and geckopy.gather_kcats.""" +import geckopy +from geckopy import gather_kcats +from geckopy.databases import load_brenda_data + +_NAMES = ("EcStats", "ec_stats_from_brenda", "review_assignment") + + +def test_review_functions_are_exported_from_both_packages(): + for name in _NAMES: + assert name in geckopy.__all__ + assert name in gather_kcats.__all__ + assert getattr(geckopy, name) is getattr(gather_kcats, name) + + +def test_ec_stats_keep_the_case_of_the_ec_number(tmp_path): + (tmp_path / "kcat.tsv").write_text( + "# BRENDA release test\n" + "ec_code\tsubstrate\torganism\tkcat_max\tkcat_median\tn\treferences\n" + "3.5.1.n3\tm1\torg\t40\t4\t2\t*\n" + "1.1.1.1\tm1\torg\t10\t2\t2\t*\n" + "1.1.1.1\tm2\torg\t6\t3\t2\t*\n", + encoding="utf-8", + ) + (tmp_path / "sa.tsv").write_text( + "# BRENDA release test\n" + "ec_code\tsubstrate\torganism\tsa_max\tsa_median\tn\treferences\n", + encoding="utf-8", + ) + (tmp_path / "mw.tsv").write_text( + "# BRENDA release test\nec_code\tsubstrate\torganism\tmw\tn\treferences\n", + encoding="utf-8", + ) + stats = geckopy.ec_stats_from_brenda(load_brenda_data(tmp_path, filter_outlier_kcats=False)) + assert sorted(stats) == ["1.1.1.1", "3.5.1.n3"] + assert stats["1.1.1.1"].kcat_max == 10.0 + assert stats["1.1.1.1"].kcat_median == 2.5 + assert stats["3.5.1.n3"].values == frozenset({4.0, 40.0}) + findings = geckopy.review_assignment( + ["R1"], [40.0], ["brenda"], ["3.5.1.n3"], ec_stats=stats, + ) + assert [f.checks for f in findings] == [("ec-maximum",)] From 2f1cdc1409ca561cc546e10c37d5e3193cdab582 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:44 +0100 Subject: [PATCH 07/27] fill_enz_concs and constrain_enz_concs raise for a negative concentration --- .../limit_proteins/constrain_enz_concs.py | 25 +++++++++++-------- src/geckopy/limit_proteins/fill_enz_concs.py | 18 ++++++------- tests/test_constrain_enz_concs.py | 18 +++++++++++++ tests/test_fill_enz_concs.py | 19 ++++++++++++++ 4 files changed, 61 insertions(+), 19 deletions(-) diff --git a/src/geckopy/limit_proteins/constrain_enz_concs.py b/src/geckopy/limit_proteins/constrain_enz_concs.py index f7ce43c..eb0591a 100644 --- a/src/geckopy/limit_proteins/constrain_enz_concs.py +++ b/src/geckopy/limit_proteins/constrain_enz_concs.py @@ -62,7 +62,9 @@ def constrain_enz_concs( Raises ------ ValueError - If the ``prot_pool`` metabolite is missing, or if any + If ``model.ec.concs`` holds a negative concentration for a + processed enzyme (unless ``remove_constraints`` is set), if the + ``prot_pool`` metabolite is missing, or if any ``usage_prot_`` reaction is missing for an enzyme in ``model.ec.enzymes`` (e.g. because the model is gecko-light, which has no usage reactions). @@ -94,7 +96,18 @@ def constrain_enz_concs( if enz in restrict_set ) - for i, enz in index_iter: + index_list = list(index_iter) + if not remove_constraints: + negative = [ + enz for i, enz in index_list if model.ec.concs[i] < 0 + ] + if negative: + raise ValueError( + "model.ec.concs has a negative concentration for " + + ", ".join(negative) + "." + ) + + for i, enz in index_list: rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{enz}") rxn.upper_bound = 1000.0 if remove_constraints: @@ -102,12 +115,4 @@ def constrain_enz_concs( conc = model.ec.concs[i] if np.isnan(conc): continue - if conc < 0: - # A negative bound would make the usage reaction infeasible - # (ub < lb == 0). Skip it and leave the enzyme unconstrained. - logger.warning( - "constrain_enz_concs: negative concentration %g for %r; " - "leaving it unconstrained.", float(conc), enz, - ) - continue rxn.upper_bound = float(conc) diff --git a/src/geckopy/limit_proteins/fill_enz_concs.py b/src/geckopy/limit_proteins/fill_enz_concs.py index 4829b39..ad86e17 100644 --- a/src/geckopy/limit_proteins/fill_enz_concs.py +++ b/src/geckopy/limit_proteins/fill_enz_concs.py @@ -59,12 +59,13 @@ def fill_enz_concs( Raises ------ + ValueError + If a protein that is an enzyme of the model has a negative + abundance. IndexError If ``data_col`` is out of range for the supplied abundances. """ n_enz = model.ec.n_enzymes - model.ec.concs = np.full(n_enz, np.nan, dtype=float) - if prot_data.abundances.ndim == 2: if not 0 <= data_col < prot_data.abundances.shape[1]: raise IndexError( @@ -80,7 +81,9 @@ def fill_enz_concs( ) column = prot_data.abundances + concs = np.full(n_enz, np.nan, dtype=float) if n_enz == 0 or len(prot_data.uniprot_ids) == 0: + model.ec.concs = concs return enzyme_to_idx = {enz: i for i, enz in enumerate(model.ec.enzymes)} @@ -90,11 +93,8 @@ def fill_enz_concs( continue value = float(conc) if value < 0: - # A measured concentration can't be negative; skip it so it - # doesn't become an infeasible bound in constrain_enz_concs. - logger.warning( - "fill_enz_concs: negative concentration %g for %r; skipped.", - value, uid, + raise ValueError( + f"prot_data has a negative abundance ({value:g}) for {uid!r}." ) - continue - model.ec.concs[idx] = value + concs[idx] = value + model.ec.concs = concs diff --git a/tests/test_constrain_enz_concs.py b/tests/test_constrain_enz_concs.py index 9089435..45908ba 100644 --- a/tests/test_constrain_enz_concs.py +++ b/tests/test_constrain_enz_concs.py @@ -215,3 +215,21 @@ def test_integration_with_fill_enz_concs(): assert model.reactions.get_by_id("usage_prot_P1").upper_bound == 4.0 assert model.reactions.get_by_id("usage_prot_P2").upper_bound == 1000.0 assert model.reactions.get_by_id("usage_prot_P3").upper_bound == 9.0 + + +# --------------------------------------------------------------------------- # +# Negative concentration +# --------------------------------------------------------------------------- # + +def test_negative_conc_raises_before_changing_the_model(): + model = _ec_model_with_pool(["P1", "P2"], concs=[3.0, -0.5]) + with pytest.raises(ValueError, match="negative concentration"): + constrain_enz_concs(model) + assert model.reactions.get_by_id("usage_prot_P1").upper_bound == 1000.0 + assert model.reactions.get_by_id("usage_prot_P2").upper_bound == 1000.0 + + +def test_negative_conc_is_not_read_when_removing_constraints(): + model = _ec_model_with_pool(["P1"], concs=[-0.5]) + constrain_enz_concs(model, remove_constraints=True) + assert model.reactions.get_by_id("usage_prot_P1").upper_bound == 1000.0 diff --git a/tests/test_fill_enz_concs.py b/tests/test_fill_enz_concs.py index 02c5b7b..069c6be 100644 --- a/tests/test_fill_enz_concs.py +++ b/tests/test_fill_enz_concs.py @@ -191,3 +191,22 @@ def test_integration_with_load_pax_db(tmp_path): fill_enz_concs(model, prot) # P2 -> index 0 (60), P1 -> index 1 (20). np.testing.assert_array_equal(model.ec.concs, [60.0, 20.0]) + + +# --------------------------------------------------------------------------- # +# Negative abundance +# --------------------------------------------------------------------------- # + +def test_negative_abundance_of_an_enzyme_raises(): + model = _ec_model(["P1", "P2"], initial_concs=[1.0, 2.0]) + pd = _prot_data([("P1", 5.0), ("P2", -0.5)]) + with pytest.raises(ValueError, match="negative abundance"): + fill_enz_concs(model, pd) + np.testing.assert_array_equal(model.ec.concs, [1.0, 2.0]) + + +def test_negative_abundance_of_a_protein_outside_the_model_is_not_read(): + model = _ec_model(["P1"]) + pd = _prot_data([("P1", 5.0), ("X1", -1.0)]) + fill_enz_concs(model, pd) + np.testing.assert_array_equal(model.ec.concs, [5.0]) From bd7a608b341fa11403f9cfd9ac8618772e47f3a2 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:44 +0100 Subject: [PATCH 08/27] load_flux_data raises for a ragged row and reads nested parentheses in headers --- src/geckopy/databases/flux_data.py | 67 +++++++++++++++++++----------- tests/test_load_flux_data.py | 26 ++++++++++++ 2 files changed, 68 insertions(+), 25 deletions(-) diff --git a/src/geckopy/databases/flux_data.py b/src/geckopy/databases/flux_data.py index 3665544..4971512 100644 --- a/src/geckopy/databases/flux_data.py +++ b/src/geckopy/databases/flux_data.py @@ -4,7 +4,6 @@ from __future__ import annotations import logging -import re from dataclasses import dataclass, field from pathlib import Path from typing import Optional @@ -15,7 +14,6 @@ # ` ()` -> capture met_name and rxn_id. -_EXCH_HEADER_RE = re.compile(r"^(.*)\s*\(([^()]+)\)\s*$") @dataclass @@ -81,6 +79,11 @@ def load_flux_data(path: str | Path) -> FluxData: * ``evotuneRMSEweight`` (float): weights for evotune kcat tuning. * ``source`` (str): free-text describing the data origin. + An exchange column header has the form `` ()``, split at + the parenthesis that matches the final closing one, so a reaction id + may contain parentheses (``X (EX_a(e))`` gives ``X`` and + ``EX_a(e)``). + Both columns are dropped from the matrix and surfaced on the returned ``FluxData`` as ``evotune_rmse_weight`` and ``source``. @@ -101,7 +104,8 @@ def load_flux_data(path: str | Path) -> FluxData: FileNotFoundError If ``path`` does not exist. ValueError - If the file has no header row or no exchange columns. + If the file has no header row, no exchange columns, or a row + with a different number of cells than the header. """ path = Path(path) if not path.is_file(): @@ -123,6 +127,12 @@ def load_flux_data(path: str | Path) -> FluxData: header = rows[0] data = rows[1:] + for line_no, r in enumerate(data, start=2): + if len(r) != len(header): + raise ValueError( + f"{path}: row {line_no} has {len(r)} cell(s), the header " + f"has {len(header)}." + ) # Pull out optional columns first; collect indices to drop. optional_indices: list[int] = [] @@ -134,15 +144,12 @@ def load_flux_data(path: str | Path) -> FluxData: if tuning_idx is not None: tuning_weights = np.array( - [_parse_float(r[tuning_idx]) if tuning_idx < len(r) else float("nan") - for r in data], + [_parse_float(r[tuning_idx]) for r in data], dtype=float, ) optional_indices.append(tuning_idx) if source_idx is not None: - source_values = [ - r[source_idx] if source_idx < len(r) else "" for r in data - ] + source_values = [r[source_idx] for r in data] optional_indices.append(source_idx) # Build the kept-column index list (preserving order). @@ -158,22 +165,13 @@ def load_flux_data(path: str | Path) -> FluxData: exch_headers = kept_header[3:] exch_mets, exch_rxn_ids = _split_exch_headers(exch_headers) - conds = [r[kept[0]] if kept[0] < len(r) else "" for r in data] - p_tot = np.array( - [_parse_float(r[kept[1]]) if kept[1] < len(r) else float("nan") - for r in data], - dtype=float, - ) - gr_rate = np.array( - [_parse_float(r[kept[2]]) if kept[2] < len(r) else float("nan") - for r in data], - dtype=float, - ) + conds = [r[kept[0]] for r in data] + p_tot = np.array([_parse_float(r[kept[1]]) for r in data], dtype=float) + gr_rate = np.array([_parse_float(r[kept[2]]) for r in data], dtype=float) exch_fluxes = np.full((len(data), len(exch_headers)), np.nan, dtype=float) for i, r in enumerate(data): for ci, j in enumerate(kept[3:]): - if j < len(r): - exch_fluxes[i, ci] = _parse_float(r[j]) + exch_fluxes[i, ci] = _parse_float(r[j]) return FluxData( conds=conds, @@ -206,6 +204,25 @@ def _parse_float(s: str) -> float: return float("nan") +def _split_header(header: str) -> Optional[tuple[str, str]]: + """Split ``' ()'`` at the parenthesis that matches the final + closing one; return None for a header of another form.""" + if not header.endswith(")"): + return None + depth = 0 + for k in range(len(header) - 1, -1, -1): + if header[k] == ")": + depth += 1 + elif header[k] == "(": + depth -= 1 + if depth == 0: + rxn = header[k + 1:-1].strip() + if not rxn: + return None + return header[:k].strip(), rxn + return None + + def _split_exch_headers(headers: list[str]) -> tuple[list[str], list[str]]: """Split each header `' ()'` into (met, rxn) pairs. @@ -217,10 +234,10 @@ def _split_exch_headers(headers: list[str]) -> tuple[list[str], list[str]]: rxns: list[str] = [] no_paren: list[str] = [] for h in headers: - m = _EXCH_HEADER_RE.match(h.strip()) - if m: - mets.append(m.group(1).strip()) - rxns.append(m.group(2).strip()) + parts = _split_header(h.strip()) + if parts is not None: + mets.append(parts[0]) + rxns.append(parts[1]) else: mets.append(h.strip()) rxns.append(h.strip()) diff --git a/tests/test_load_flux_data.py b/tests/test_load_flux_data.py index baa5188..7c32541 100644 --- a/tests/test_load_flux_data.py +++ b/tests/test_load_flux_data.py @@ -289,3 +289,29 @@ def test_loaded_data_can_feed_apply_flux_data_constraints(tmp_path): apply_flux_data_constraints(model, fd) assert model.reactions.get_by_id("EX_glc").lower_bound == -5.0 assert model.reactions.get_by_id("biomass").upper_bound == pytest.approx(0.4) + + +# --------------------------------------------------------------------------- # +# Malformed rows and headers +# --------------------------------------------------------------------------- # + +def test_ragged_row_raises(tmp_path): + p = tmp_path / "ragged.tsv" + _write_tsv(p, [ + ["condition", "Ptot", "grRate", "Glc (EX_a)", "O2 (EX_b)"], + ["c1", "0.5", "0.4", "-1", "2"], + ["c2", "0.5", "0.4", "-1"], + ]) + with pytest.raises(ValueError, match="row 3 has 4 cell"): + load_flux_data(p) + + +def test_nested_parentheses_in_reaction_id(tmp_path): + p = tmp_path / "nested.tsv" + _write_tsv(p, [ + ["condition", "Ptot", "grRate", "X (EX_a(e))", "O2 (EX_b)"], + ["c1", "0.5", "0.4", "-1", "2"], + ]) + fd = load_flux_data(p) + assert fd.exch_mets == ["X", "O2"] + assert fd.exch_rxn_ids == ["EX_a(e)", "EX_b"] From cb47f39f66ffc443e6fb5b06292836f14e37a938 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:44 +0100 Subject: [PATCH 09/27] enzyme_usage warns about enzymes without usage reaction; report_enzyme_usage rejects a negative top_abs_usage --- src/geckopy/utilities/enzyme_usage.py | 14 ++++++++++++++ src/geckopy/utilities/report_enzyme_usage.py | 9 +++++++++ tests/test_enzyme_usage.py | 13 +++++++++++++ tests/test_report_enzyme_usage.py | 11 +++++++++++ 4 files changed, 47 insertions(+) diff --git a/src/geckopy/utilities/enzyme_usage.py b/src/geckopy/utilities/enzyme_usage.py index 328d462..ac780dc 100644 --- a/src/geckopy/utilities/enzyme_usage.py +++ b/src/geckopy/utilities/enzyme_usage.py @@ -5,6 +5,7 @@ """ from __future__ import annotations +import logging from dataclasses import dataclass, field from typing import TYPE_CHECKING, Mapping, Union @@ -17,6 +18,8 @@ from ..ec_model.constants import USAGE_PREFIX +logger = logging.getLogger(__name__) + @dataclass class EnzymeUsageResult: @@ -83,6 +86,9 @@ def enzyme_usage( When False, enzymes with zero absolute usage are dropped from the result. + An enzyme without a ``usage_prot_`` reaction in the model is + left out of the result, with a logged warning. + Raises ------ NotImplementedError @@ -98,9 +104,11 @@ def enzyme_usage( abs_usage_list: list[float] = [] ub_list: list[float] = [] + left_out: list[str] = [] for enzyme in model.ec.enzymes: rxn_id = f"{USAGE_PREFIX}{enzyme}" if rxn_id not in cobra_rxn_ids: + left_out.append(enzyme) continue flux = float(_lookup_flux(fluxes, rxn_id)) rxn_ub = float(model.reactions.get_by_id(rxn_id).upper_bound) @@ -108,6 +116,12 @@ def enzyme_usage( abs_usage_list.append(abs(flux)) ub_list.append(rxn_ub) + if left_out: + logger.warning( + "enzyme_usage: the model has no usage reaction for the enzymes " + "%s; they are left out.", ", ".join(left_out), + ) + abs_usage_arr = np.array(abs_usage_list, dtype=float) ub_arr = np.array(ub_list, dtype=float) with np.errstate(divide="ignore", invalid="ignore"): diff --git a/src/geckopy/utilities/report_enzyme_usage.py b/src/geckopy/utilities/report_enzyme_usage.py index 9526e5c..65af600 100644 --- a/src/geckopy/utilities/report_enzyme_usage.py +++ b/src/geckopy/utilities/report_enzyme_usage.py @@ -108,7 +108,16 @@ def report_enzyme_usage( Returns ------- EnzymeUsageReport + + Raises + ------ + ValueError + If ``top_abs_usage`` is negative or NaN. """ + if not top_abs_usage >= 0: + raise ValueError( + f"top_abs_usage must be a non-negative number, got {top_abs_usage!r}." + ) pool_rxn = model.reactions.get_by_id(POOL_EXCHANGE_ID) total_pool = float(pool_rxn.upper_bound) diff --git a/tests/test_enzyme_usage.py b/tests/test_enzyme_usage.py index b6b2888..4b3869d 100644 --- a/tests/test_enzyme_usage.py +++ b/tests/test_enzyme_usage.py @@ -189,3 +189,16 @@ def test_fluxes_carried_through_to_result(): result = enzyme_usage(model, sol.fluxes) # Same Series object (or equivalent values). assert result.fluxes is sol.fluxes + + +# --------------------------------------------------------------------------- # +# Enzyme without a usage reaction +# --------------------------------------------------------------------------- # + +def test_enzyme_without_usage_reaction_is_left_out_with_a_warning(caplog): + model = _build_toy(enzyme_ub=5.0) + model.ec.enzymes.append("F") + with caplog.at_level("WARNING"): + result = enzyme_usage(model, {"usage_prot_E": 3.0}) + assert result.prot_id == ["E"] + assert "no usage reaction for the enzymes F" in caplog.text diff --git a/tests/test_report_enzyme_usage.py b/tests/test_report_enzyme_usage.py index 2c41082..3a6eda9 100644 --- a/tests/test_report_enzyme_usage.py +++ b/tests/test_report_enzyme_usage.py @@ -314,3 +314,14 @@ def test_top_abs_usage_inf_returns_all(): ud = enzyme_usage(model, sol.fluxes) rep = report_enzyme_usage(model, ud, top_abs_usage=float("inf")) assert len(rep.top_abs_usage) == 1 + + +# --------------------------------------------------------------------------- # +# Invalid top_abs_usage +# --------------------------------------------------------------------------- # + +def test_negative_top_abs_usage_raises(): + model = _build_single_enzyme_model() + ud = enzyme_usage(model, model.optimize().fluxes) + with pytest.raises(ValueError, match="top_abs_usage"): + report_enzyme_usage(model, ud, top_abs_usage=-2) From 7897ec71188f5b8cad9aea676619617e97e45d3e Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:45 +0100 Subject: [PATCH 10/27] Add tests for the shared contract of calculate_f_factor, relax_proteomics_greedy, apply_complex_data and add_new_rxns_to_ec --- tests/test_add_new_rxns_to_ec.py | 22 ++++++++++++++++++++++ tests/test_apply_complex_data.py | 17 +++++++++++++++++ tests/test_calculate_f_factor.py | 11 +++++++++++ tests/test_relax_proteomics_greedy.py | 9 +++++++++ 4 files changed, 59 insertions(+) diff --git a/tests/test_add_new_rxns_to_ec.py b/tests/test_add_new_rxns_to_ec.py index 7818750..8ade448 100644 --- a/tests/test_add_new_rxns_to_ec.py +++ b/tests/test_add_new_rxns_to_ec.py @@ -318,3 +318,25 @@ def test_multiple_new_rxns_added(tmp_path): assert sorted(result.enz_added) == ["Y", "Z"] assert "R_AY" in model.ec.rxns assert "R_BZ" in model.ec.rxns + + +# --------------------------------------------------------------------------- # +# Genes that already have enzymes +# --------------------------------------------------------------------------- # + +def test_reaction_with_genes_that_already_have_enzymes_is_added(tmp_path): + adapter = _adapter(tmp_path) + model = _build_base_model(adapter) + A_c = model.metabolites.get_by_id("A_c") + B_c = cobra.Metabolite("B_c", compartment="c") + model.add_metabolites([B_c]) + n_enzymes = len(model.ec.enzymes) + rxn = _make_rxn( + "R_AX", metabolites={A_c: -1.0, B_c: 1.0}, gpr="g_X", + ) + + result = add_new_rxns_to_ec(model, [rxn], []) + + assert result.enz_added == [] + assert len(model.ec.enzymes) == n_enzymes + assert "R_AX" in model.ec.rxns diff --git a/tests/test_apply_complex_data.py b/tests/test_apply_complex_data.py index 29f5cb7..8f26a34 100644 --- a/tests/test_apply_complex_data.py +++ b/tests/test_apply_complex_data.py @@ -283,3 +283,20 @@ def test_missing_path_and_no_complex_data_raises(tmp_path): with pytest.raises(FileNotFoundError): apply_complex_data(ec_model) + + +def test_partial_zero_stoichiometry_counts_the_subunit_as_one(): + """A subunit with stoichiometry 0 has no stated stoichiometry: it is a + subunit with stoichiometry 1, so the complex still matches.""" + ec_model = _ectestgem_ec_model() + data = [ComplexPortalEntry( + complex_id="R2Compl", + name="x", + species="t", + gene_names=["G1", "G2"], + protein_ids=["P1", "P2"], + stoichiometry=[2, 0], + )] + apply_complex_data(ec_model, complex_data=data, apply=False) + assert _stoich(ec_model, "R2_EXP_1", "P1") == 2.0 + assert _stoich(ec_model, "R2_EXP_1", "P2") == 1.0 diff --git a/tests/test_calculate_f_factor.py b/tests/test_calculate_f_factor.py index cae1d40..cc23daf 100644 --- a/tests/test_calculate_f_factor.py +++ b/tests/test_calculate_f_factor.py @@ -178,3 +178,14 @@ def test_integration_with_pax_db_loader(tmp_path): f = calculate_f_factor(model, pd_data) # P1+P2 = 80, total = 110. assert f == pytest.approx(80.0 / 110.0) + + +# --------------------------------------------------------------------------- # +# Explicit empty enzyme list +# --------------------------------------------------------------------------- # + +def test_explicit_empty_enzyme_list_means_no_enzymes(): + model = _ec_model(["P1", "P2"]) + pd = _prot_data([("P1", 100.0), ("P2", 100.0)]) + assert calculate_f_factor(model, pd, enzymes=[]) == 0.0 + assert calculate_f_factor(model, pd) == pytest.approx(1.0) diff --git a/tests/test_relax_proteomics_greedy.py b/tests/test_relax_proteomics_greedy.py index 80a62aa..a7fdc1a 100644 --- a/tests/test_relax_proteomics_greedy.py +++ b/tests/test_relax_proteomics_greedy.py @@ -99,3 +99,12 @@ def test_max_iterations_raises(): relax_proteomics_greedy( model, minimal_growth=1e9, max_iterations=1, ) + + +def test_empty_enzyme_set_makes_no_enzyme_eligible(): + model, _ = _ec_with_tight_concs() + result = relax_proteomics_greedy( + model, minimal_growth=0.001, enzyme_set=[], max_iterations=20, + ) + assert result.relaxed == {} + assert result.trace == [] From b7ce0ee618ff6d9dc2f8dab8d035e0b6adeefa2d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:37:45 +0100 Subject: [PATCH 11/27] Add the CHANGELOG entries of round 5 group 3 --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e50d57..1ebae89 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,11 @@ ## Unreleased +* `fill_enz_concs` raises a `ValueError` for a negative abundance of a protein that is an enzyme of the model, as GECKO's `fillEnzConcs` does with the error `fillEnzConcs:negativeAbundance`; it used to skip the value with a warning. +* `constrain_enz_concs` raises a `ValueError` for a negative concentration in `model.ec.concs` before it changes any bound, as GECKO's `constrainEnzConcs` does with the error `constrainEnzConcs:negativeConcentration`; it used to skip the enzyme with a warning. `remove_constraints=True` does not read the concentrations. +* `load_flux_data` raises a `ValueError` for a row with a different number of cells than the header, as GECKO's `loadFluxData` does with the error `loadFluxData:raggedRow`; it used to pad the missing cells with NaN. A header `X (EX_a(e))` gives the metabolite `X` and the reaction id `EX_a(e)`: the reaction id is the text between the parenthesis that matches the final one and that one; a header with nested parentheses used to give the whole header for both. +* `enzyme_usage` logs a warning that lists the enzymes without a `usage_prot_` reaction, which it leaves out of the result, as GECKO's `getEnzymeUsage` does with the warning `getEnzymeUsage:missingUsageReaction`. +* `report_enzyme_usage` raises a `ValueError` for a negative `top_abs_usage`; it used to list all enzymes except the last ones, counted from the end of the ranking. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. From 40ce95fd6253532d202ac7ede410ad6ad1410ee7 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:28 +0100 Subject: [PATCH 12/27] read_kcat_source skips a row of empty quoted fields and raises for a field with an opening quote and no closing one --- CHANGELOG.md | 1 + src/geckopy/gather_kcats/kcat_sources.py | 19 +++++-- tests/test_kcat_contract.py | 71 ++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 5 deletions(-) create mode 100644 tests/test_kcat_contract.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e50d57..ed3b7e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `read_kcat_source` and `load_kcat_source` skip a row whose fields are all empty quoted fields (`""""`), as they skip a row of tabs, and raise `ValueError` for a field that opens a double quote and does not close it. Before, the first was reported as an invalid row and the second was read as a row with a kcat that is not a number. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/gather_kcats/kcat_sources.py b/src/geckopy/gather_kcats/kcat_sources.py index 41ef74c..f3437b1 100644 --- a/src/geckopy/gather_kcats/kcat_sources.py +++ b/src/geckopy/gather_kcats/kcat_sources.py @@ -344,8 +344,10 @@ def read_kcat_source(path: str | Path) -> tuple[pd.DataFrame, dict[str, str]]: The file is read as UTF-8 (a leading byte-order mark is ignored). The ``#`` lines before the column header are the header; ``#`` lines and - blank lines further down are skipped. Fields wrapped in double quotes, - as a spreadsheet may save them, are unwrapped. ``kcat`` and + blank lines further down are skipped, as are rows without a filled field + (a row of empty quoted fields included). Fields wrapped in double quotes, + as a spreadsheet may save them, are unwrapped; a field that opens a + double quote and does not close it raises ``ValueError``. ``kcat`` and ``log10_sd`` must be plain numbers (``12``, ``0.5``, ``1e-05``); other text, such as ``5,5`` with a decimal comma, is read as NaN, which :func:`load_kcat_source` reports. A row with more filled fields than @@ -389,9 +391,16 @@ def read_kcat_source(path: str | Path) -> tuple[pd.DataFrame, dict[str, str]]: raise ValueError(f"{path}, line {j + 1} has more fields than the {len(columns)} columns.") fields = fields[: len(columns)] fields += [""] * (len(columns) - len(fields)) - if not any(f.strip() for f in fields): - continue # a row of tabs only - rows.append([_unquote(f) for f in fields]) + for f in fields: + t = f.strip() + if t.startswith('"') and not (len(t) >= 2 and t.endswith('"')): + raise ValueError( + f"{path}, line {j + 1} has a field with an opening double quote and no closing one." + ) + values = [_unquote(f) for f in fields] + if not any(values): + continue # a row without a filled field + rows.append(values) line_numbers.append(j + 1) keep = [i for i, c in enumerate(columns) if c] records = pd.DataFrame([[r[i] for i in keep] for r in rows], columns=[columns[i] for i in keep], dtype=object) diff --git a/tests/test_kcat_contract.py b/tests/test_kcat_contract.py new file mode 100644 index 0000000..ef145aa --- /dev/null +++ b/tests/test_kcat_contract.py @@ -0,0 +1,71 @@ +"""Tests for the behaviour GECKO's kcatContractTests check as well.""" +import shutil +from pathlib import Path + +import cobra +import numpy as np +import pandas as pd +import pytest + +from geckopy import ModelAdapter, make_ec_model +from geckopy.databases import UniprotDB, load_dlkcat_ignore_lists +from geckopy.gather_kcats import ( + assign_standard_kcat, + kcats_from_ec_tree, + read_dlkcat_output, + read_kcat_source, + write_dlkcat_input, +) +from geckopy.gather_kcats.ec_tree import ECTree +from geckopy.get_enzyme_data.ec_from_gem import fill_eccodes_from_gem + +EXAMPLE_DIR = Path(__file__).parents[1] / "examples" / "ecTestGEM" + + +def _project(tmp_path: Path, extra_reactions: bool = False, light: bool = False): + project = tmp_path / "ecTestGEM" + shutil.copytree(EXAMPLE_DIR, project) + brenda = project / "data" / "brenda" + brenda.mkdir() + for name in ("kcat.tsv", "sa.tsv", "mw.tsv"): + shutil.copy(project / "data" / name, brenda / name) + adapter = ModelAdapter.from_folder(project) + conv = cobra.io.read_sbml_model(str(adapter.params.conv_gem)) + if extra_reactions: + m1 = conv.metabolites.get_by_id("m1c") + m2 = conv.metabolites.get_by_id("m2c") + for rxn_id, name, stoich in (("R6", "R6", {m1: -1.0, m2: 1.0}), + ("R7", "R7", {m1: -1.0, m2: -1.0}), + ("R8", "Biomass Pseudoreaction", {m1: -1.0, m2: 1.0})): + rxn = cobra.Reaction(rxn_id, name=name) + rxn.lower_bound = 0.0 + rxn.upper_bound = 1000.0 + rxn.add_metabolites(stoich) + conv.add_reactions([rxn]) + model = make_ec_model(conv, adapter, gecko_light=light) + fill_eccodes_from_gem(model) + return model, adapter, project + + +def _uniprot() -> UniprotDB: + mws = [100.0, 200.0, 300.0] + return UniprotDB(ids=["P000", "P001", "P002"], genes=["g0", "g1", "g2"], eccodes=["1.1.1.1"] * 3, + mw=np.array(mws), sequences=[""] * 3) + + +def test_read_kcat_source_quoted_rows(tmp_path): + """A row of empty quoted fields is blank and a quoted comment is a comment; a + field that opens a quote and does not close it is an error.""" + file = tmp_path / "odd.tsv" + head = ["# source: odd", "rxn\tenzymes\tkcat\tsource", "R3\tP4\t25\todd"] + file.write_text("\n".join(head + ['"# a quoted comment"', '""\t""\t""\t""', "R5\tP5\t30\todd"]) + "\n", + encoding="utf-8") + records, _ = read_kcat_source(file) + assert list(records["kcat"]) == [25.0, 30.0] + assert list(records["line"]) == [3, 6] + file.write_text("\n".join(head + ['"# a comment that was quoted']) + "\n", encoding="utf-8") + with pytest.raises(ValueError, match="line 4"): + read_kcat_source(file) + file.write_text("\n".join(head[:2] + ['"R3\tP4\t25\todd']) + "\n", encoding="utf-8") + with pytest.raises(ValueError, match="no closing one"): + read_kcat_source(file) From a066de990f97b75bbabb9b6fa8a7cc6edad31db2 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:28 +0100 Subject: [PATCH 13/27] read_dlkcat_output drops infinite kcats; write_dlkcat_input writes the rows in the order of model.ec.rxns --- CHANGELOG.md | 2 + .../gather_kcats/read_dlkcat_output.py | 11 ++++-- .../gather_kcats/write_dlkcat_input.py | 8 ++-- tests/test_kcat_contract.py | 37 +++++++++++++++++++ 4 files changed, 51 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ed3b7e7..8156adb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,8 @@ ## Unreleased * `read_kcat_source` and `load_kcat_source` skip a row whose fields are all empty quoted fields (`""""`), as they skip a row of tabs, and raise `ValueError` for a field that opens a double quote and does not close it. Before, the first was reported as an invalid row and the second was read as a row with a kcat that is not a number. +* `read_dlkcat_output` drops a row whose kcat is infinite (`inf`), as it drops `NA` and `None`, so the kcat list holds finite numbers only. +* `write_dlkcat_input` and `extract_enzyme_substrate_pairs` write the rows in the order of `model.ec.rxns` whatever the order of `ec_rxns`, and a reaction named twice once, as GECKO's `writeDLKcatInput` does with its logical vector. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/gather_kcats/read_dlkcat_output.py b/src/geckopy/gather_kcats/read_dlkcat_output.py index f91e46d..2795fc8 100644 --- a/src/geckopy/gather_kcats/read_dlkcat_output.py +++ b/src/geckopy/gather_kcats/read_dlkcat_output.py @@ -9,6 +9,7 @@ from pathlib import Path from typing import TYPE_CHECKING +import numpy as np import pandas as pd if TYPE_CHECKING: @@ -43,8 +44,9 @@ def read_dlkcat_output( The input file is the 6-column TSV that DLKcat writes (see ``run_dlkcat``), with one header line and columns ``rxn_id, gene, substrate, smiles, sequence, kcat``. Rows whose - kcat column is non-numeric (``"NA"``, empty, etc.) are silently - dropped before validation against the model. + kcat column is not a finite number (``"NA"``, ``"None"``, empty, + ``"inf"``, ``"1.2.3"``) are dropped before the reactions are validated + against the model. The output DataFrame matches the schema produced by ``fuzzy_kcat_matching`` so downstream functions @@ -87,8 +89,8 @@ def read_dlkcat_output( ValueError If the file contains no numeric kcat values; if ``strict=True`` and any substrate name from the file is not - in ``model.metabolites.`` (case-insensitive); or if any - ``rxn_id`` is not in ``model.ec.rxns``. + in ``model.metabolites.`` (case-insensitive); or if the + ``rxn_id`` of a row that is kept is not in ``model.ec.rxns``. """ file_path = Path(file_path) if not file_path.is_file(): @@ -104,6 +106,7 @@ def read_dlkcat_output( ) numeric_kcat = pd.to_numeric(df["kcat"], errors="coerce") + numeric_kcat = numeric_kcat.where(np.isfinite(numeric_kcat)) valid_mask = numeric_kcat.notna() if not valid_mask.any(): diff --git a/src/geckopy/gather_kcats/write_dlkcat_input.py b/src/geckopy/gather_kcats/write_dlkcat_input.py index fc17cc9..2cd7737 100644 --- a/src/geckopy/gather_kcats/write_dlkcat_input.py +++ b/src/geckopy/gather_kcats/write_dlkcat_input.py @@ -82,7 +82,9 @@ def write_dlkcat_input( ``load_dlkcat_ignore_lists``). ec_rxns Optional iterable of reaction IDs (must be in - ``model.ec.rxns``). ``None`` means all. + ``model.ec.rxns``). ``None`` means all; an empty iterable selects + no reaction and writes an empty file. The rows follow the order of + ``model.ec.rxns``, whatever the order of ``ec_rxns``. only_with_smiles When True (default), drop rows whose substrate has no SMILES. When False, write ``"None"`` in the SMILES column @@ -150,7 +152,7 @@ def extract_enzyme_substrate_pairs( Loaded ``DLKcatIgnoreLists``. ec_rxns Optional iterable of reaction IDs (must be in ``model.ec.rxns``). - ``None`` means all. + ``None`` means all; the rows follow the order of ``model.ec.rxns``. only_with_smiles Drop rows with no SMILES when True; otherwise write ``"None"``. @@ -172,7 +174,7 @@ def extract_enzyme_substrate_pairs( f"{len(unknown)} reaction ID(s) in ec_rxns are not present " f"in model.ec.rxns (examples: {preview})" ) - ec_rxn_indices = [index_by_id[rid] for rid in ec_rxns_list] + ec_rxn_indices = sorted({index_by_id[rid] for rid in ec_rxns_list}) if not ec_rxn_indices: return pd.DataFrame(columns=_PAIR_COLUMNS) diff --git a/tests/test_kcat_contract.py b/tests/test_kcat_contract.py index ef145aa..a3f5506 100644 --- a/tests/test_kcat_contract.py +++ b/tests/test_kcat_contract.py @@ -69,3 +69,40 @@ def test_read_kcat_source_quoted_rows(tmp_path): file.write_text("\n".join(head[:2] + ['"R3\tP4\t25\todd']) + "\n", encoding="utf-8") with pytest.raises(ValueError, match="no closing one"): read_kcat_source(file) + + +def test_read_dlkcat_output_kcat_fields(tmp_path): + """A kcat field that is not a finite number drops the row before the reactions + are checked; a file without a number raises.""" + model, _, _ = _project(tmp_path) + file = tmp_path / "DLKcat_out.tsv" + head = "rxns\tgenes\tsubstrates\tsmiles\tsequence\tkcat" + + def write(*rows): + file.write_text("\n".join([head] + [f"{r}\tG1\t{s}\tsmiles\tseq\t{k}" for r, s, k in rows]) + "\n", + encoding="utf-8") + + write(("R2_EXP_1", "m1", "12.5"), ("R2_EXP_2", "m1", "inf"), ("R2_EXP_2", "m1", "1.2.3")) + out = read_dlkcat_output(model, file) + assert list(out["rxn_id"]) == ["R2_EXP_1"] and list(out["kcat"]) == [12.5] + write(("R2_EXP_1", "m1", "None"), ("R2_EXP_2", "m1", "None")) + with pytest.raises(ValueError, match="no numeric kcat"): + read_dlkcat_output(model, file) + write(("R2_EXP_1", "m1", "12.5"), ("NOTAREACTION_EXP_1", "m1", "NA")) + assert list(read_dlkcat_output(model, file)["rxn_id"]) == ["R2_EXP_1"] + write(("R2_EXP_1", "m1", "12.5"), ("NOTAREACTION_EXP_1", "m1", "7")) + with pytest.raises(ValueError, match="not present in model.ec.rxns"): + read_dlkcat_output(model, file) + + +def test_write_dlkcat_input_selection(tmp_path): + """An empty selection writes an empty file; the rows follow model.ec.rxns.""" + model, _, _ = _project(tmp_path) + ignore = load_dlkcat_ignore_lists() + target = tmp_path / "DLKcat_in.tsv" + df = write_dlkcat_input(model, target, ignore, ec_rxns=[], overwrite=True, only_with_smiles=False) + assert df.empty and target.read_text(encoding="utf-8") == "" + df = write_dlkcat_input(model, target, ignore, ec_rxns=["R5", "R3", "R5"], overwrite=True, + only_with_smiles=False) + wanted = [r for r in model.ec.rxns if r in {"R3", "R5"}] + assert list(dict.fromkeys(df["rxn_id"])) == wanted From 90e6dd1dfc2fe194e69205c511ca5b83cd9d2a11 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:29 +0100 Subject: [PATCH 14/27] ECTree raises for k = NaN; kcats_from_ec_tree takes an empty ec_rxns as no reaction --- CHANGELOG.md | 2 ++ src/geckopy/gather_kcats/ec_tree.py | 4 ++-- src/geckopy/gather_kcats/kcats_from.py | 9 +++++---- tests/test_kcat_contract.py | 23 +++++++++++++++++++++++ 4 files changed, 32 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8156adb..eb1e184 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,8 @@ * `read_kcat_source` and `load_kcat_source` skip a row whose fields are all empty quoted fields (`""""`), as they skip a row of tabs, and raise `ValueError` for a field that opens a double quote and does not close it. Before, the first was reported as an invalid row and the second was read as a row with a kcat that is not a number. * `read_dlkcat_output` drops a row whose kcat is infinite (`inf`), as it drops `NA` and `None`, so the kcat list holds finite numbers only. * `write_dlkcat_input` and `extract_enzyme_substrate_pairs` write the rows in the order of `model.ec.rxns` whatever the order of `ec_rxns`, and a reaction named twice once, as GECKO's `writeDLKcatInput` does with its logical vector. +* `ECTree` raises `ValueError` for `k = NaN`, as GECKO's `ecTreeEstimate` does; it returned NaN estimates. +* `kcats_from_ec_tree` with an empty `ec_rxns` returns no kcats instead of raising, as an empty selection means no reaction, and counts a reaction named twice once in the `ec_rxns` parameter of the source file. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/gather_kcats/ec_tree.py b/src/geckopy/gather_kcats/ec_tree.py index a5fce75..a58dd6f 100644 --- a/src/geckopy/gather_kcats/ec_tree.py +++ b/src/geckopy/gather_kcats/ec_tree.py @@ -21,8 +21,8 @@ class ECTree: """EC-tree estimates from a BRENDA measurement table (``ec``, ``log10_kcat``).""" def __init__(self, measurements: pd.DataFrame, k: float = 5.0): - if k <= 0: - raise ValueError("k must be positive") + if not k > 0: + raise ValueError("k must be a positive number") self.k = float(k) ec = measurements["ec"].astype(str) logk = measurements["log10_kcat"].to_numpy(dtype=float) diff --git a/src/geckopy/gather_kcats/kcats_from.py b/src/geckopy/gather_kcats/kcats_from.py index f338343..0e90894 100644 --- a/src/geckopy/gather_kcats/kcats_from.py +++ b/src/geckopy/gather_kcats/kcats_from.py @@ -191,8 +191,9 @@ def kcats_from_ec_tree( Shrinkage constant: the number of measurements at which an EC's own median and its class count equally. ec_rxns - ``model.ec.rxns`` to estimate kcats for (default: all). With the - default ``save_as`` the file then only holds those reactions. + ``model.ec.rxns`` to estimate kcats for (default: all). An empty + list selects no reaction. With the default ``save_as`` the file then + only holds those reactions. measurements BRENDA measurement table (default: ``measurements.tsv`` of the adapter's BRENDA folder, see @@ -209,7 +210,7 @@ def kcats_from_ec_tree( (``n``) and the deepest EC level with measurements (``depth``). """ adapter = _adapter(model, adapter, "kcats_from_ec_tree needs the BRENDA folder and project folder") - ec_rxns = _rxn_subset(ec_rxns) + ec_rxns = None if ec_rxns is None else list(ec_rxns) if save: path = kcat_source_file(save_as, model=model, adapter=adapter) # fails before estimating _warn_subset(ec_rxns, model, path) @@ -232,7 +233,7 @@ def kcats_from_ec_tree( if save: parameters = {"k": k} if ec_rxns is not None: - parameters["ec_rxns"] = len(ec_rxns) + parameters["ec_rxns"] = len(set(ec_rxns)) path = save_kcat_source( model, kcat_list, save_as, adapter=adapter, generator="kcats_from_ec_tree", parameters=parameters, data=measurements.attrs.get("release", ""), diff --git a/tests/test_kcat_contract.py b/tests/test_kcat_contract.py index a3f5506..1921dcb 100644 --- a/tests/test_kcat_contract.py +++ b/tests/test_kcat_contract.py @@ -106,3 +106,26 @@ def test_write_dlkcat_input_selection(tmp_path): only_with_smiles=False) wanted = [r for r in model.ec.rxns if r in {"R3", "R5"}] assert list(dict.fromkeys(df["rxn_id"])) == wanted + + +def test_ec_tree_rejects_nan_k(): + table = pd.DataFrame({"ec": ["1.1.1.1", "1.1.1.1"], "log10_kcat": [1.0, 2.0]}) + with pytest.raises(ValueError, match="positive"): + ECTree(table, k=float("nan")) + + +def test_kcats_from_ec_tree_subset(tmp_path): + """An empty ec_rxns selects no reaction; a reaction named twice gives one row.""" + model, adapter, _ = _project(tmp_path) + ecs = sorted({e.strip() for c in model.ec.eccodes for e in str(c or "").split(";") if e.strip()}) + table = pd.DataFrame([(ec, float(i + 1)) for i, ec in enumerate(ecs) for _ in range(2)], + columns=["ec", "log10_kcat"]) + table["kcat"] = 10.0 ** table["log10_kcat"] + table.attrs["release"] = "BRENDA release test" + with_ec = [r for r, c in zip(model.ec.rxns, model.ec.eccodes) if str(c or "").strip()] + assert len(with_ec) >= 2 + none = kcats_from_ec_tree(model, measurements=table, adapter=adapter, ec_rxns=[], save=False) + assert none.empty + twice = kcats_from_ec_tree(model, measurements=table, adapter=adapter, save=False, + ec_rxns=[with_ec[0], with_ec[1], with_ec[1]]) + assert list(twice["rxn_id"]) == with_ec[:2] From 594e1633487ebd82f7b26c7ea248fcf2bd8c168b Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 19:44:29 +0100 Subject: [PATCH 15/27] assign_standard_kcat leaves out reactions without a substrate or without a product as exchange reactions --- CHANGELOG.md | 1 + src/geckopy/gather_kcats/get_standard_kcat.py | 17 +++++++++++++---- tests/test_kcat_contract.py | 12 ++++++++++++ 3 files changed, 26 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eb1e184..a250e10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ * `write_dlkcat_input` and `extract_enzyme_substrate_pairs` write the rows in the order of `model.ec.rxns` whatever the order of `ec_rxns`, and a reaction named twice once, as GECKO's `writeDLKcatInput` does with its logical vector. * `ECTree` raises `ValueError` for `k = NaN`, as GECKO's `ecTreeEstimate` does; it returned NaN estimates. * `kcats_from_ec_tree` with an empty `ec_rxns` returns no kcats instead of raising, as an empty selection means no reaction, and counts a reaction named twice once in the `ec_rxns` parameter of the source file. +* `assign_standard_kcat` leaves out every reaction without a substrate or without a product (RAVEN's `getExchangeRxns`), such as `A + B =>`; it left out single-metabolite boundary reactions only. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. diff --git a/src/geckopy/gather_kcats/get_standard_kcat.py b/src/geckopy/gather_kcats/get_standard_kcat.py index b57cf13..d8760f0 100644 --- a/src/geckopy/gather_kcats/get_standard_kcat.py +++ b/src/geckopy/gather_kcats/get_standard_kcat.py @@ -57,9 +57,10 @@ def assign_standard_kcat( (identified by ``model.ec.source == "standard"``) are stripped before re-applying. - Filtering: exchange, transport (same metabolite name across - compartments), pseudo (name contains "pseudoreaction"), SLIME - (name contains "SLIME rxn"), spontaneous (via + Filtering: exchange (a reaction without a substrate or without a + product), transport (same metabolite name across compartments), pseudo + (name contains "pseudoreaction", in any letter case), SLIME (name + contains "SLIME rxn", in any letter case), spontaneous (via ``adapter.get_spontaneous_reactions(model)``), and reactions listed in ``adapter.params.path/data/pseudoRxns.tsv`` are excluded from receiving the standard pseudoenzyme. @@ -237,6 +238,14 @@ def _find_reactions_missing_enzyme(model: "EcModel") -> list[str]: return no_gpr + has_gpr_no_ec +def _is_exchange(rxn: cobra.Reaction) -> bool: + """Whether a reaction only consumes or only produces metabolites, as + RAVEN's ``getExchangeRxns`` counts it: ``A + B <=>`` is an exchange + reaction as well as ``A <=>``.""" + coefficients = list(rxn.metabolites.values()) + return all(c <= 0 for c in coefficients) or all(c >= 0 for c in coefficients) + + def _detect_transport_reactions(model: "EcModel") -> set[str]: """A transport reaction has the same metabolite name appearing in different compartments.""" @@ -257,7 +266,7 @@ def _classify_reactions_to_ignore( ignore: set[str] = set(custom_rxns) for rxn in model.reactions: - if rxn.boundary: + if _is_exchange(rxn): ignore.add(rxn.id) name_lower = (rxn.name or "").lower() if "pseudoreaction" in name_lower or "slime rxn" in name_lower: diff --git a/tests/test_kcat_contract.py b/tests/test_kcat_contract.py index 1921dcb..f2e260b 100644 --- a/tests/test_kcat_contract.py +++ b/tests/test_kcat_contract.py @@ -129,3 +129,15 @@ def test_kcats_from_ec_tree_subset(tmp_path): twice = kcats_from_ec_tree(model, measurements=table, adapter=adapter, save=False, ec_rxns=[with_ec[0], with_ec[1], with_ec[1]]) assert list(twice["rxn_id"]) == with_ec[:2] + + +def test_standard_kcat_exchange_and_pseudoreactions(tmp_path): + """R7 consumes two metabolites and produces none: an exchange reaction. R8 is + named like a pseudoreaction in another letter case. R6 is an ordinary reaction + without a GPR.""" + model, _, _ = _project(tmp_path, extra_reactions=True) + model.ec.kcat[:] = 10.0 + assign_standard_kcat(model, _uniprot()) + assert "R6" in model.ec.rxns + assert "R7" not in model.ec.rxns + assert "R8" not in model.ec.rxns From a4d73c5739b63ce8dd57d8819d9183a793de8416 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 21:55:32 +0100 Subject: [PATCH 16/27] save_ec_model keeps an existing metaData description The default description is written only when the model has none, as saveEcModel does. --- CHANGELOG.md | 1 + src/geckopy/utilities/save_ec_model.py | 7 ++++--- tests/test_save_ec_model.py | 8 ++++++++ 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f816a95..193062a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `save_ec_model` keeps the description of a model that has one (`metaData` of the file it was loaded from) and writes `Enzyme-constrained model of ` only when there is none; it used to overwrite it. GECKO's `saveEcModel` does the same. * `best_parsimonious` skips points with a NaN `rmse`, as MATLAB's `bestParsimonious` does with `min(rmses)`: a NaN at the first point of the frontier made every comparison false and raised `ValueError`. It still raises when every `rmse` is NaN. * `screen_kcat_leverage` keeps rows of equal rank in the order of the tie groups (a stable sort), as MATLAB's `screenKcatLeverage` does, and `compute_excarbon` warns when no metabolite of the model has a formula (every exchange then gets the carbon weight 1). * `cmaes_kcat_tuning` runs `max_generations` generations on a flat objective instead of stopping after the first (the flat-fitness criteria of `cma` do not end the search), as MATLAB's `cmaesKcatTuning` does. diff --git a/src/geckopy/utilities/save_ec_model.py b/src/geckopy/utilities/save_ec_model.py index f49c97a..7e94bb6 100644 --- a/src/geckopy/utilities/save_ec_model.py +++ b/src/geckopy/utilities/save_ec_model.py @@ -139,8 +139,8 @@ def _build_metadata(model: EcModel) -> dict: format the file follows (as ``project_format`` in ``model_adapter.toml``, written as text as GECKO's ``saveEcModel`` writes it); ``generator``, the release that wrote the file; and a - description string of the form ``"Enzyme-constrained model of "``. - GECKO's ``saveEcModel`` writes the same ``gecko_format`` and + description string of the form ``"Enzyme-constrained model of "``, + written only when the model has no description of its own. GECKO's ``saveEcModel`` writes the same ``gecko_format`` and ``generator``. """ try: @@ -152,5 +152,6 @@ def _build_metadata(model: EcModel) -> dict: meta["date"] = meta.get("date") or date.today().isoformat() meta["gecko_format"] = "4" meta["generator"] = f"geckopy {geckopy_version}" - meta["description"] = f"Enzyme-constrained model of {model.id}" + if not meta.get("description"): + meta["description"] = f"Enzyme-constrained model of {model.id}" return meta diff --git a/tests/test_save_ec_model.py b/tests/test_save_ec_model.py index 9c89601..d90195e 100644 --- a/tests/test_save_ec_model.py +++ b/tests/test_save_ec_model.py @@ -201,6 +201,14 @@ def test_doc_top_level_keys(tmp_path): }.issubset(doc.keys()) +def test_existing_description_is_kept(tmp_path): + adapter = _adapter(tmp_path) + model = _build_full_ec_model(adapter) + model.notes["metaData"] = {"description": "my own description"} + path = save_ec_model(model) + assert _read_yaml(path)["metaData"]["description"] == "my own description" + + def test_metadata_includes_date_version_description(tmp_path): adapter = _adapter(tmp_path) model = _build_full_ec_model(adapter) From 739f020b99288b6787211de8c4c4d58d4d2dd590 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:00:56 +0100 Subject: [PATCH 17/27] review_assignment recognises the sources custom and standard without regard to case --- src/geckopy/gather_kcats/review.py | 7 +++++-- tests/test_review_public_api.py | 11 +++++++++++ 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/geckopy/gather_kcats/review.py b/src/geckopy/gather_kcats/review.py index 177ef90..1ebbf24 100644 --- a/src/geckopy/gather_kcats/review.py +++ b/src/geckopy/gather_kcats/review.py @@ -113,6 +113,9 @@ def review_assignment( flagged rows are a minority of the model's leverage. include Names from :data:`OPTIONAL_CHECKS` to report as well. + sources + The source labels ``custom`` and ``standard`` are recognised + without regard to letter case. Returns ------- @@ -176,7 +179,7 @@ def review_assignment( if kcats[i] < slow or kcats[i] > fast: checks.append("magnitude") detail.append(f"{kcats[i]:g} 1/s is outside {slow:g}-{fast:g}") - if (str(sources[i]) == "custom" and st is not None + if (str(sources[i]).lower() == "custom" and st is not None and round(float(kcats[i]), 6) in st.values): checks.append("custom-duplicate") detail.append("custom value equals a database value for its own EC") @@ -184,7 +187,7 @@ def review_assignment( if "tied-isozymes" in include: checks.append("tied-isozymes") detail.append(f"{len(idx)} isozyme copies share this value") - if str(sources[i]) == "standard": + if str(sources[i]).lower() == "standard": if "standard-fallback" not in include: continue checks.append("standard-fallback") diff --git a/tests/test_review_public_api.py b/tests/test_review_public_api.py index fe704a4..90cab47 100644 --- a/tests/test_review_public_api.py +++ b/tests/test_review_public_api.py @@ -40,3 +40,14 @@ def test_ec_stats_keep_the_case_of_the_ec_number(tmp_path): ["R1"], [40.0], ["brenda"], ["3.5.1.n3"], ec_stats=stats, ) assert [f.checks for f in findings] == [("ec-maximum",)] + + +def test_review_recognises_custom_and_standard_without_regard_to_case(): + stats = {"1.1.1.1": geckopy.EcStats(n_kcat=1, kcat_max=2.0, kcat_median=2.0, values=frozenset({2.0}))} + found = geckopy.review_assignment( + ["R1", "R2"], [2.0, 20000.0], ["Custom", "STANDARD"], ["1.1.1.1", ""], ec_stats=stats, + include=("standard-fallback",), + ) + by_rxn = {f.rxn_id: f.checks for f in found} + assert "custom-duplicate" in by_rxn["R1"] + assert "standard-fallback" in by_rxn["R2"] From d42bdcebacf6311bd1f85cf16ffa8fe98b998fb6 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:00:57 +0100 Subject: [PATCH 18/27] find_met_smiles skips metabolites without a name with a warning --- src/geckopy/databases/pubchem.py | 22 +++++++++++++++++----- tests/test_pubchem.py | 19 +++++++++++++++++++ 2 files changed, 36 insertions(+), 5 deletions(-) diff --git a/src/geckopy/databases/pubchem.py b/src/geckopy/databases/pubchem.py index edcee5a..fc9ef0d 100644 --- a/src/geckopy/databases/pubchem.py +++ b/src/geckopy/databases/pubchem.py @@ -73,9 +73,11 @@ def find_met_smiles( For every metabolite in the model: 1. Skip if the name starts with ``prot_`` (enzyme pseudometabolite). - 2. Skip if ``annotation['smiles']`` is already set. - 3. Look up the metabolite name in the local TSV cache. - 4. If still not found, query PubChem; cache the result. + 2. Skip, with a warning, if the metabolite has no name. + 3. Skip if ``annotation['smiles']`` is already set. + 4. Look up the metabolite name in the local TSV cache. When the cache + holds a name more than once, the last row is used. + 5. If still not found, query PubChem; cache the result. The same metabolite name may appear on many ``cobra.Metabolite`` instances (different compartments). Lookup is done once per unique @@ -105,14 +107,24 @@ def find_met_smiles( # Build map: metabolite name -> list of metabolite objects sharing it. name_to_mets: dict[str, list] = {} + unnamed = [] for met in model.metabolites: - name = met.name or met.id - if not name or _PROT_PATTERN.match(name): + name = met.name + if not name: + unnamed.append(met.id) + continue + if _PROT_PATTERN.match(name): continue if met.annotation.get("smiles"): continue name_to_mets.setdefault(name, []).append(met) + if unnamed: + logger.warning( + "find_met_smiles: %d metabolite(s) have no name and get no SMILES: %s", + len(unnamed), ", ".join(unnamed), + ) + if not name_to_mets: logger.info("find_met_smiles: nothing to look up.") return diff --git a/tests/test_pubchem.py b/tests/test_pubchem.py index 4ba33d2..3a4e41d 100644 --- a/tests/test_pubchem.py +++ b/tests/test_pubchem.py @@ -368,3 +368,22 @@ def test_met_smiles_accepts_the_shape_a_loaded_model_carries(): del met.annotation["smiles"] assert met_smiles(met) == "" + + +def test_duplicate_cache_name_uses_the_last_row(monkeypatch, tmp_path): + cache_path = tmp_path / "smilesDB.tsv" + cache_path.write_text("glucose\tFIRST\nglucose\tLAST\n") + model = _build_model([("glc_c", "glucose")]) + find_met_smiles(model, cache_path=cache_path) + assert _smiles(model, "glc_c") == "LAST" + + +def test_metabolite_without_name_is_skipped_with_a_warning(monkeypatch, tmp_path, caplog): + cache_path = tmp_path / "smilesDB.tsv" + cache_path.write_text("m1_c\tFROM_ID\nglucose\tCCO\n") + model = _build_model([("m1_c", ""), ("glc_c", "glucose")]) + with caplog.at_level(logging.WARNING): + find_met_smiles(model, cache_path=cache_path) + assert _smiles(model, "m1_c") == "" + assert _smiles(model, "glc_c") == "CCO" + assert any("m1_c" in r.getMessage() and "no name" in r.getMessage() for r in caplog.records) From 18397ec1a9371e530997431be06664e3e400ee3a Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:01:05 +0100 Subject: [PATCH 19/27] Changelog: source case in review_assignment, unnamed metabolites in find_met_smiles --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 72fecf4..27e653b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,8 @@ ## Unreleased +* `review_assignment` recognises the source labels `custom` and `standard` without regard to letter case, as MATLAB's `reviewKcatAssignment` does: a source `Custom` is checked for a duplicate of a database value and `STANDARD` is the fallback. +* `find_met_smiles` skips a metabolite without a name with a logged warning and gives it no SMILES; it used to look the metabolite up by its id. When the cache file holds a name more than once, the last row is used (as before, now stated in the docstring). * `review_assignment`, `EcStats` and `ec_stats_from_brenda` are exported from `geckopy.gather_kcats` and `geckopy`, and listed in `__all__`. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. From 7b25cd994df6d1d08d3cd81eafe66c2c9b03512b Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:02:41 +0100 Subject: [PATCH 20/27] pfba_enzymes multiplies a given fraction_of_optimum by 0.999999 --- src/geckopy/utilities/pfba_enzymes.py | 27 ++++++++++++++++----------- tests/test_pfba_enzymes.py | 13 ++++++++++--- 2 files changed, 26 insertions(+), 14 deletions(-) diff --git a/src/geckopy/utilities/pfba_enzymes.py b/src/geckopy/utilities/pfba_enzymes.py index 39581f7..66be1c9 100644 --- a/src/geckopy/utilities/pfba_enzymes.py +++ b/src/geckopy/utilities/pfba_enzymes.py @@ -33,22 +33,23 @@ _PFBA_OBJ_NAME = "_pfba_enzymes_objective" -# RAVEN's solveLP fixes the objective at this fraction of its optimum in the flux-minimising -# stage (getPfbaEnzymes calls it with minFlux 1). A solver reports the optimum to within its -# feasibility tolerance, so a second solve that demands the full value can be infeasible. -_PFBA_FRACTION = 0.999999 +# RAVEN's solveLP multiplies the fraction of the optimum by this margin in the +# flux-minimising stage (getPfbaEnzymes calls it with minFlux 1). A solver reports the +# optimum to within its feasibility tolerance, so a second solve that demands the full +# value can be infeasible. +_PFBA_MARGIN = 0.999999 def pfba_enzymes( model: "EcModel", objective: Union[cobra.Reaction, str, None] = None, *, - fraction_of_optimum: float = _PFBA_FRACTION, + fraction_of_optimum: float = 1.0, ) -> cobra.Solution: """Parsimonious FBA minimising the L1 norm of ``usage_prot_*`` fluxes. Sets the objective (if given), fixes it at - ``fraction_of_optimum * optimal_value``, then minimises + ``fraction_of_optimum * 0.999999 * optimal_value``, then minimises ``sum(usage_prot_*.forward_variable)``. Usage reactions are forward-only in the GECKO 3 layout (``lb=0``); if a user has flipped one to allow reverse flux, its ``reverse_variable`` is @@ -69,10 +70,12 @@ def pfba_enzymes( used as-is. fraction_of_optimum Fraction of the objective's optimal value to fix as a - constraint before minimising enzyme usage. The default, - 0.999999, is the fraction MATLAB's ``getPfbaEnzymes`` uses - (RAVEN's ``solveLP``); a lower value trades some objective - value for potentially lower enzyme usage. + constraint before minimising enzyme usage (default 1.0). The + fraction is multiplied by 0.999999, the margin of RAVEN's + ``solveLP`` that MATLAB's ``getPfbaEnzymes`` applies, so the + default fixes the objective at 0.999999 of its optimum and + 0.9 fixes it at 0.8999991. A lower value trades some + objective value for potentially lower enzyme usage. Returns ------- @@ -105,7 +108,9 @@ def pfba_enzymes( f"Objective is already {_PFBA_OBJ_NAME}; " "pfba_enzymes appears to be nested." ) - cobra.util.fix_objective_as_constraint(m, fraction=fraction_of_optimum) + cobra.util.fix_objective_as_constraint( + m, fraction=fraction_of_optimum * _PFBA_MARGIN, + ) usage_rxns = [r for r in m.reactions if r.id.startswith(USAGE_PREFIX)] variables = [] diff --git a/tests/test_pfba_enzymes.py b/tests/test_pfba_enzymes.py index fd5aa2f..2259839 100644 --- a/tests/test_pfba_enzymes.py +++ b/tests/test_pfba_enzymes.py @@ -45,11 +45,18 @@ def test_returns_solution(): assert sol.status == "optimal" -def test_objective_preserved_at_full_fraction(): +def test_full_fraction_keeps_0_999999_of_the_optimum(): model = _ectestgem_ec_model_with_kcats() fba_optimum = model.slim_optimize() sol = pfba_enzymes(model, fraction_of_optimum=1.0) - assert sol.fluxes["R3"] == pytest.approx(fba_optimum, abs=1e-6) + assert sol.fluxes["R3"] / fba_optimum == pytest.approx(0.999999, abs=2e-7) + + +def test_given_fraction_is_multiplied_by_the_margin(): + model = _ectestgem_ec_model_with_kcats() + fba_optimum = model.slim_optimize() + sol = pfba_enzymes(model, fraction_of_optimum=0.9) + assert sol.fluxes["R3"] / fba_optimum == pytest.approx(0.9 * 0.999999, abs=2e-7) def test_default_fixes_the_objective_at_0_999999_of_the_optimum(): @@ -65,7 +72,7 @@ def test_objective_preserved_at_partial_fraction(): model = _ectestgem_ec_model_with_kcats() fba_optimum = model.slim_optimize() sol = pfba_enzymes(model, fraction_of_optimum=0.5) - assert sol.fluxes["R3"] >= 0.5 * fba_optimum - 1e-6 + assert sol.fluxes["R3"] >= 0.5 * 0.999999 * fba_optimum - 1e-6 def test_total_usage_le_plain_fba(): From ce7ad7dac35b300e72a6245cc2d05434ad4ed2dc Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:02:41 +0100 Subject: [PATCH 21/27] apply_kcat_constraints takes cobra reaction ids and ec ids on a gecko-light model --- src/geckopy/ec_model/pipeline/apply_kcat.py | 43 ++++++++++++++++----- tests/test_apply_kcat.py | 23 ++++++++++- 2 files changed, 56 insertions(+), 10 deletions(-) diff --git a/src/geckopy/ec_model/pipeline/apply_kcat.py b/src/geckopy/ec_model/pipeline/apply_kcat.py index 989251c..221a8d1 100644 --- a/src/geckopy/ec_model/pipeline/apply_kcat.py +++ b/src/geckopy/ec_model/pipeline/apply_kcat.py @@ -65,7 +65,9 @@ def apply_kcat_constraints( An EcModel produced by ``make_ec_model``, with at least some kcat values assigned to ``model.ec.kcat``. Mutated in place. update_rxns - Reaction IDs (from ``model.ec.rxns``) to update. ``None`` or an + Reaction IDs (from ``model.ec.rxns``) to update. A gecko-light + model also takes the IDs of its cobra reactions (see + ``_apply_kcat_constraints_light``). ``None`` or an empty sequence updates all entries, as ``applyKcatConstraints`` does with an empty ``updateRxns``; a caller that builds the list from a filter must skip the call when the filter finds nothing. @@ -216,10 +218,13 @@ def _apply_kcat_constraints_light( model A gecko-light EcModel with ``ec.gecko_light is True``. update_rxns - ec.rxns IDs (with the ``###_`` prefix) to update. ``None`` or an - empty sequence means update every cobra reaction backed by at - least one ec.rxns row. A reaction is updated as a whole: naming one of its - isozyme rows recomputes the cheapest cost over all its rows. + IDs to update: ec.rxns IDs (with the ``###_`` prefix) or cobra + reaction IDs. ``None`` or an empty sequence means update every + cobra reaction backed by at least one ec.rxns row. A reaction + is updated as a whole: naming one of its isozyme rows, or the + cobra reaction, recomputes the cheapest cost over all its rows. + An ID that is both an ec.rxns ID and a cobra reaction ID, or + neither, raises ``ValueError``. """ ec_rxn_ids = model.ec.rxns @@ -230,14 +235,34 @@ def _apply_kcat_constraints_light( selected_idx = np.arange(len(ec_rxn_ids)) else: id_to_idx = {rxn_id: i for i, rxn_id in enumerate(ec_rxn_ids)} - unknown = [r for r in update_rxns if r not in id_to_idx] + rows_of_cobra: dict[str, list[int]] = defaultdict(list) + for i, rxn_id in enumerate(ec_rxn_ids): + rows_of_cobra[split_light_rxn_id(rxn_id)[1]].append(i) + cobra_ids = {r.id for r in model.reactions} + ambiguous = [r for r in update_rxns if r in id_to_idx and r in cobra_ids] + if ambiguous: + raise ValueError( + f"update_rxns contains IDs that are both ec.rxns IDs and " + f"reaction IDs: {ambiguous[:5]}" + ) + unknown = [ + r for r in update_rxns + if r not in id_to_idx and r not in rows_of_cobra + ] if unknown: raise ValueError( - f"update_rxns contains IDs not present in ec.rxns: " - f"{unknown[:5]}" + f"update_rxns contains IDs present neither in ec.rxns nor " + f"as reactions of the model: {unknown[:5]}" ) selected_idx = np.array( - [id_to_idx[r] for r in update_rxns], dtype=int, + [ + i + for r in update_rxns + for i in ( + [id_to_idx[r]] if r in id_to_idx else rows_of_cobra[r] + ) + ], + dtype=int, ) if selected_idx.size == 0: diff --git a/tests/test_apply_kcat.py b/tests/test_apply_kcat.py index 73e87f5..c18e1f1 100644 --- a/tests/test_apply_kcat.py +++ b/tests/test_apply_kcat.py @@ -324,10 +324,31 @@ def test_light_setting_kcat_to_zero_and_reapplying_clears_old_constraint(): def test_light_update_rxns_unknown_id_raises(): ec = _ectestgem_light_model() - with pytest.raises(ValueError, match="not present in ec.rxns"): + with pytest.raises(ValueError, match="neither in ec.rxns nor"): apply_kcat_constraints(ec, update_rxns=["999_nonsense"]) +def test_light_update_rxns_takes_cobra_ids_and_ec_ids(): + ec = _ectestgem_light_model() + ec.ec.kcat[ec.ec.rxns.index("001_R3")] = 10.0 + r3 = ec.reactions.get_by_id("R3") + apply_kcat_constraints(ec, update_rxns=["R3"]) + coef = _get_s_coef(r3, "prot_pool") + assert coef < 0 + apply_kcat_constraints(ec, update_rxns=["001_R5"]) + assert _get_s_coef(r3, "prot_pool") == coef + ec.ec.kcat[ec.ec.rxns.index("001_R3")] = 0.0 + apply_kcat_constraints(ec, update_rxns=["001_R3"]) + assert _get_s_coef(r3, "prot_pool") == 0.0 + + +def test_light_update_rxns_id_of_both_kinds_raises(): + ec = _ectestgem_light_model() + ec.reactions.get_by_id("R5").id = "001_R3" + with pytest.raises(ValueError, match="both"): + apply_kcat_constraints(ec, update_rxns=["001_R3"]) + + def test_light_warns_when_no_valid_kcats(): """When every selected reaction has kcat=0 across all isozymes, emit the same warning the full branch does and leave the model unchanged.""" From 44a78153737abcb68cb940326604fb974bf842c5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:02:41 +0100 Subject: [PATCH 22/27] constrain_enz_concs resets the lower bound and the pool coefficient of a usage reaction --- src/geckopy/limit_proteins/constrain_enz_concs.py | 15 ++++++++++----- tests/test_constrain_enz_concs.py | 15 +++++++++++++++ 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/src/geckopy/limit_proteins/constrain_enz_concs.py b/src/geckopy/limit_proteins/constrain_enz_concs.py index eb0591a..a1e1442 100644 --- a/src/geckopy/limit_proteins/constrain_enz_concs.py +++ b/src/geckopy/limit_proteins/constrain_enz_concs.py @@ -35,9 +35,10 @@ def constrain_enz_concs( ``usage_prot_`` reaction's upper bound is set to that concentration. Enzymes with NaN concentrations have their usage reaction reset to the default ``upper_bound = 1000``, drawing - freely from ``prot_pool``. Only the reaction's upper bound is - changed; its stoichiometry (as set by - ``add_protein_usage_reactions``) is left untouched. + freely from ``prot_pool``. Each processed usage reaction is first + reset to ``lower_bound = 0``, ``upper_bound = 1000`` and a + ``prot_pool`` coefficient of -1, so a prior lower bound above the + new upper bound does not conflict with it. Ported from GECKO MATLAB: src/geckomat/limit_proteins/constrainEnzConcs.m. @@ -50,7 +51,8 @@ def constrain_enz_concs( machinery already in place. Mutated in place. remove_constraints When True, all usage reactions are reset to their default - ``upper_bound = 1000``, regardless of ``ec.concs``. + ``upper_bound = 1000`` (with ``lower_bound = 0`` and a + ``prot_pool`` coefficient of -1), regardless of ``ec.concs``. ``ec.concs`` itself is left unchanged. restrict_to If supplied, only the listed enzyme uniprot IDs are @@ -107,9 +109,12 @@ def constrain_enz_concs( + ", ".join(negative) + "." ) + pool_met = model.metabolites.get_by_id(POOL_ID) for i, enz in index_list: rxn = model.reactions.get_by_id(f"{USAGE_PREFIX}{enz}") - rxn.upper_bound = 1000.0 + rxn.bounds = (0.0, 1000.0) + if rxn.metabolites.get(pool_met) != -1.0: + rxn.add_metabolites({pool_met: -1.0}, combine=False) if remove_constraints: continue conc = model.ec.concs[i] diff --git a/tests/test_constrain_enz_concs.py b/tests/test_constrain_enz_concs.py index 45908ba..25045d1 100644 --- a/tests/test_constrain_enz_concs.py +++ b/tests/test_constrain_enz_concs.py @@ -233,3 +233,18 @@ def test_negative_conc_is_not_read_when_removing_constraints(): model = _ec_model_with_pool(["P1"], concs=[-0.5]) constrain_enz_concs(model, remove_constraints=True) assert model.reactions.get_by_id("usage_prot_P1").upper_bound == 1000.0 + + +# --------------------------------------------------------------------------- # +# Prior bounds and coefficient +# --------------------------------------------------------------------------- # + +def test_prior_lower_bound_above_the_new_upper_bound_is_reset(): + model = _ec_model_with_pool(["P1"], concs=[0.01]) + rxn = model.reactions.get_by_id("usage_prot_P1") + rxn.lower_bound = 0.05 + pool = model.metabolites.get_by_id("prot_pool") + rxn.add_metabolites({pool: -2.0}, combine=False) + constrain_enz_concs(model) + assert rxn.bounds == (0.0, 0.01) + assert rxn.metabolites[pool] == -1.0 From 75dfa1546a9e67c064d58c2c6d0128dc64268e99 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 22:13:32 +0100 Subject: [PATCH 23/27] Add the CHANGELOG entries for pfba_enzymes, apply_kcat_constraints and constrain_enz_concs --- CHANGELOG.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1ebae89..d9e9e0a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,8 @@ ## Unreleased +* `apply_kcat_constraints` takes cobra reaction ids as well as ec reaction ids (with the `###_` prefix) in `update_rxns` on a gecko-light model, as GECKO's `applyKcatConstraints` does for `model.rxns` ids; naming the cobra reaction updates it as a whole. An id that is both kinds or neither raises a `ValueError`. +* `constrain_enz_concs` resets each processed usage reaction to a lower bound of 0 and a `prot_pool` coefficient of -1 before it sets the upper bound, as GECKO's `constrainEnzConcs` does; a prior lower bound above the new upper bound used to raise. * `fill_enz_concs` raises a `ValueError` for a negative abundance of a protein that is an enzyme of the model, as GECKO's `fillEnzConcs` does with the error `fillEnzConcs:negativeAbundance`; it used to skip the value with a warning. * `constrain_enz_concs` raises a `ValueError` for a negative concentration in `model.ec.concs` before it changes any bound, as GECKO's `constrainEnzConcs` does with the error `constrainEnzConcs:negativeConcentration`; it used to skip the enzyme with a warning. `remove_constraints=True` does not read the concentrations. * `load_flux_data` raises a `ValueError` for a row with a different number of cells than the header, as GECKO's `loadFluxData` does with the error `loadFluxData:raggedRow`; it used to pad the missing cells with NaN. A header `X (EX_a(e))` gives the metabolite `X` and the reaction id `EX_a(e)`: the reaction id is the text between the parenthesis that matches the final one and that one; a header with nested parentheses used to give the whole header for both. @@ -10,7 +12,7 @@ * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. * `ec_fseof` follows `ecFSEOF.m`. The target flux is fixed at each level of the scan (`equality=True` of the raven-toolbox `fseof`, which `ecFSEOF.m` sets with `setParam(model, 'eq', ...)`) instead of lower-bounded. `targets` lists the reaction targets first and then the transport targets (MATLAB's `rxnTargets` and `transportTargets`), each group sorted by score, and has a new boolean column `transport`; a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting `prot_*`. This changes the row order of `targets` (and of `amplification`, `knockout` and `gene_targets`, which are derived from it) for callers that take the top N: transport targets now come after all reaction targets. -* `pfba_enzymes` fixes the objective at 0.999999 of its optimum by default, as MATLAB's `getPfbaEnzymes` does through RAVEN's `solveLP`; the default was 1.0, so `objective / optimum` of the solution was 1.0 against MATLAB's 0.999999. Pass `fraction_of_optimum=1.0` for the old behaviour. +* `pfba_enzymes` fixes the objective at 0.999999 of its optimum by default, as MATLAB's `getPfbaEnzymes` does through RAVEN's `solveLP`; a given `fraction_of_optimum` is multiplied by 0.999999 as well, so `fraction_of_optimum=0.9` fixes the objective at 0.8999991 of its optimum and the default 1.0 at 0.999999. * `model.tolerance` reaches HiGHS on the scipy solver. optlang's scipy interface ignored every tolerance, so a model set to `1e-9` (as the Gurobi and GLPK runs of the parity scenarios are) was solved at HiGHS's `1e-7`; the flux ranges of `ec_fva` on ecYeastGEM then differed from Gurobi's by up to 0.2 percent. The feasibility tolerance becomes `primal_feasibility_tolerance` and the optimality tolerance `dual_feasibility_tolerance`; both stay at HiGHS's default until set, and survive a copy or a pickle of the model. HiGHS ends some LPs with numerical difficulties (status 4) at `1e-9`; those are solved again at its default, and a solve that still ends so has the status `numeric` instead of raising `KeyError`. * A scipy `EcModel` is copied and pickled into the sparse problem. optlang pickles a model as JSON and restores it by adding the variables and constraints to a dense problem, so every worker process of `ec_fva` (and `cobra.Model.copy`, `copy.deepcopy`) spent tens of minutes rebuilding the matrix of a genome-scale ecModel; ecYeastGEM now restores in 1.6 s and copies in 2.2 s. The solver of a repaired or converted scipy model is a `geckopy.solver_support._SparseModel`. * `load_ec_model` reads the file on glpk and converts the ecModel to scipy when scipy is cobra's default solver, as `load_conventional_gem` does, and setting `model.solver = "scipy"` on an `EcModel` that is not on scipy builds the sparse problem of `geckopy.solver_support.switch_to_scipy`. optlang's own clone grows the dense coefficient matrix row by row: loading ecYeastGEM (8001 reactions) on scipy did not finish in 25 minutes and is now a matter of seconds. The scenarios that load a genome-scale ecModel on scipy spent most of their hours there. `geckopy.solver_support.is_scipy_interface` tells whether a solver name or interface is scipy. From ae3722c5886be8dfe6ee418178082ca285f094d4 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 1 Nov 2026 10:23:23 +0100 Subject: [PATCH 24/27] geckopy.kcat_tuning.evotune exports the parsimony and isozyme-tying helpers --- CHANGELOG.md | 1 + src/geckopy/kcat_tuning/evotune/__init__.py | 6 ++++++ 2 files changed, 7 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index bbbfb09..ee330ec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `geckopy.kcat_tuning.evotune` exports `best_parsimonious`, `n_changed`, `parsimony_frontier` and `isozyme_tie_map`, the helpers behind the parsimony sweep and the tying of isozyme copies, which GECKO has as `bestParsimonious`, `nChanged`, `parsimonyFrontier` and `isozymeTieMap`. * `read_kcat_source` and `load_kcat_source` skip a row whose fields are all empty quoted fields (`""""`), as they skip a row of tabs, and raise `ValueError` for a field that opens a double quote and does not close it. Before, the first was reported as an invalid row and the second was read as a row with a kcat that is not a number. * `read_dlkcat_output` drops a row whose kcat is infinite (`inf`), as it drops `NA` and `None`, so the kcat list holds finite numbers only. * `write_dlkcat_input` and `extract_enzyme_substrate_pairs` write the rows in the order of `model.ec.rxns` whatever the order of `ec_rxns`, and a reaction named twice once, as GECKO's `writeDLKcatInput` does with its logical vector. diff --git a/src/geckopy/kcat_tuning/evotune/__init__.py b/src/geckopy/kcat_tuning/evotune/__init__.py index 2a6e00f..7243d87 100644 --- a/src/geckopy/kcat_tuning/evotune/__init__.py +++ b/src/geckopy/kcat_tuning/evotune/__init__.py @@ -25,6 +25,7 @@ recommend_target_impact_share, sweep_tunable_mask_size, ) +from .parsimony import best_parsimonious, n_changed, parsimony_frontier from .priors import ( UNLABELLED_GROUP, build_sigma0_log, @@ -32,6 +33,7 @@ classify_kcat_sources, ) from .simulate import ConditionSimResult, simulate_evotune_dataset +from .tying import isozyme_tie_map from .tuning import ( EvotuneResult, cmaes_kcat_tuning, @@ -49,13 +51,17 @@ "ConditionSimResult", "MaskSizePoint", "evotune_distance", + "best_parsimonious", "build_sigma0_log", "classify_kcat_source", "classify_kcat_sources", "cmaes_kcat_tuning", "compute_excarbon", "dataset_rmse", + "isozyme_tie_map", "load_evotune_data", + "n_changed", + "parsimony_frontier", "recommend_target_impact_share", "screen_kcat_leverage", "select_tunable_mask", From 2e2c12cef39d1d993a135532f9fd53a12a479391 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 1 Nov 2026 10:26:09 +0100 Subject: [PATCH 25/27] save_ec_model: write no enzymes entry for an ec reaction without enzymes Tests cover the round trip, the saved text and the light tutorial ecModel. The format page and the docstrings describe the optional enzymes entry. --- CHANGELOG.md | 1 + docs/yaml_format.md | 2 +- src/geckopy/utilities/load_ec_model.py | 4 +- src/geckopy/utilities/save_ec_model.py | 4 +- tests/test_load_ec_model.py | 2 +- tests/test_save_ec_model.py | 72 ++++++++++++++++++++++++++ 6 files changed, 81 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c19cb2..af20b66 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `save_ec_model` writes an ec reaction without enzymes with no `enzymes` entry, as GECKO's `saveEcModel` does (it wrote `enzymes: !!omap []`). `load_ec_model` reads a missing `enzymes` entry as no enzymes, and still reads `enzymes: !!omap` with no entries and `enzymes: !!omap []` from files written by earlier versions of GECKO, RAVEN and geckopy. * `get_conc_control_coeffs` takes the control coefficient by finite difference by default, as GECKO's `getConcControlCoeffs` does: it raises one usage bound by `fold_change` (2), re-solves and divides the rise of the optimum by the rise of the bound. The reduced-cost path, which reads the dual of the usage reaction and is unreliable at a degenerate vertex, is `method="shadow_price"`. The finite-difference path needs one solve per protein, so each iteration of `flexibilize_enz_concs`, which calls it, takes minutes on a genome-scale ecModel. * `get_subset_ec_model` works on the scipy solver. It raised `LookupError: Constraint prot_P1 ... not in solver` at the next `optimize` when it had removed two or more constraints: optlang's scipy interface calls the generic removal inside its per-constraint loop, with the whole list each time, so the second pass finds the constraints gone. The sparse scipy model removes constraints correctly, and the function removes genes, reactions and metabolites on glpk and converts the subset back to scipy, because optlang rewrites every constraint expression for each removed variable (65 minutes for the HT29-GEM subset of the light Human-GEM ecModel). `geckopy.solver_support.edit_on_fast_solver` is the context manager that does this for any structural edit of a scipy model. * `switch_to_scipy` and so `make_ec_model`, `load_conventional_gem` and `load_ec_model` on scipy create each constraint with its whole expression instead of setting the coefficients one at a time, which made the scipy interface remove and add the rows of a constraint for every call (the light Human-GEM ecModel took 2 h 17 min to build on scipy against 1.5 min on Gurobi). `apply_kcat_constraints` edits the model on glpk when it updates 100 or more ec reactions of a scipy model, for the same reason. diff --git a/docs/yaml_format.md b/docs/yaml_format.md index 283c58a..516ada8 100644 --- a/docs/yaml_format.md +++ b/docs/yaml_format.md @@ -122,7 +122,7 @@ entry here. Fields: | `source` | string | optional | Where the kcat came from (`"brenda"`, `"dlkcat"`, `"manual"`, ...). Default `""` | | `notes` | string | optional | Free-form note. Default `""` | | `eccodes` | string \| list | optional | One EC code as a string (`"1.1.2.4"`), or a list when several apply. Default `""` | -| `enzymes` | mapping `enzyme_id -> stoich` | yes | Which enzymes catalyse this reaction, and the subunit count of each. Every key must appear in `ec-enzymes` | +| `enzymes` | mapping `enzyme_id -> stoich` | optional | Which enzymes catalyse this reaction, and the subunit count of each. Every key must appear in `ec-enzymes`. A reaction without enzymes (such as a `standard` kcat row) has no `enzymes` entry; the loader also reads the `enzymes: !!omap` and `enzymes: !!omap []` forms that earlier versions of GECKO, RAVEN and raven-toolbox wrote for it | Example: diff --git a/src/geckopy/utilities/load_ec_model.py b/src/geckopy/utilities/load_ec_model.py index 18cc53a..4700a9f 100644 --- a/src/geckopy/utilities/load_ec_model.py +++ b/src/geckopy/utilities/load_ec_model.py @@ -23,7 +23,9 @@ 4. **Tolerated input**: a file that the YAML reader rejects is read again after the two forms that GECKO's ``loadEcModel`` (through RAVEN's ``readYAMLmodel``) accepts are rewritten as valid YAML: ``enzymes: !!omap`` - with no entries, as GECKO writes it for a reaction without enzymes, and a + with no entries, as earlier GECKO and RAVEN versions wrote it for a + reaction without enzymes (current writers omit the ``enzymes`` entry, which + is read as no enzymes), and a value in double quotes that holds unescaped double quotes, such as ``""text""``, which RAVEN 2 wrote for a value that has quotes of its own. The value loaded is the one MATLAB loads: the outer pair of quotes is diff --git a/src/geckopy/utilities/save_ec_model.py b/src/geckopy/utilities/save_ec_model.py index f49c97a..ccc4ad2 100644 --- a/src/geckopy/utilities/save_ec_model.py +++ b/src/geckopy/utilities/save_ec_model.py @@ -73,7 +73,9 @@ def save_ec_model( Returns ------- Path - The absolute path of the file that was written. + The absolute path of the file that was written. An ec reaction + without enzymes (a zero row in ``model.ec.rxn_enz_mat``) is written + without an ``enzymes`` entry. Raises ------ diff --git a/tests/test_load_ec_model.py b/tests/test_load_ec_model.py index 01febe0..cad9903 100644 --- a/tests/test_load_ec_model.py +++ b/tests/test_load_ec_model.py @@ -554,7 +554,7 @@ def read_with_warning(path): # Input that MATLAB's readYAMLmodel reads and a YAML parser rejects # --------------------------------------------------------------------------- # -# The ec reactions as GECKO writes them: R1_EXP_1 has no enzymes (`enzymes: !!omap` +# The ec reactions as earlier GECKO versions wrote them: R1_EXP_1 has no enzymes (`enzymes: !!omap` # with no entries, as for a `standard` row), R3 is the last row before the next section. # The notes of R1_EXP_1 and R3 hold quotes that RAVEN 2 did not escape. _GECKO_FORM = """\ diff --git a/tests/test_save_ec_model.py b/tests/test_save_ec_model.py index 9c89601..92c6d84 100644 --- a/tests/test_save_ec_model.py +++ b/tests/test_save_ec_model.py @@ -1,4 +1,6 @@ """Tests for save_ec_model.""" +import os +import re from datetime import date from pathlib import Path @@ -413,3 +415,73 @@ def test_round_trip_with_gecko_light_true(tmp_path): path = save_ec_model(model) loaded = load_ec_model(path) assert loaded.ec.gecko_light is True + + +# --------------------------------------------------------------------------- # +# ec reaction without enzymes +# --------------------------------------------------------------------------- # + +def _build_model_with_enzymeless_rxn(adapter: ModelAdapter) -> EcModel: + """`_build_full_ec_model` with R2 as an ec reaction without enzymes.""" + model = _build_full_ec_model(adapter) + mat = model.ec.rxn_enz_mat.tolil() + mat[1, 1] = 0.0 + model.ec.rxn_enz_mat = mat.tocsr() + model.ec.rxn_enz_mat.eliminate_zeros() + model.ec.source[1] = "standard" + return model + + +def test_rxn_without_enzymes_has_no_enzymes_entry(tmp_path): + adapter = _adapter(tmp_path) + path = save_ec_model(_build_model_with_enzymeless_rxn(adapter)) + text = path.read_text(encoding="utf-8") + ec_rxns = text.split("- ec-rxns:")[1].split("- ec-enzymes:")[0] + assert ec_rxns.count("enzymes") == 1 + r2 = re.split(r'- id: "?R2"?\n', ec_rxns)[1] + assert "enzymes" not in r2 + assert "standard" in r2 + assert "!!omap []" not in text + doc = _read_yaml(path) + assert "enzymes" not in doc["ec-rxns"][1] + + +def test_round_trip_rxn_without_enzymes(tmp_path): + adapter = _adapter(tmp_path) + model = _build_model_with_enzymeless_rxn(adapter) + loaded = load_ec_model(save_ec_model(model)) + assert loaded.ec.rxns == ["R1", "R2"] + assert loaded.ec.source == ["brenda", "standard"] + np.testing.assert_array_equal(loaded.ec.kcat, model.ec.kcat) + np.testing.assert_array_equal( + loaded.ec.rxn_enz_mat.toarray(), [[1.0, 0.0], [0.0, 0.0]], + ) + assert loaded.ec.enzymes == ["E1", "E2"] + + +_LIGHT_ECMODEL = Path(os.environ.get("LIGHT_ECMODEL_YML", "/nonexistent")) + + +@pytest.mark.skipif( + not _LIGHT_ECMODEL.is_file(), + reason="set LIGHT_ECMODEL_YML to GECKO's tutorials/light_ecModel/models/ecModel.yml", +) +def test_light_tutorial_model_loses_only_the_empty_enzymes_entries(tmp_path): + """The light tutorial ecModel, written by GECKO 3, has `enzymes: !!omap` with + no entries for its reactions without enzymes. Loaded and saved again, the + file has no `enzymes` entry for them and the ec data is unchanged.""" + loaded = load_ec_model(_LIGHT_ECMODEL) + n_zero = int((loaded.ec.rxn_enz_mat.getnnz(axis=1) == 0).sum()) + assert n_zero == 52 + path = save_ec_model(loaded, tmp_path / "out.yml", adapter=None) + saved = path.read_text(encoding="utf-8") + saved_rxns = saved.split("- ec-rxns:")[1].split("- ec-enzymes:")[0] + assert "!!omap []" not in saved_rxns + assert saved_rxns.count("- id:") == len(loaded.ec.rxns) + assert saved_rxns.count("enzymes:") == len(loaded.ec.rxns) - n_zero + again = load_ec_model(path) + assert again.ec.rxns == loaded.ec.rxns + np.testing.assert_array_equal(again.ec.kcat, loaded.ec.kcat) + assert again.ec.source == loaded.ec.source + assert again.ec.enzymes == loaded.ec.enzymes + assert (again.ec.rxn_enz_mat != loaded.ec.rxn_enz_mat).nnz == 0 From 720498f08933a36e8d7bc0970504298d30a5d500 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 1 Nov 2026 13:23:51 +0100 Subject: [PATCH 26/27] ec_fseof orders targets by the score rounded to 9 significant digits --- CHANGELOG.md | 1 + src/geckopy/utilities/ec_fseof.py | 13 +++++++++---- tests/test_ec_fseof.py | 27 +++++++++++++++++++++++++++ 3 files changed, 37 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e6be01..228aa58 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Unreleased +* `ec_fseof` ranks the targets by the score `|slope| * |correlation|` rounded to 9 significant digits, with scores that are equal after rounding in model order, as GECKO's `ecFSEOF` does; scores that differ by solver noise used to decide the order. * `geckopy.kcat_tuning.evotune` exports `best_parsimonious`, `n_changed`, `parsimony_frontier` and `isozyme_tie_map`, the helpers behind the parsimony sweep and the tying of isozyme copies, which GECKO has as `bestParsimonious`, `nChanged`, `parsimonyFrontier` and `isozymeTieMap`. * `read_kcat_source` and `load_kcat_source` skip a row whose fields are all empty quoted fields (`""""`), as they skip a row of tabs, and raise `ValueError` for a field that opens a double quote and does not close it. Before, the first was reported as an invalid row and the second was read as a row with a kcat that is not a number. * `read_dlkcat_output` drops a row whose kcat is infinite (`inf`), as it drops `NA` and `None`, so the kcat list holds finite numbers only. diff --git a/src/geckopy/utilities/ec_fseof.py b/src/geckopy/utilities/ec_fseof.py index 9bebcee..2ba1283 100644 --- a/src/geckopy/utilities/ec_fseof.py +++ b/src/geckopy/utilities/ec_fseof.py @@ -116,7 +116,7 @@ def ec_fseof( rows filtered out of both ``scan`` and ``targets``. ``targets`` has an added boolean column ``transport`` and lists the reaction targets first, then the transport targets, each group in the - order of the score ``|slope| * |correlation|`` (equal scores in + order of the score ``|slope| * |correlation|`` rounded to 9 significant digits (equal scores in model order); a reaction is a transport target when a metabolite name occurs more than once among its metabolites, not counting the ``prot_*`` ones. A caller that takes the top N rows gets the @@ -219,14 +219,19 @@ def _group_transport_targets(result: FSEOFResult, model: "EcModel") -> FSEOFResu """Add the ``transport`` column and list reaction targets before transport targets. As in ``ecFSEOF.m``, a reaction is a transport reaction when a metabolite name repeats - among its metabolites after the ``prot_*`` metabolites are dropped. The order within a - group is the order of ``result.targets``. + among its metabolites after the ``prot_*`` metabolites are dropped. Within a group the + reactions are ordered by their score rounded to 9 significant digits, equal scores in + model order. """ targets = result.targets.copy() targets["transport"] = [ _is_transport(model.reactions.get_by_id(rid)) for rid in targets["reaction"] ] - targets = targets.sort_values("transport", kind="stable", ignore_index=True) + position = {rxn.id: i for i, rxn in enumerate(model.reactions)} + targets["_model_order"] = [position[rid] for rid in targets["reaction"]] + targets["_score"] = [float(f"{x:.9g}") for x in targets["score"]] + targets = targets.sort_values(["transport", "_score", "_model_order"], ascending=[True, False, True], ignore_index=True) + targets = targets.drop(columns=["_model_order", "_score"]) return FSEOFResult(scan=result.scan, enforced=result.enforced, targets=targets) diff --git a/tests/test_ec_fseof.py b/tests/test_ec_fseof.py index 79048fa..6c988fc 100644 --- a/tests/test_ec_fseof.py +++ b/tests/test_ec_fseof.py @@ -415,3 +415,30 @@ def spy(*args, **kwargs): result = ec_fseof(_build_fseof_model(_adapter(tmp_path)), "EX_PROD", n_steps=4) assert seen["equality"] is True assert "transport" in result.targets.columns + + +def test_scores_equal_to_nine_digits_keep_model_order(): + """Scores that differ only below 9 significant digits count as equal and are listed in + model order.""" + import importlib + + import pandas as pd + + module = importlib.import_module("geckopy.utilities.ec_fseof") + model = cobra.Model("m") + met = cobra.Metabolite("a_c", name="a", compartment="c") + reactions = [] + for rid in ("R1", "R2", "R3"): + rxn = cobra.Reaction(rid) + rxn.add_metabolites({met: -1.0}) + reactions.append(rxn) + model.add_reactions(reactions) + result = FSEOFResult( + scan=pd.DataFrame(), + enforced=[1.0], + targets=pd.DataFrame({"reaction": ["R3", "R2", "R1"], "score": [2.0 + 1e-12, 2.0, 2.5]}), + ) + + grouped = module._group_transport_targets(result, model).targets + + assert list(grouped["reaction"]) == ["R1", "R2", "R3"] From e6c7a8a578d641cd5eaa8de9a7cd5b85345b8a22 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sun, 1 Nov 2026 14:13:33 +0100 Subject: [PATCH 27/27] load_ec_model: comment on the criterion of the reader --- src/geckopy/utilities/load_ec_model.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/geckopy/utilities/load_ec_model.py b/src/geckopy/utilities/load_ec_model.py index 4700a9f..7d5a918 100644 --- a/src/geckopy/utilities/load_ec_model.py +++ b/src/geckopy/utilities/load_ec_model.py @@ -147,9 +147,9 @@ def count_turned(message, category, filename, lineno, file=None, line=None): ec_model.ec = ec if on_fast_solver: switch_to_scipy(ec_model) - # raven-toolbox turns the reactions with a negative lower bound; GECKO - # turns those whose stoichiometry has the GECKO 3 sign, which also - # covers a usage reaction fixed to zero. + # raven-toolbox turns the reactions whose stoichiometry has the GECKO 3 + # sign; the same rule is applied here to a model that was not read from + # YAML. turned += len(upgrade_ec_model(ec_model)) if turned: warnings.warn(