diff --git a/CHANGELOG.md b/CHANGELOG.md index 448ee4c..6935470 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ ### Breaking changes * `get_model_from_homology` with `complex_policy="drop"` removes each complex that has a subunit without an ortholog, and removes the reaction only when no complex remains (was: the reaction was removed as soon as one subunit had no ortholog). For the rule `(A and B) or (C and D)` without an ortholog of `B`, the draft has `C and D`. +* `manipulation.contract_model` is the new name of `remove_duplicate_reactions`, after RAVEN's `contractModel`, and the old name is removed. `contract_model` takes `dist_reverse` (default `True`), the `distReverse` option of `contractModel`: with `False`, reactions with the same metabolites and a different reversibility are duplicates. ### Data diff --git a/src/raven_toolbox/manipulation/__init__.py b/src/raven_toolbox/manipulation/__init__.py index 69c0919..a51ade2 100644 --- a/src/raven_toolbox/manipulation/__init__.py +++ b/src/raven_toolbox/manipulation/__init__.py @@ -15,10 +15,10 @@ from .replace import replace_metabolite from .simplify import ( constrain_reversible_reactions, + contract_model, find_duplicate_reactions, group_linear_reactions, remove_dead_end_reactions, - remove_duplicate_reactions, remove_no_flux_reactions, remove_zero_interval_reactions, simplify_model, @@ -34,6 +34,7 @@ "change_reaction_equations", "close_model", "constrain_reversible_reactions", + "contract_model", "convert_to_irreversible", "copy_to_compartment", "expand_model", @@ -44,7 +45,6 @@ "merge_models", "remove_bad_reactions", "remove_dead_end_reactions", - "remove_duplicate_reactions", "remove_genes", "remove_metabolites", "remove_no_flux_reactions", diff --git a/src/raven_toolbox/manipulation/replace.py b/src/raven_toolbox/manipulation/replace.py index 49a4178..f5bfff6 100644 --- a/src/raven_toolbox/manipulation/replace.py +++ b/src/raven_toolbox/manipulation/replace.py @@ -9,7 +9,7 @@ import cobra -from raven_toolbox.manipulation.simplify import remove_duplicate_reactions +from raven_toolbox.manipulation.simplify import contract_model __all__ = ["replace_metabolite"] @@ -109,7 +109,7 @@ def replace_metabolite( rxn.add_metabolites({old: -coef, rep: coef}) model.remove_metabolites([old]) - return remove_duplicate_reactions(model, restrict_to_metabolites=[rep]) + return contract_model(model, restrict_to_metabolites=[rep]) reps = [m for m in model.metabolites if m.name == replacement] if not reps: @@ -142,7 +142,7 @@ def replace_metabolite( survivor, *duplicates = group # Track the metabolite that actually ends up holding this group's # stoichiometry (not necessarily r itself -- the first-encountered - # member of the group survives, same as remove_duplicate_reactions' + # member of the group survives, same as contract_model's # own survivor rule, and r can be the one that gets merged away). # Reactions touching it are what may have become duplicates. survivors.append(survivor) @@ -155,4 +155,4 @@ def replace_metabolite( if to_delete: model.remove_metabolites(to_delete) - return remove_duplicate_reactions(model, restrict_to_metabolites=survivors) + return contract_model(model, restrict_to_metabolites=survivors) diff --git a/src/raven_toolbox/manipulation/simplify.py b/src/raven_toolbox/manipulation/simplify.py index 0dd2961..a585a40 100644 --- a/src/raven_toolbox/manipulation/simplify.py +++ b/src/raven_toolbox/manipulation/simplify.py @@ -12,8 +12,9 @@ model (exchange reactions are opened by default), found by FVA (minimum and maximum flux both zero); RAVEN ``deleteMinMax``. It also removes blocked reactions that are not topological dead ends. -* ``remove_duplicate_reactions``: removes all but one of each set of reactions with identical - stoichiometry and reversibility (bounds merged to the widest); RAVEN ``deleteDuplicates`` +* ``contract_model``: removes all but one of each set of reactions with identical + stoichiometry and reversibility (bounds merged to the widest); RAVEN ``contractModel`` and + the ``deleteDuplicates`` mode of ``simplifyModel`` (detection only: ``find_duplicate_reactions``). * ``constrain_reversible_reactions``: does not remove reactions; it tightens bounds and makes reversible reactions that can only carry flux in one direction irreversible (by FVA); RAVEN @@ -190,12 +191,12 @@ def remove_no_flux_reactions( return blocked -def _signature(rxn, met_index: dict[str, int]): - """Stoichiometry (reversible reactions written in the ``sort_model`` direction) and - reversibility, as RAVEN ``contractModel`` compares reactions.""" +def _signature(rxn, met_index: dict[str, int], dist_reverse: bool = True): + """Stoichiometry (reversible reactions written in the ``sort_model`` direction) and, with + ``dist_reverse``, reversibility, as RAVEN ``contractModel`` compares reactions.""" sign = reversible_orientation(rxn, met_index) mets = frozenset((m.id, sign * c) for m, c in rxn.metabolites.items()) - return (mets, rxn.reversibility) + return (mets, rxn.reversibility) if dist_reverse else (mets,) def _stoich_signature(rxn, *, ignore_direction: bool) -> frozenset: @@ -220,7 +221,7 @@ def find_duplicate_reactions( ) -> list[list[cobra.Reaction]]: """Return groups of reactions that share identical stoichiometry. - This is the detection-only counterpart to :func:`remove_duplicate_reactions`. + This is the detection-only counterpart to :func:`contract_model`. Only the stoichiometry is compared; bounds, objective coefficients, GPRs and annotations are ignored, as in yeast-GEM ``findDuplicatedRxns``. The groups are reactions that could be merged. @@ -247,15 +248,16 @@ def find_duplicate_reactions( return [g for g in groups.values() if len(g) >= 2] -def remove_duplicate_reactions( +def contract_model( model: cobra.Model, *, + dist_reverse: bool = True, restrict_to_metabolites: Iterable[cobra.Metabolite] | None = None, ) -> list[str]: - """Remove all but one of each set of duplicate reactions. + """Contract a model by merging identical reactions (RAVEN ``contractModel``). - Reactions are duplicates when they have identical stoichiometry and the same - reversibility, as in RAVEN ``contractModel``: reversible reactions are compared in the + Reactions are duplicates when they have identical stoichiometry and, with ``dist_reverse``, + the same reversibility, as in RAVEN ``contractModel``: reversible reactions are compared in the direction ``sort_model(sort_reversible=True)`` gives them, so ``A <=> B`` and ``B <=> A`` are duplicates. Bounds and objective coefficients may differ; the survivor takes the lowest lower bound, the highest upper bound and the highest objective coefficient of its @@ -269,13 +271,18 @@ def remove_duplicate_reactions( in ``contractModel``; an isozyme recorded on another duplicate is therefore kept. If every reaction in the set has an ``_EXP_`` suffix (the naming of :func:`~raven_toolbox.manipulation.expand.expand_model`), the suffix is removed - from the id of the survivor, so that ``expand_model`` followed by ``remove_duplicate_reactions`` + from the id of the survivor, so that ``expand_model`` followed by ``contract_model`` restores the original reaction id. Parameters ---------- model : cobra.Model Model modified in place. + dist_reverse : bool + ``True`` (default) treats reactions with the same metabolites and a different + reversibility as different reactions. ``False`` compares the stoichiometry only, so an + irreversible and a reversible reaction with the same metabolites are duplicates, as with + ``distReverse=false`` in ``contractModel``. restrict_to_metabolites: If given, only reactions that contain at least one of these metabolites are checked for duplicates; every other reaction is left unchanged even @@ -289,7 +296,8 @@ def remove_duplicate_reactions( met_index = {m.id: i for i, m in enumerate(model.metabolites)} groups: dict = {} for i, rxn in enumerate(model.reactions): - key = _signature(rxn, met_index) if touching is None or rxn in touching else (object(), i) + key = (_signature(rxn, met_index, dist_reverse) + if touching is None or rxn in touching else (object(), i)) groups.setdefault(key, []).append(rxn) removed: list[str] = [] @@ -552,7 +560,7 @@ def simplify_model( delete_no_flux : bool Run :func:`remove_no_flux_reactions`. The default is ``False``. delete_duplicates : bool - Run :func:`remove_duplicate_reactions` on the whole model. The default is ``False``. + Run :func:`contract_model` on the whole model. The default is ``False``. group_linear : bool Run :func:`group_linear_reactions`, which clears all gene-reaction rules and removes all genes. The default is ``False``. @@ -567,7 +575,7 @@ def simplify_model( Passed to the no-flux mode. The default is ``True``. """ if delete_duplicates: - remove_duplicate_reactions(model) + contract_model(model) if delete_zero_interval: remove_zero_interval_reactions(model, reserved=reserved) if delete_dead_end: diff --git a/tests/parity/test_determinism.py b/tests/parity/test_determinism.py index ab6b13c..135ce54 100644 --- a/tests/parity/test_determinism.py +++ b/tests/parity/test_determinism.py @@ -21,7 +21,7 @@ import pytest from raven_toolbox.gapfilling import connect_blocked_reactions -from raven_toolbox.manipulation import remove_duplicate_reactions +from raven_toolbox.manipulation import contract_model pytestmark = pytest.mark.parity @@ -107,7 +107,7 @@ def test_duplicate_removal_keeps_the_same_reaction_every_run(): {m: c for m, c in model.reactions.R1.metabolites.items()} ) model.add_reactions([clone]) - remove_duplicate_reactions(model) + contract_model(model) kept.append(tuple(sorted(r.id for r in model.reactions))) assert all("R1_dup1" not in run or "R1_dup2" not in run for run in kept), ( @@ -115,7 +115,7 @@ def test_duplicate_removal_keeps_the_same_reaction_every_run(): "selection does -- fix the fixture, not the assertion" ) assert len(set(kept)) == 1, ( - f"remove_duplicate_reactions kept different reactions across {RUNS} " + f"contract_model kept different reactions across {RUNS} " f"runs: {set(kept)}" ) diff --git a/tests/test_manipulation_simplify.py b/tests/test_manipulation_simplify.py index a7f7c50..35dd7fd 100644 --- a/tests/test_manipulation_simplify.py +++ b/tests/test_manipulation_simplify.py @@ -5,9 +5,9 @@ from raven_toolbox.manipulation import ( add_reactions_from_equations, constrain_reversible_reactions, + contract_model, group_linear_reactions, remove_dead_end_reactions, - remove_duplicate_reactions, remove_no_flux_reactions, remove_zero_interval_reactions, simplify_model, @@ -72,7 +72,7 @@ def test_dead_end_ignores_bounds_on_its_own_matching_raven(): assert {"R_in", "R1", "R_locked", "R_out"} == {r.id for r in m.reactions} -# --- remove_duplicate_reactions ------------------------------------------- +# --- contract_model ------------------------------------------- def _two_rxn_model(eq1, b1, eq2, b2, c2=0.0): m = cobra.Model("t") @@ -90,14 +90,14 @@ def _two_rxn_model(eq1, b1, eq2, b2, c2=0.0): def test_duplicates_removed(): m = _two_rxn_model("a --> b", (0, 1000), "a --> b", (0, 1000)) - assert remove_duplicate_reactions(m) == ["R2"] + assert contract_model(m) == ["R2"] assert [r.id for r in m.reactions] == ["R1"] def test_duplicates_with_different_bounds_merge_to_the_widest(): m = _two_rxn_model("a <=> b", (-5, 10), "a <=> b", (-10, 20), c2=1.0) with pytest.warns(UserWarning) as rec: - removed = remove_duplicate_reactions(m) + removed = contract_model(m) assert removed == ["R2"] assert m.reactions.R1.bounds == (-10, 20) assert m.reactions.R1.objective_coefficient == 1.0 @@ -106,17 +106,24 @@ def test_duplicates_with_different_bounds_merge_to_the_widest(): def test_duplicates_differing_in_reversibility_are_kept(): m = _two_rxn_model("a --> b", (0, 10), "a <=> b", (-10, 10)) - assert remove_duplicate_reactions(m) == [] + assert contract_model(m) == [] + + +def test_dist_reverse_false_merges_reactions_that_differ_in_reversibility(): + m = _two_rxn_model("a --> b", (0, 10), "a <=> b", (-10, 10)) + with pytest.warns(UserWarning): + assert contract_model(m, dist_reverse=False) == ["R2"] + assert [r.id for r in m.reactions] == ["R1"] def test_reversible_duplicates_in_opposite_directions_merge(): m = _two_rxn_model("a <=> b", (-10, 10), "b <=> a", (-10, 10)) - assert remove_duplicate_reactions(m) == ["R2"] + assert contract_model(m) == ["R2"] def test_irreversible_opposite_directions_are_kept(): m = _two_rxn_model("a --> b", (0, 10), "b --> a", (0, 10)) - assert remove_duplicate_reactions(m) == [] + assert contract_model(m) == [] def test_duplicates_ignore_reserved_like_contract_model(): @@ -349,7 +356,7 @@ def test_simplify_model_runs_modes_in_raven_order(monkeypatch): calls = [] for name in ( - "remove_duplicate_reactions", "remove_zero_interval_reactions", + "contract_model", "remove_zero_interval_reactions", "remove_dead_end_reactions", "remove_no_flux_reactions", "constrain_reversible_reactions", "group_linear_reactions", ): @@ -359,7 +366,7 @@ def test_simplify_model_runs_modes_in_raven_order(monkeypatch): delete_duplicates=True, group_linear=True, constrain_reversible=True, ) assert calls == [ - "remove_duplicate_reactions", "remove_zero_interval_reactions", + "contract_model", "remove_zero_interval_reactions", "remove_dead_end_reactions", "remove_no_flux_reactions", "constrain_reversible_reactions", "group_linear_reactions", ]