Merge pull request #684 from carrascomj/fix-nan-growth

refactor: allow configuring a minimal growth rate to determine viability

HISTORY.rst

@@ -3,6 +3,8 @@ History
 
 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.
 
 0.10.2 (2020-03-24)
 -------------------
@@ -184,14 +186,14 @@ Next Release
 
 0.8.9 (2018-12-11)
 ------------------
-* Compress JSON and SQLite storage of results using gzip by default. JSON 
-  continues to work either compressed or uncompressed. At the moment we 
-  offer no database migration, please contact us if you need help in 
+* Compress JSON and SQLite storage of results using gzip by default. JSON
+  continues to work either compressed or uncompressed. At the moment we
+  offer no database migration, please contact us if you need help in
   migrating a large existing SQLite database rather than just re-computing it.
 
 0.8.8 (2018-12-10)
 ------------------
-* Adjust the reversibility index test to not use name matching and increase 
+* Adjust the reversibility index test to not use name matching and increase
   the threshold slightly. Also adjust the description of the test.
 * Adjust tests to the change in the ``add_boundary`` interface.
 * Identify blocked reactions using the cobrapy built-in function.
@@ -202,16 +204,16 @@ Next Release
 * Add a test that checks if reactions are annotated with reactome identifiers.
 * Add a feature that allows identifying specific metabolites by matching
   annotation information against the metabolite shortlist for a given MNX ID.
-* Change every usage of SBO key to lower case to conform to the identifiers.org 
+* Change every usage of SBO key to lower case to conform to the identifiers.org
   namespace for the Systems Biology Ontology.
-* Remove that metabolite SBO terms are used when identifying transport 
+* Remove that metabolite SBO terms are used when identifying transport
   reactions as this may lead to false positives.
-* Return metabolite IDs when finding duplicate metabolites to avoid 
+* Return metabolite IDs when finding duplicate metabolites to avoid
   serialization errors.
 * Identify transport reactions first by formula then by annotation.
 * For the diff report, run pytest in different processes to avoid accidentally
   overwriting the results of the former with the results of the later runs.
-* In the diff report, fix a typo that allowed the diff button to depart the 
+* In the diff report, fix a typo that allowed the diff button to depart the
   defined colour scheme (blue -> red) to cyan.
 * Fix the snapshot report not showing environment information.
 * Allow ``memote run`` to skip commits where the model was not
@@ -253,7 +255,7 @@ Next Release
   category are grouped logically.
 * Updated the documentation to include a newer flowchart, up-to-date getting
   started and custom test sections.
-* Update code to account for breaking changes in the most recent version of 
+* Update code to account for breaking changes in the most recent version of
   cobrapy (0.13.0) and subsequently unpin cobrapy dependency (set to >=0.13.0).
 
 0.8.1 (2018-06-27)

docs/experimental_data.rst

@@ -49,6 +49,27 @@ The file simply states that files relating to media, essentiality, and growth
 experiments can be found at directories stated under ``path:`` relative to the
 location of the file.
 
+Minimal growth rate
+-------------------
+Both growth and essentiality experiments require a minimal biomass function
+value to determine if the model predicts viability. This value can be set via
+the ``minimal_growth_rate`` option in the ``experiments.yml`` file.
+For instance:
+
+.. code-block:: yaml
+
+    version: "0.1"
+    medium:
+      path: "media/"
+    essentiality:
+      path: "essentiality/"
+    growth:
+      path: "growth/"
+    minimal_growth_rate: 0.05
+
+By default, ``minimal_growth_rate`` is set to 10% of the biomass function value
+under the default constraints in the model.
+
 Media
 =====
 

src/memote/experimental/config.py

@@ -21,6 +21,7 @@
 
 import json
 import logging
+from math import isnan
 from io import open
 from os.path import dirname, join, isabs
 
@@ -120,6 +121,7 @@ def load_essentiality(self, model):
             return
         path = self.get_path(data,
                              join("data", "experimental", "essentiality"))
+        minimal_growth_rate = self.get_minimal_growth_rate(model)
         for exp_id, exp in iteritems(experiments):
             if exp is None:
                 exp = dict()
@@ -129,7 +131,9 @@ def load_essentiality(self, model):
             elif not isabs(filename):
                 filename = join(path, filename)
             experiment = EssentialityExperiment(
-                identifier=exp_id, obj=exp, filename=filename)
+                identifier=exp_id, obj=exp,
+                filename=filename, minimal_growth_rate=minimal_growth_rate
+            )
             if experiment.medium is not None:
                 assert experiment.medium in self.media, \
                     "Experiment '{}' has an undefined medium '{}'.".format(
@@ -149,6 +153,7 @@ def load_growth(self, model):
             return
         path = self.get_path(data,
                              join("data", "experimental", "growth"))
+        minimal_growth_rate = self.get_minimal_growth_rate(model)
         for exp_id, exp in iteritems(experiments):
             if exp is None:
                 exp = dict()
@@ -158,7 +163,9 @@ def load_growth(self, model):
             elif not isabs(filename):
                 filename = join(path, filename)
             growth = GrowthExperiment(
-                identifier=exp_id, obj=exp, filename=filename)
+                identifier=exp_id, obj=exp,
+                filename=filename, minimal_growth_rate=minimal_growth_rate
+            )
             if growth.medium is not None:
                 assert growth.medium in self.media, \
                     "Growth-experiment '{}' has an undefined medium '{}'." \
@@ -176,3 +183,30 @@ def get_path(self, obj, default):
         if not isabs(path):
             path = join(self._base, path)
         return path
+
+    def get_minimal_growth_rate(self, model, threshold=0.1):
+        """Calculate min growth default value or return input value.
+
+        This value is used to determine if a model is capable of growth under
+        certain experimental conditions.
+
+        Parameters
+        ----------
+        model : cobra.Model
+        threshold : float, optional
+            If no input is provided by the user the default value is set to a
+            coefficient `threshold` times the growth under default constraints
+            (default: 0.1).
+
+        """
+        minimal_growth_rate = self.config.get("minimal_growth_rate")
+        if minimal_growth_rate is None:
+            minimal_growth_rate = model.slim_optimize() * threshold
+            if isnan(minimal_growth_rate):
+                LOGGER.error(
+                    "Threshold set to {} due to infeasible "
+                    "solution (NaN produced) with default "
+                    "constraints.".format(model.tolerance)
+                )
+                minimal_growth_rate = model.tolerance
+        return minimal_growth_rate

src/memote/experimental/essentiality.py

@@ -87,10 +87,9 @@ def evaluate(self, model):
             if self.objective is not None:
                 model.objective = self.objective
             model.add_cons_vars(self.constraints)
-            max_val = model.slim_optimize()
             essen = single_gene_deletion(
                 model, gene_list=self.data["gene"], processes=1)
         essen["gene"] = [list(g)[0] for g in essen.index]
-        essen["essential"] = (essen["growth"] < (max_val * 0.1)) \
+        essen["essential"] = (essen["growth"] < self.minimal_growth_rate) \
             | essen["growth"].isna()
         return essen

src/memote/experimental/experiment.py

@@ -31,13 +31,16 @@
 class Experiment(ExperimentalBase):
     """Base class for experimental data."""
 
-    def __init__(self, obj, **kwargs):
+    def __init__(self, obj, minimal_growth_rate, **kwargs):
         """
         Initialize an experiment.
 
         Parameters
         ----------
         obj : dict
+        minimal_growth_rate : float
+            minimum value of biomass function for the model to be considered as
+            growing. Not used by Medium. Default: None
         kwargs
 
         """
@@ -51,6 +54,7 @@ def __init__(self, obj, **kwargs):
             raise NotImplementedError(
                 "This feature is not implemented yet. Please let us know that "
                 "you want it at https://github.com/opencobra/memote/issues.")
+        self.minimal_growth_rate = minimal_growth_rate
 
     def evaluate(self, model):
         """Abstract base method for evaluating experimental data."""

src/memote/experimental/experimental_base.py

@@ -50,7 +50,9 @@ class ExperimentalBase(object):
         "YES"
     }
 
-    def __init__(self, identifier, obj, filename, **kwargs):
+    def __init__(
+        self, identifier, obj, filename, **kwargs
+    ):
         """
         Initialize a medium.
 

src/memote/experimental/growth.py

@@ -64,15 +64,14 @@ def load(self, dtype_conversion=None):
         super(GrowthExperiment, self).load(dtype_conversion=dtype_conversion)
         self.data["growth"] = self.data["growth"].isin(self.TRUTHY)
 
-    def evaluate(self, model, threshold=0.1):
+    def evaluate(self, model):
         """Evaluate in silico growth rates."""
         with model:
             if self.medium is not None:
                 self.medium.apply(model)
             if self.objective is not None:
                 model.objective = self.objective
             model.add_cons_vars(self.constraints)
-            threshold *= model.slim_optimize()
             growth = list()
             for row in self.data.itertuples(index=False):
                 with model:
@@ -81,7 +80,9 @@ def evaluate(self, model, threshold=0.1):
                         exchange.lower_bound = -row.uptake
                     else:
                         exchange.upper_bound = row.uptake
-                    growth.append(model.slim_optimize() >= threshold)
+                    growth.append(
+                        model.slim_optimize() >= self.minimal_growth_rate
+                    )
         return DataFrame({
             "exchange": self.data["exchange"],
             "growth": growth

src/memote/experimental/schemata/configuration.json

@@ -118,6 +118,12 @@
     },
     "growth": {
       "$ref": "#/definitions/experiments"
+    },
+    "minimal_growth_rate": {
+      "type": [
+        "number",
+        "null"
+      ]
     }
   },
   "required": [

tests/test_for_experimental/data/growth.yml

@@ -3,6 +3,7 @@ medium:
   path: "media/"
   definitions:
     m9_minimal:
+minimal_growth_rate: 0.002
 growth:
   path: "growth"
   experiments:

tests/test_for_experimental/test_for_experimental_base.py

@@ -49,7 +49,16 @@ def klass(request):
     {"label": "bird"}
 ])
 def test_init(klass, obj):
-    exp = klass(identifier="that", obj=obj, filename=join(DATA_PATH, "."))
+    kwargs = {
+        "identifier": "that",
+        "obj": obj,
+        "filename": join(DATA_PATH, "."),
+    }
+    if issubclass(klass, Experiment):
+        # minimal_growth_rate is required and only used by classes derived
+        # from Experiment (GrowthExperiment and EssentialityExperiment)
+        kwargs["minimal_growth_rate"] = 0.002
+    exp = klass(**kwargs)
     for key, value in iteritems(obj):
         assert getattr(exp, key) == value