diff --git a/CHANGELOG.md b/CHANGELOG.md index abfd2b9..39e0c31 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ Milestones in the raven-toolbox port. For function-level status see ## Unreleased +* `read_yaml_model` turns a `usage_prot_*` reaction forward when its `prot_pool` coefficient is positive and `prot_pool_exchange` when it is negative, as GECKO's `upgradeEcModel` does, and changes the sign of the objective coefficient of a reaction it turns (was: a negative lower bound decided, and the coefficient stayed, which reversed the objective of the model). * `write_yaml_model` writes an ec reaction without enzymes with no `enzymes` entry (was: `enzymes: !!omap []`). `read_yaml_model` reads a missing `enzymes` entry, a content-less `enzymes: !!omap` and `enzymes: !!omap []` as an empty enzyme map, so files written by earlier versions of GECKO, RAVEN and raven-toolbox still load. * `init.fill_tasks` with `resolve_ties=True` on Gurobi restores `MIPGap` and `MIPGapAbs` of `reference_model` after each gap-fill (was: 0.0 and 0.4 stayed set). * `tasks.find_task_essential_reactions` writes its `cache_path` as JSON (was: pickle) and ignores, with a warning, an existing file that is not a JSON cache of the current version, such as a pickle. `checkpoint_every` below 1 raises `ValueError` (was: a write after every task). diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index f2d9312..0557a7b 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -34,8 +34,8 @@ * very old RAVEN files written as a bare ``-`` sequence of single-key mappings rather than one big mapping; * MATLAB GECKO ecModels whose ``usage_prot_*`` and ``prot_pool_exchange`` - reactions use the older reverse-direction convention (negative lower bound, - swapped stoichiometry signs); these are flipped on load so consumers always + reactions use the older reverse-direction convention (swapped + stoichiometry signs, negative flux); these are flipped on load so consumers always see the forward convention. """ from __future__ import annotations @@ -597,19 +597,27 @@ def _flip_legacy_prot_direction(model: cobra.Model) -> None: The current convention (in both geckopy and recent MATLAB GECKO) treats them as ordinary forward reactions with positive flux. When a loaded model still uses the older convention we flip the affected - reactions in place so consumers never have to handle two shapes. + reactions in place so consumers never have to handle two shapes. The + objective coefficient of a flipped reaction changes sign with its flux. - The signature we look for is any ``usage_prot_*`` or - ``prot_pool_exchange`` reaction whose lower bound is negative. + The signature we look for is the sign of the ``prot_pool`` coefficient: + positive in a ``usage_prot_*`` reaction, negative in + ``prot_pool_exchange``. The bounds do not decide, so a reaction with + the current stoichiometry and a negative lower bound is left as it is. """ + try: + pool = model.metabolites.get_by_id("prot_pool") + except KeyError: + return flipped: list[str] = [] for rxn in model.reactions: - if not ( - rxn.id.startswith("usage_prot_") - or rxn.id == "prot_pool_exchange" - ): + if rxn.id == "prot_pool_exchange": + legacy = rxn.metabolites.get(pool, 0.0) < 0 + elif rxn.id.startswith("usage_prot_"): + legacy = rxn.metabolites.get(pool, 0.0) > 0 + else: continue - if rxn.lower_bound >= -1e-9: + if not legacy: continue rxn.add_metabolites( {met: -2.0 * coef for met, coef in rxn.metabolites.items()}, @@ -618,6 +626,7 @@ def _flip_legacy_prot_direction(model: cobra.Model) -> None: lb, ub = rxn.lower_bound, rxn.upper_bound rxn.lower_bound = -ub rxn.upper_bound = -lb + rxn.objective_coefficient = -rxn.objective_coefficient flipped.append(rxn.id) if flipped: warnings.warn( diff --git a/tests/test_io_yaml_ec_data.py b/tests/test_io_yaml_ec_data.py index 786c296..ccc4fb7 100644 --- a/tests/test_io_yaml_ec_data.py +++ b/tests/test_io_yaml_ec_data.py @@ -341,6 +341,7 @@ def test_legacy_reverse_direction_prot_flipped(tmp_path): "metabolites": {"prot_P1": -1, "prot_pool": 1}, # swapped signs "lower_bound": -1000.0, "upper_bound": 0.0, + "objective_coefficient": 2.0, } ) with pytest.warns(UserWarning, match="reverse-direction"): @@ -352,6 +353,27 @@ def test_legacy_reverse_direction_prot_flipped(tmp_path): coefs = {m.id: c for m, c in r.metabolites.items()} assert coefs["prot_P1"] == 1.0 assert coefs["prot_pool"] == -1.0 + # The objective coefficient changes sign with the flux. + assert r.objective_coefficient == -2.0 + + +def test_forward_prot_reaction_with_negative_lower_bound_not_flipped(tmp_path): + """The sign of the ``prot_pool`` coefficient decides, not the bounds.""" + doc = _minimal_model_doc() + doc["metabolites"].append({"id": "prot_P1", "name": "prot", "compartment": "c"}) + doc["metabolites"].append({"id": "prot_pool", "name": "pool", "compartment": "c"}) + doc["reactions"].append( + { + "id": "usage_prot_P1", + "metabolites": {"prot_P1": 1, "prot_pool": -1}, + "lower_bound": -1000.0, + "upper_bound": 1000.0, + } + ) + model = read_yaml_model(_write_yaml(doc, tmp_path / "m.yml")) + r = model.reactions.get_by_id("usage_prot_P1") + assert (r.lower_bound, r.upper_bound) == (-1000.0, 1000.0) + assert {m.id: c for m, c in r.metabolites.items()}["prot_pool"] == -1.0 def test_legacy_bare_sequence_root_merged(tmp_path):