Merge pull request #210 from jonls/prodcheck

gapcheck: Implement alternative production check

docs/commands.rst

@@ -376,7 +376,25 @@ every compound to also be consumed by another reaction (in other words, for
 the purpose of this analysis there are implicit sinks for every compound in
 the model). This means that even if this command reports that no compounds are
 blocked, it may still not be possible for the model to be viable under the
-steady-state assumption of FBA.
+steady-state assumption of FBA. The option ``--no-implicit-sinks`` can be used
+to perform the gap check without implicit sinks.
+
+The gap check is performed with the medium that is defined in the model. It
+may be useful to run the gap check with every compound in the medium available.
+This can easily be done by specifying the ``--unrestricted-medium`` option
+which removes all limits on the exchange reactions during the check.
+
+There are some additional gap checking methods that can be enabled with the
+``--method`` option. The method ``sinkcheck`` can be used to find compounds
+that cannot be synthesized from scratch. The standard gap check will report
+compounds as produced if they can participate in a reaction, even if the
+compound itself cannot be synthesized from precursors in the medium. To find
+such compounds use the ``sinkcheck``. This check will generally indicate more
+compounds as blocked. Lastly, the method ``gapfind`` can be used. This method
+should produce the same result as the default method but is implemented in an
+alternative way that is specified in [Kumar07]_. This method is *not* used by
+default because it tends to result in difficulties for the solver when used
+with larger models.
 
 GapFill (``gapfill``)
 ---------------------

psamm/commands/gapcheck.py

@@ -13,25 +13,43 @@
 # You should have received a copy of the GNU General Public License
 # along with PSAMM.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Copyright 2014-2016  Jon Lund Steffensen <jon_steffensen@uri.edu>
+# Copyright 2014-2017  Jon Lund Steffensen <jon_steffensen@uri.edu>
+# Copyright 2016  Keith Dufault-Thompson <keitht547@my.uri.edu>
 
 from __future__ import unicode_literals
 
 import logging
+from six import iteritems
 
 from ..command import Command, MetabolicMixin, SolverCommandMixin
 from ..gapfill import gapfind, GapFillError
+from ..lpsolver import lp
 
 logger = logging.getLogger(__name__)
 
 
 class GapCheckCommand(MetabolicMixin, SolverCommandMixin, Command):
     """Check for compound production gaps in model."""
+
     @classmethod
     def init_parser(cls, parser):
+        parser.add_argument(
+            '--method', choices=['gapfind', 'prodcheck', 'sinkcheck'],
+            default='prodcheck',
+            help='Method to use for gap checking (default: prodcheck)')
+        parser.add_argument(
+            '--no-implicit-sinks', action='store_true',
+            help='Do not include implicit sinks when gap-filling')
+        parser.add_argument(
+            '--unrestricted-medium', action='store_true',
+            help='Remove all limits on exchange reactions while gap checking')
+        parser.add_argument(
+            '--exclude-extracellular', action='store_true',
+            help='Exclude extracellular compounds from the result')
         parser.add_argument(
             '--epsilon', type=float, default=1e-5,
             help='Threshold for compound production')
+
         super(GapCheckCommand, cls).init_parser(parser)
 
     def run(self):
@@ -43,28 +61,161 @@ def run(self):
             elif compound.id not in compound_name:
                 compound_name[compound.id] = compound.id
 
-        solver = self._get_solver(integer=True)
         extracellular_comp = self._model.extracellular_compartment
         epsilon = self._args.epsilon
         v_max = float(self._model.default_flux_limit)
 
-        # Run GapFind on model
+        implicit_sinks = not self._args.no_implicit_sinks
+
+        if self._args.unrestricted_medium:
+            # Allow all exchange reactions with no flux limits
+            for reaction in self._mm.reactions:
+                if self._mm.is_exchange(reaction):
+                    del self._mm.limits[reaction].bounds
+
         logger.info('Searching for blocked compounds...')
-        result = gapfind(self._mm, solver=solver, epsilon=epsilon, v_max=v_max)
-
-        try:
-            blocked = set(compound for compound in result
-                          if compound.compartment != extracellular_comp)
-        except GapFillError as e:
-            self._log_epsilon_and_fail(epsilon, e)
-
-        if len(blocked) > 0:
-            logger.info('Blocked compounds: {}'.format(len(blocked)))
-            for compound in sorted(blocked):
-                name = compound_name.get(compound.name, compound.name)
-                print('{}\t{}'.format(compound, name))
+
+        if self._args.method == 'gapfind':
+            # Run GapFind on model
+            solver = self._get_solver(integer=True)
+            try:
+                blocked = sorted(gapfind(
+                    self._mm, solver=solver, epsilon=epsilon, v_max=v_max,
+                    implicit_sinks=implicit_sinks))
+            except GapFillError as e:
+                self._log_epsilon_and_fail(epsilon, e)
+        elif self._args.method == 'prodcheck':
+            # Run production check on model
+            solver = self._get_solver()
+            blocked = self.run_reaction_production_check(
+                self._mm, solver=solver, threshold=self._args.epsilon,
+                implicit_sinks=implicit_sinks)
+        elif self._args.method == 'sinkcheck':
+            # Run sink check on model
+            solver = self._get_solver()
+            blocked = self.run_sink_check(
+                self._mm, solver=solver, threshold=self._args.epsilon,
+                implicit_sinks=implicit_sinks)
         else:
-            logger.info('No blocked compounds found')
+            self.argument_error('Invalid method: {}'.format(self._args.method))
+
+        # Show result
+        count = 0
+        for compound in blocked:
+            if (self._args.exclude_extracellular and
+                    compound.compartment == extracellular_comp):
+                continue
+
+            count += 1
+            name = compound_name.get(compound.name, compound.name)
+            print('{}\t{}'.format(compound, name))
+
+        logger.info('Blocked compounds: {}'.format(count))
+
+    def run_sink_check(self, model, solver, threshold, implicit_sinks=True):
+        """Run sink production check method."""
+        prob = solver.create_problem()
+
+        # Create flux variables
+        v = prob.namespace()
+        for reaction_id in model.reactions:
+            lower, upper = model.limits[reaction_id]
+            v.define([reaction_id], lower=lower, upper=upper)
+
+        # Build mass balance constraints
+        massbalance_lhs = {compound: 0 for compound in model.compounds}
+        for spec, value in iteritems(model.matrix):
+            compound, reaction_id = spec
+            massbalance_lhs[compound] += v(reaction_id) * value
+
+        mass_balance_constrs = {}
+        for compound, lhs in iteritems(massbalance_lhs):
+            if implicit_sinks:
+                # The constraint is merely >0 meaning that we have implicit
+                # sinks for all compounds.
+                prob.add_linear_constraints(lhs >= 0)
+            else:
+                # Save these constraints so we can temporarily remove them
+                # to create a sink.
+                c, = prob.add_linear_constraints(lhs == 0)
+                mass_balance_constrs[compound] = c
+
+        for compound, lhs in sorted(iteritems(massbalance_lhs)):
+            if not implicit_sinks:
+                mass_balance_constrs[compound].delete()
+
+            prob.set_objective(lhs)
+            result = prob.solve(lp.ObjectiveSense.Maximize)
+            if not result:
+                logger.warning('Failed to solve for compound: {}'.format(
+                    compound))
+
+            if result.get_value(lhs) < threshold:
+                yield compound
+
+            if not implicit_sinks:
+                # Restore mass balance constraint.
+                c, = prob.add_linear_constraints(lhs == 0)
+                mass_balance_constrs[compound] = c
+
+    def run_reaction_production_check(self, model, solver, threshold,
+                                      implicit_sinks=True):
+        """Run reaction production check method."""
+        prob = solver.create_problem()
+
+        # Create flux variables
+        v = prob.namespace()
+        for reaction_id in model.reactions:
+            lower, upper = model.limits[reaction_id]
+            v.define([reaction_id], lower=lower, upper=upper)
+
+        # Build mass balance constraints
+        massbalance_lhs = {compound: 0 for compound in model.compounds}
+        for spec, value in iteritems(model.matrix):
+            compound, reaction_id = spec
+            massbalance_lhs[compound] += v(reaction_id) * value
+
+        # Create production variables and apply constraints
+        for compound, lhs in iteritems(massbalance_lhs):
+            if implicit_sinks:
+                # The constraint is merely >0 meaning that we have implicit
+                # sinks for all compounds.
+                prob.add_linear_constraints(lhs >= 0)
+            else:
+                prob.add_linear_constraints(lhs == 0)
+
+        confirmed_production = set()
+        for reaction in model.reactions:
+            if all(c in confirmed_production for c, _ in
+                   model.get_reaction_values(reaction)):
+                continue
+
+            prob.set_objective(v(reaction))
+            for sense in (lp.ObjectiveSense.Maximize,
+                          lp.ObjectiveSense.Minimize):
+                result = prob.solve(sense)
+                if not result:
+                    self.fail(
+                        'Failed to solve for compound, reaction: {}, {}:'
+                        ' {}'.format(compound, reaction, result.status))
+
+                flux = result.get_value(v(reaction))
+                for compound, value in model.get_reaction_values(reaction):
+                    if compound in confirmed_production:
+                        continue
+
+                    production = 0
+                    if sense == lp.ObjectiveSense.Maximize and flux > 0:
+                        production = float(value) * flux
+                    elif sense == lp.ObjectiveSense.Minimize and flux < 0:
+                        production = float(value) * flux
+
+                    if production >= threshold:
+                        confirmed_production.add(compound)
+
+        for compound in sorted(model.compounds):
+            if compound not in confirmed_production:
+                yield compound
 
     def _log_epsilon_and_fail(self, epsilon, exc):
         msg = ('Finding blocked compounds failed with epsilon set to {}. Try'

psamm/tests/test_command.py

@@ -371,9 +371,38 @@ def test_run_fva_with_tfba(self):
         self.run_solver_command(
             FluxVariabilityCommand, ['--loop-removal=tfba'], {'integer': True})
 
-    def test_run_gapcheck(self):
+    def test_run_gapcheck_prodcheck(self):
         self.run_solver_command(
-            GapCheckCommand, requirements={'integer': True})
+            GapCheckCommand, ['--method=prodcheck'])
+
+    def test_run_gapcheck_prodcheck_without_implicit_sinks(self):
+        self.run_solver_command(
+            GapCheckCommand, ['--method=prodcheck', '--no-implicit-sinks'])
+
+    def test_run_gapcheck_prodcheck_without_extracellular(self):
+        self.run_solver_command(
+            GapCheckCommand, ['--method=prodcheck', '--exclude-extracellular'])
+
+    def test_run_gapcheck_prodcheck_with_unrestricted_medium(self):
+        self.run_solver_command(
+            GapCheckCommand, ['--method=prodcheck', '--unrestricted-medium'])
+
+    def test_run_gapcheck_sinkcheck(self):
+        self.run_solver_command(
+            GapCheckCommand, ['--method=sinkcheck'])
+
+    def test_run_gapcheck_sinkcheck_without_implicit_sinks(self):
+        self.run_solver_command(
+            GapCheckCommand, ['--method=sinkcheck', '--no-implicit-sinks'])
+
+    def test_run_gapcheck_gapfind(self):
+        self.run_solver_command(
+            GapCheckCommand, ['--method=gapfind'], {'integer': True})
+
+    def test_run_gapcheck_gapfind_without_implicit_sinks(self):
+        self.run_solver_command(
+            GapCheckCommand, ['--method=gapfind', '--no-implicit-sinks'],
+            {'integer': True})
 
     def test_run_gapfill(self):
         self.run_solver_command(