Rename remove_duplicate_reactions to contract_model

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
 

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",

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)

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_<digits>`` 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:

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,15 +107,15 @@ 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), (
         "no duplicates were removed, so this test would pass whatever the "
         "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)}"
     )
 

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",
     ]