Merge pull request #688 from carrascomj/fix/stoichiometry

fix: detect minimal uncoservable sets

HISTORY.rst

@@ -5,6 +5,7 @@ Next Release
 ------------
 * Fix an issue where experimental growth was incorrectly not reported.
 * Allow the user to set a threshold value for growth in experimental data.
+* Fix and enable the consistency test for minimal unconservable sets.
 
 0.10.2 (2020-03-24)
 -------------------

src/memote/suite/tests/test_consistency.py

@@ -27,7 +27,7 @@
 import memote.support.consistency_helpers as con_helpers
 
 
-@annotate(title="Stoichiometric Consistency", format_type="count")
+@annotate(title="Stoichiometric Consistency", format_type="percent")
 def test_stoichiometric_consistency(model):
     """
     Expect that the stoichiometry is consistent.
@@ -48,19 +48,34 @@ def test_stoichiometric_consistency(model):
     doi: 10.1093/bioinformatics/btn425
     Should the model be inconsistent, then the list of unconserved metabolites
     is computed using the algorithm described in section 3.2 of the same
-    publication.
+    publication. In addition, the list of min unconservable sets is computed
+    using the algorithm described in section 3.3.
 
     """
     ann = test_stoichiometric_consistency.annotation
     is_consistent = consistency.check_stoichiometric_consistency(
         model)
-    ann["data"] = [] if is_consistent else get_ids(
-        consistency.find_unconserved_metabolites(model))
-    ann["metric"] = len(ann["data"]) / len(model.metabolites)
+    ann["data"] = {
+        "unconserved_metabolites":  [] if is_consistent else get_ids(
+            consistency.find_unconserved_metabolites(model)),
+        "minimal_unconservable_sets": [] if is_consistent else [
+            get_ids(mets)
+            for mets in
+            consistency.find_inconsistent_min_stoichiometry(model)],
+    }
+    ann["metric"] = len(ann["data"]["unconserved_metabolites"]) / len(
+        model.metabolites
+    )
     ann["message"] = wrapper.fill(
         """This model contains {} ({:.2%}) unconserved
-        metabolites: {}""".format(
-            len(ann["data"]), ann["metric"], truncate(ann["data"])))
+        metabolites: {}; and {} minimal unconservable sets: {}""".format(
+            len(ann["data"]["unconserved_metabolites"]),
+            ann["metric"],
+            truncate(ann["data"]["unconserved_metabolites"]),
+            len(ann["data"]["minimal_unconservable_sets"]),
+            truncate(ann["data"]["minimal_unconservable_sets"]),
+        )
+    )
     assert is_consistent, ann["message"]
 
 

src/memote/support/consistency.py

@@ -180,8 +180,6 @@ def find_unconserved_metabolites(model):
             " Solver status is '{}' (only optimal expected).".format(status))
 
 
-# FIXME: The results of this function are currently inconsistent with
-# published results.
 def find_inconsistent_min_stoichiometry(model, atol=1e-13):
     """
     Detect inconsistent minimal net stoichiometries.
@@ -219,23 +217,28 @@ def find_inconsistent_min_stoichiometry(model, atol=1e-13):
     metabolites = sorted(internal_mets, key=get_id)
     stoich, met_index, rxn_index = con_helpers.stoichiometry_matrix(
         metabolites, reactions)
-    left_ns = con_helpers.nullspace(stoich.T)
-    # deal with numerical instabilities
-    left_ns[np.abs(left_ns) < atol] = 0.0
-    LOGGER.info("nullspace has dimension %d", left_ns.shape[1])
+    left_ns = con_helpers.nullspace(stoich.T, atol)
+    if left_ns.size != 0:
+        (problem, indicators) = con_helpers.create_milp_problem(
+            left_ns, metabolites, Model, Variable, Constraint, Objective)
+        # We clone the existing configuration in order to apply non-default
+        # settings, for example, for solver verbosity or timeout.
+        problem.configuration = model.problem.Configuration.clone(
+            config=model.solver.configuration,
+            problem=problem
+        )
+        LOGGER.info(
+            "nullspace has dimension %d", left_ns.shape[1]
+            if len(left_ns.shape) > 1 else 0)
+        LOGGER.debug(str(problem))
+    else:
+        LOGGER.info("Left nullspace is empty!")
+        return {(met,) for met in unconserved_mets}
     inc_minimal = set()
-    (problem, indicators) = con_helpers.create_milp_problem(
-        left_ns, metabolites, Model, Variable, Constraint, Objective)
-    # We clone the existing configuration in order to apply non-default
-    # settings, for example, for solver verbosity or timeout.
-    problem.configuration = model.problem.Configuration.clone(
-        config=model.solver.configuration,
-        problem=problem
-    )
-    LOGGER.debug(str(problem))
     cuts = list()
     for met in unconserved_mets:
-        row = met_index[met]
+        # always add the met as an uncoserved set if there is no left nullspace
+        row = met_index[met] if left_ns.size != 0 else None
         if (left_ns[row] == 0.0).all():
             LOGGER.debug("%s: singleton minimal unconservable set", met.id)
             # singleton set!

src/memote/support/consistency_helpers.py

@@ -24,6 +24,7 @@
 
 import numpy as np
 import sympy
+
 from numpy.linalg import svd
 from six import iteritems, itervalues
 from builtins import zip, dict
@@ -115,15 +116,13 @@ def rank(matrix, atol=1e-13, rtol=0):
     rtol : float
         The relative tolerance for a zero singular value.  Singular values less
         than the relative tolerance times the largest singular value are
-        considered to be zero.
+        considered to be zero
 
     Notes
     -----
     If both `atol` and `rtol` are positive, the combined tolerance is the
     maximum of the two; that is::
-
         tol = max(atol, rtol * smax)
-
     Singular values smaller than ``tol`` are considered to be zero.
 
     Returns
@@ -168,9 +167,7 @@ def nullspace(matrix, atol=1e-13, rtol=0.0):  # noqa: D402
     -----
     If both `atol` and `rtol` are positive, the combined tolerance is the
     maximum of the two; that is::
-
         tol = max(atol, rtol * smax)
-
     Singular values smaller than ``tol`` are considered to be zero.
 
     Returns
@@ -274,9 +271,11 @@ def create_milp_problem(kernel, metabolites, Model, Variable, Constraint,
         k_var = Variable("k_{}".format(met.id), type="binary")
         k_vars.append(k_var)
         ns_problem.add([y_var, k_var])
-        # This constraint is equivalent to 0 <= y[i] <= k[i].
+        # These following constraints are equivalent to 0 <= y[i] <= k[i].
         ns_problem.add(Constraint(
             y_var - k_var, ub=0, name="switch_{}".format(met.id)))
+        ns_problem.add(Constraint(
+            y_var, lb=0, name="switch2_{}".format(met.id)))
     ns_problem.update()
     # add nullspace constraints
     for (j, column) in enumerate(kernel.T):

src/memote/support/matrix.py

@@ -58,8 +58,8 @@ def number_independent_conservation_relations(model):
     s_matrix, _, _ = con_helpers.stoichiometry_matrix(
         model.metabolites, model.reactions
     )
-    ln_matrix = con_helpers.nullspace(s_matrix.T)
-    return ln_matrix.shape[1]
+    left_ns = con_helpers.nullspace(s_matrix.T)
+    return left_ns.shape[1] if len(left_ns) > 1 else 0
 
 
 def matrix_rank(model):

tests/test_for_support/test_for_consistency.py

@@ -582,10 +582,9 @@ def test_find_unconserved_metabolites(model, inconsistent):
     assert set([met.id for met in unconserved_mets]) == set(inconsistent)
 
 
-@pytest.mark.xfail(reason="Bug in current implementation.")
 @pytest.mark.parametrize("model, inconsistent", [
     ("textbook", []),
-    ("figure_1", [("A'",), ("B'",), ("C'",)]),
+    ("figure_1", [("A'", "B'", "C'",)]),
     ("equation_8", [("A",), ("B",), ("C",)]),
     ("figure_2", [("X",)]),
 ], indirect=["model"])