From 138db0d13e0ac5ae7cbe53f6c1eff4fa908b5427 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Fri, 7 Aug 2015 12:08:42 -0400 Subject: [PATCH 01/46] native: Add property to access the root FilePathContext This can be used to obtain information on where the model was loaded from. --- psamm/datasource/native.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/psamm/datasource/native.py b/psamm/datasource/native.py index 0f3a722d..aba2b58a 100644 --- a/psamm/datasource/native.py +++ b/psamm/datasource/native.py @@ -239,6 +239,10 @@ def parse_compounds(self): self._context, self._model['compounds']): yield compound + @property + def context(self): + return self._context + def parse_compound(compound_def, context=None): """Parse a structured compound definition as obtained from a YAML file From 4e298d80896e74571bacf5e7fe3b2b414d41f058 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Fri, 7 Aug 2015 12:10:14 -0400 Subject: [PATCH 02/46] command: Log model name after loading the model --- psamm/command.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/psamm/command.py b/psamm/command.py index 4bd64e92..7e28a59e 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -72,6 +72,11 @@ def __init__(self, model, args): self._model = model self._args = args + name = self._model.get_name() + if name is None: + name = str(self._model.context) + logger.info('Model: {}'.format(name)) + # Create metabolic model database = DictDatabase() for reaction in model.parse_reactions(): From b76488e62424cebf828d06fe4d11e258a0a2381a Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Fri, 7 Aug 2015 12:14:03 -0400 Subject: [PATCH 03/46] util: Add function to describe the current Git repository version --- psamm/util.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/psamm/util.py b/psamm/util.py index 606768af..eff8f491 100644 --- a/psamm/util.py +++ b/psamm/util.py @@ -19,6 +19,7 @@ import re import math +import subprocess from six import iteritems @@ -47,6 +48,25 @@ def flush(self): This is a noop.""" +def git_try_describe(repo_path): + """Try to describe the current commit of a Git repository. + + Return a string containing a string with the commit ID and/or a base tag, + if successful. Otherwise, return None. + """ + try: + p = subprocess.Popen(['git', 'describe', '--always', '--dirty'], + cwd=repo_path, stdout=subprocess.PIPE) + output, _ = p.communicate() + except: + return None + else: + if p.returncode == 0: + return output.strip() + + return None + + def convex_cardinality_relaxed(f, epsilon=1e-5): """Transform L1-norm optimization function into cardinality optimization. From 0a10f32dbb60cd2148a523a0801e952749869134 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Fri, 7 Aug 2015 12:14:25 -0400 Subject: [PATCH 04/46] command: Log Git repository version if available --- psamm/command.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/psamm/command.py b/psamm/command.py index 7e28a59e..45316e47 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -47,6 +47,7 @@ from .datasource import sbml from . import fluxanalysis, fluxcoupling, massconsistency, fastcore from .lpsolver import generic +from . import util from six import add_metaclass, iteritems @@ -77,6 +78,10 @@ def __init__(self, model, args): name = str(self._model.context) logger.info('Model: {}'.format(name)) + version = util.git_try_describe(self._model.context.basepath) + if version is not None: + logger.info('Model Git version: {}'.format(version)) + # Create metabolic model database = DictDatabase() for reaction in model.parse_reactions(): From 16a136a69b8e43c3eb949db98a0fdb85adb289d0 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Sun, 9 Aug 2015 12:54:37 -0400 Subject: [PATCH 05/46] Fix omissions from test coverage report Properly configure tox environment and tests to be excluded from coverage. --- tox.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tox.ini b/tox.ini index b6faa4e7..96a6063c 100644 --- a/tox.ini +++ b/tox.ini @@ -20,7 +20,7 @@ deps = cplex: {env:CPLEX_PYTHON_PACKAGE} qsoptex: python-qsoptex>=0.4 commands = - coverage run -p --branch --omit={envdir},psamm/tests \ + coverage run -p --branch --omit={envdir}/*,env/*,psamm/tests/*,setup.py \ ./setup.py test [testenv:coverage] From 20fcc4c37f1fdb0b577d02acaccf167dd74a3391 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 17 Aug 2015 16:13:03 -0400 Subject: [PATCH 06/46] command: Use CommandError for errors running a command This particular exception is caught in main() and reported through the parser which prints out the usage information. This should be used to report errors that originate from misspecified command line arguments (e.g. specifying a float outside the valid range, specifying a non-existent model reaction, etc.) --- psamm/command.py | 48 +++++++++++++++++++++++++++++++----------------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/psamm/command.py b/psamm/command.py index 45316e47..72536e8d 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -55,6 +55,16 @@ logger = logging.getLogger(__name__) +class CommandError(Exception): + """Error from running a command. + + This should be raised from a ``Command.run()`` if any arguments are + misspecified. When the command is run and the ``CommandError`` is raised, + the caller will exit with an error code and print appropriate usage + information. + """ + + @add_metaclass(abc.ABCMeta) class Command(object): """Represents a command in the interface, operating on a model. @@ -337,11 +347,12 @@ def run(self): else: maximized_reaction = self._model.get_biomass_reaction() if maximized_reaction is None: - raise ValueError('The maximized reaction was not specified') + raise CommandError('The maximized reaction was not specified') if not not self._mm.has_reaction(maximized_reaction): - raise ValueError(('The biomass reaction is not a valid model' + - ' reaction: {}').format(maximized_reaction)) + raise CommandError( + 'The biomass reaction is not a valid model' + ' reaction: {}'.format(maximized_reaction)) logger.info('Flux balance on induced model maximizing {}'.format( maximized_reaction)) @@ -404,10 +415,10 @@ def run(self): else: reaction = self._model.get_biomass_reaction() if reaction is None: - raise ValueError('The biomass reaction was not specified') + raise CommandError('The biomass reaction was not specified') if not self._mm.has_reaction(reaction): - raise ValueError('Specified reaction is not in model: {}'.format( + raise CommandError('Specified reaction is not in model: {}'.format( reaction)) if self._args.no_tfba: @@ -582,7 +593,7 @@ def run(self): max_reaction = self._model.get_biomass_reaction() if max_reaction is None: - raise ValueError('The biomass reaction was not specified') + raise CommandError('The biomass reaction was not specified') fba_fluxes = dict(fluxanalysis.flux_balance( self._mm, max_reaction, tfba=False, solver=solver)) @@ -708,10 +719,10 @@ def run(self): else: reaction = self._model.get_biomass_reaction() if reaction is None: - raise ValueError('The biomass reaction was not specified') + raise CommandError('The biomass reaction was not specified') if not self._mm.has_reaction(reaction): - raise ValueError('Specified reaction is not in model: {}'.format( + raise CommandError('Specified reaction is not in model: {}'.format( reaction)) enable_tfba = not self._args.no_tfba @@ -991,15 +1002,15 @@ def run(self): else: reaction = self._model.get_biomass_reaction() if reaction is None: - raise ValueError('The biomass reaction was not specified') + raise CommandError('The biomass reaction was not specified') if not self._mm.has_reaction(reaction): - raise ValueError('Specified reaction is not in model: {}'.format( + raise CommandError('Specified reaction is not in model: {}'.format( reaction)) threshold = self._args.threshold if threshold < 0.0 or threshold > 1.0: - raise ValueError( + raise CommandError( 'Invalid threshold, must be in [0;1]: {}'.format(threshold)) if not self._args.tfba: @@ -1139,20 +1150,20 @@ def run_tfba(model, reaction): else: reaction = self._model.get_biomass_reaction() if reaction is None: - raise ValueError('The biomass reaction was not specified') + raise CommandError('The biomass reaction was not specified') if not self._mm.has_reaction(reaction): - raise ValueError('Specified reaction is not in model: {}'.format( + raise CommandError('Specified reaction is not in model: {}'.format( reaction)) varying_reaction = self._args.varying if not self._mm.has_reaction(varying_reaction): - raise ValueError('Specified reaction is not in model: {}'.format( + raise CommandError('Specified reaction is not in model: {}'.format( varying_reaction)) steps = self._args.steps if steps <= 0: - raise ValueError('Invalid number of steps: {}\n'.format(steps)) + raise CommandError('Invalid number of steps: {}\n'.format(steps)) # Run FBA on model at different fixed flux values flux_min = self._mm.limits[varying_reaction].lower @@ -1163,7 +1174,7 @@ def run_tfba(model, reaction): flux_max = self._args.maximum if flux_min > flux_max: - raise ValueError('Invalid flux range: {}, {}\n'.format( + raise CommandError('Invalid flux range: {}, {}\n'.format( flux_min, flux_max)) for i in xrange(steps): @@ -1374,7 +1385,10 @@ def main(command_class=None): # Instantiate command with model and run command = args.command(model, args) - command.run() + try: + command.run() + except CommandError as e: + parser.error(str(e)) if __name__ == '__main__': From 4827ef44f930f3b30853ea774c80589a2c915b73 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 18 Aug 2015 15:03:59 -0400 Subject: [PATCH 07/46] util: Add helper type for possible relative parameters --- psamm/tests/test_util.py | 55 ++++++++++++++++++++++++++++++++++ psamm/util.py | 65 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 120 insertions(+) create mode 100644 psamm/tests/test_util.py diff --git a/psamm/tests/test_util.py b/psamm/tests/test_util.py new file mode 100644 index 00000000..f14303e7 --- /dev/null +++ b/psamm/tests/test_util.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2015 Jon Lund Steffensen + +import unittest + +from psamm.util import MaybeRelative + + +class TestMaybeRelative(unittest.TestCase): + def test_init_from_float(self): + arg = MaybeRelative(24.5) + self.assertFalse(arg.relative) + self.assertIsNone(arg.reference) + self.assertAlmostEqual(float(arg), 24.5) + + def test_init_from_float_string(self): + arg = MaybeRelative('-10') + self.assertFalse(arg.relative) + self.assertIsNone(arg.reference) + self.assertAlmostEqual(float(arg), -10.0) + + def test_init_from_percentage(self): + arg = MaybeRelative('110%') + self.assertTrue(arg.relative) + self.assertIsNone(arg.reference) + with self.assertRaises(ValueError): + float(arg) + + def test_resolve_relative(self): + arg = MaybeRelative('40%') + arg.reference = 200.0 + self.assertAlmostEqual(float(arg), 80.0) + + def test_invalid(self): + with self.assertRaises(ValueError): + arg = MaybeRelative('abc') + + +if __name__ == '__main__': + unittest.main() diff --git a/psamm/util.py b/psamm/util.py index eff8f491..2074dc4b 100644 --- a/psamm/util.py +++ b/psamm/util.py @@ -48,6 +48,71 @@ def flush(self): This is a noop.""" +class MaybeRelative(object): + """Helper type for parsing possibly relative parameters. + + >>> arg = MaybeRelative('40%') + >>> arg.reference = 200.0 + >>> float(arg) + 80.0 + + >>> arg = MaybeRelative('24.5') + >>> arg.reference = 150.0 + >>> float(arg) + 24.5 + """ + + def __init__(self, s): + try: + self._value = float(s) + self._relative = False + except ValueError: + self._value = self._parse_percentage(s) + self._relative = True + + self._reference = None + + @classmethod + def _parse_percentage(cls, s): + """Parse string as a percentage (e.g. '42.0%') and return as float.""" + m = re.match('^(.+)%$', s) + if not m: + raise ValueError('Unable to parse as percentage: {}'.format( + repr(s))) + + return float(m.group(1)) / 100.0 + + @property + def relative(self): + """Whether the parsed number was relative.""" + return self._relative + + @property + def reference(self): + """The reference used for converting to absolute value.""" + return self._reference + + @reference.setter + def reference(self, value): + self._reference = value + + def __float__(self): + if self._relative: + if self._reference is None: + raise ValueError('Reference not set!') + return self._reference * self._value + else: + return self._value + + def __repr__(self): + if self._relative: + return '<{}, {:.1%} of {}>'.format( + self.__class__.__name__, self._value, self._reference) + else: + return '<{}, {}>'.format( + self.__class__.__name__, self._value) + + def git_try_describe(repo_path): """Try to describe the current commit of a Git repository. From 458ecd7719349981f419e21bbc61fdd84b8410d4 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 18 Aug 2015 15:05:19 -0400 Subject: [PATCH 08/46] command: Use MaybeRelative to parse randomsparse threshold If the parameter is relative, an initial FBA is performed and used as a reference. --- psamm/command.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/psamm/command.py b/psamm/command.py index 72536e8d..5c61020b 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -976,6 +976,9 @@ class RandomSparseNetworkCommand(SolverCommandMixin, Command): Keep deleting reactions until no more reactions can be deleted. By default this uses standard FBA (not tFBA). Since the internal fluxes are irrelevant the FBA and tFBA are equivalent for this purpose. + + The threshold can be specified as an absolute flux (e.g. '1.23') or a + relative flux of the full model flux (e.g. '40.5%'). """ name = 'randomsparse' @@ -987,10 +990,10 @@ def init_parser(cls, parser): '--tfba', help='Enable thermodynamic constraints on FBA', action='store_true', default=False) parser.add_argument( - '--reaction', help='Reaction to maximize', nargs='?') + '--reaction', help='Reaction to maximize') parser.add_argument( - 'threshold', help='Relative threshold of max reaction flux', - type=float) + 'threshold', help='Threshold of max reaction flux', + type=util.MaybeRelative) parser.add_argument( '--exchange', help='Only analyze the exchange reaction in model', action='store_true') @@ -1008,11 +1011,6 @@ def run(self): raise CommandError('Specified reaction is not in model: {}'.format( reaction)) - threshold = self._args.threshold - if threshold < 0.0 or threshold > 1.0: - raise CommandError( - 'Invalid threshold, must be in [0;1]: {}'.format(threshold)) - if not self._args.tfba: fb_problem = fluxanalysis.FluxBalanceProblem solver = self._get_solver() @@ -1021,8 +1019,13 @@ def run(self): solver = self._get_solver(integer=True) p = fb_problem(self._mm, solver) - p.solve(reaction) - flux_threshold = p.get_flux(reaction) * threshold + + threshold = self._args.threshold + if threshold.relative: + p.solve(reaction) + threshold.reference = p.get_flux(reaction) + + flux_threshold = float(threshold) logger.info('Flux threshold for {} is {}'.format( reaction, flux_threshold)) From cae7e2423dfc0a22335786a32ab2e5b799ec3304 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 18 Aug 2015 15:19:47 -0400 Subject: [PATCH 09/46] docs: Update documentation on randomsparse command --- docs/commands.rst | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/docs/commands.rst b/docs/commands.rst index 6cfb505c..7aedb010 100644 --- a/docs/commands.rst +++ b/docs/commands.rst @@ -153,12 +153,13 @@ specifying the reaction explicitly. .. code-block:: shell - $ psamm-model randomsparse 0.95 + $ psamm-model randomsparse 95% When the given reaction is the biomass reaction, this results in a smaller model which is still producing biomass within the tolerance given by the -threshold. Aggregating the results from multiple random sparse networks allows -classifying reactions as essential, semi-essential or non-essential. +threshold. The tolerance can be specified as a relative value (as above) or as +an absolute flux. Aggregating the results from multiple random sparse networks +allows classifying reactions as essential, semi-essential or non-essential. If the option ``--exchange`` is given, the model will only try to delete exchange reactions. This can be used to provide putative minimal media for @@ -166,8 +167,8 @@ the model. The output of the command is a tab-separated list of reaction IDs and a value indicating whether the reaction was eliminated (``0`` when eliminated, ``1`` -otherwise). If multiply minimal networks are desired, the command can be run -again and it will produce a different random minimal network. +otherwise). If multiple minimal networks are desired, the command can be run +again and it will sample another random minimal network. Flux coupling analysis (``fluxcoupling``) ----------------------------------------- From d30ec7296559d89d3f7645f8a1f75447e0343c4e Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 18 Aug 2015 15:35:53 -0400 Subject: [PATCH 10/46] lpsolver: Fix warning about docstring from Sphinx --- psamm/lpsolver/lp.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psamm/lpsolver/lp.py b/psamm/lpsolver/lp.py index 769cc89b..07db6ae9 100644 --- a/psamm/lpsolver/lp.py +++ b/psamm/lpsolver/lp.py @@ -337,11 +337,11 @@ def set(self, names): @abc.abstractmethod def add_linear_constraints(self, *relations): - """Add constraints to the problem + """Add constraints to the problem. Each constraint is given as a :class:`.Relation`, and the expression in that relation can be a set expression. Returns a sequence of - :class:`.Constraint`s. + :class:`Constraints <.Constraint>`. """ @abc.abstractmethod From 94ed3cddec169f6491134dafae37bdd25ccb7430 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 18 Aug 2015 15:37:14 -0400 Subject: [PATCH 11/46] command: Avoid randomsparse model copying by modifying constraints Previously, when deleting reactions from the model a copy of the model was made and the flux bounds in that copy were set to 0. This resulted in a new LP problem being generated for each iteration. Using the new interface to modifying constraints in the LP, a problem can be created once and the constraints can be modified directly. This is much faster then creating a new problem each time. --- psamm/command.py | 38 +++++++++++++++----------------------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/psamm/command.py b/psamm/command.py index 5c61020b..4ec23110 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -1030,38 +1030,35 @@ def run(self): logger.info('Flux threshold for {} is {}'.format( reaction, flux_threshold)) + essential = {reaction} + deleted = set() if self._args.exchange: - model_test = self._mm.copy() - essential = {reaction} - deleted = set() - exchange = set() + reactions = set() for reaction_id in self._mm.reactions: if self._mm.is_exchange(reaction_id): - exchange.add(reaction_id) - test_set = set(exchange) - essential + reactions.add(reaction_id) else: - model_test = self._mm.copy() - essential = {reaction} - deleted = set() - test_set = set(self._mm.reactions) - essential + reactions = set(self._mm.reactions) + + test_set = reactions - essential while len(test_set) > 0: testing_reaction = random.sample(test_set, 1)[0] test_set.remove(testing_reaction) - saved_bounds = model_test.limits[testing_reaction].bounds - model_test.limits[testing_reaction].bounds = 0, 0 + + flux_var = p.get_flux_var(testing_reaction) + c, = p.prob.add_linear_constraints(flux_var == 0) logger.info('Trying FBA without reaction {}...'.format( testing_reaction)) - p = fb_problem(model_test, solver) try: p.solve(reaction) except fluxanalysis.FluxBalanceError: logger.info( 'FBA is infeasible, marking {} as essential'.format( testing_reaction)) - model_test.limits[testing_reaction].bound = saved_bounds + c.delete() essential.add(testing_reaction) continue @@ -1069,7 +1066,7 @@ def run(self): reaction, p.get_flux(reaction))) if p.get_flux(reaction) < flux_threshold: - model_test.limits[testing_reaction].bounds = saved_bounds + c.delete() essential.add(testing_reaction) logger.info('Reaction {} was essential'.format( testing_reaction)) @@ -1077,14 +1074,9 @@ def run(self): deleted.add(testing_reaction) logger.info('Reaction {} was deleted'.format(testing_reaction)) - if self._args.exchange: - for reaction_id in sorted(exchange): - value = 0 if reaction_id in deleted else 1 - print('{}\t{}'.format(reaction_id, value)) - else: - for reaction_id in sorted(self._mm.reactions): - value = 0 if reaction_id in deleted else 1 - print('{}\t{}'.format(reaction_id, value)) + for reaction_id in sorted(reactions): + value = 0 if reaction_id in deleted else 1 + print('{}\t{}'.format(reaction_id, value)) class RobustnessCommand(SolverCommandMixin, Command): From 34c208998e5acb8d1b1a0f82c333e945cd5b26f5 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 18 Aug 2015 15:47:47 -0400 Subject: [PATCH 12/46] tests: Add placeholder test_command.py for command testing --- psamm/tests/test_command.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 psamm/tests/test_command.py diff --git a/psamm/tests/test_command.py b/psamm/tests/test_command.py new file mode 100644 index 00000000..58d8e9a7 --- /dev/null +++ b/psamm/tests/test_command.py @@ -0,0 +1,18 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2015 Jon Lund Steffensen + +from psamm import command From ce807c2e784769a0dcd7c3dd184d35cee678474b Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Wed, 19 Aug 2015 13:01:56 -0400 Subject: [PATCH 13/46] sbml: Log a warning when a reaction references a non-existent compound In strict mode this would result in a ParseError being raised but in non-strict mode no warning was previously generated. --- psamm/datasource/sbml.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/psamm/datasource/sbml.py b/psamm/datasource/sbml.py index 890e043c..b60afe21 100644 --- a/psamm/datasource/sbml.py +++ b/psamm/datasource/sbml.py @@ -22,12 +22,16 @@ from fractions import Fraction from functools import partial from itertools import count +import logging from six import itervalues, iteritems from ..reaction import Reaction, Compound +logger = logging.getLogger(__name__) + + # Level 1 namespaces SBML_NS_L1 = 'http://www.sbml.org/sbml/level1' @@ -160,12 +164,13 @@ def __init__(self, reader, root): # Skip boundary species when ignoring these continue except KeyError: + message = ('Reaction {} references non-existent' + ' species {}'.format(self._id, species_id)) if not self._reader._strict: # In non-strict mode simply skip these references + logger.warn(message) continue - raise ParseError( - 'Reaction {} references non-existent' - ' species {}'.format(self._id, species_id)) + raise ParseError(message) species_id, species_comp = ( species_entry.id, species_entry.compartment) From d366746883192b388b30bfcead2b3fb904b453ff Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Wed, 19 Aug 2015 13:09:48 -0400 Subject: [PATCH 14/46] docs: Reword install instructions to be more clear The install instructions are changed to more clearly highlight the procedure for installing psamm into a Virtualenv. The commands are presented in the natural order, which should help someone scanning the document. The use of a Virtualenv is not just presented as an option but as a recommendation. --- docs/install.rst | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/docs/install.rst b/docs/install.rst index 306429b3..01628003 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -2,32 +2,35 @@ Install ======= -The Python module can be installed using ``pip``. This will typically require -*root* permissions. +PSAMM can be installed using the Python package installer ``pip``. We recommend +that you use a `Virtualenv`_ when installing PSAMM. First, create a Virtualenv +in your project directory and activate the environment. On Linux/OSX the +following terminal commands can be used: .. code-block:: shell - $ pip install psamm + $ virtualenv env + $ source env/bin/activate -Another option that does not require *root* permissions is to use a -`Virtualenv`_. First set up a new environment in your project directory and -activate it: +Then, install PSAMM using the ``pip`` command: .. code-block:: shell - $ virtualenv env - $ . env/bin/activate + (env) $ pip install psamm + +When returning to the project from a new terminal window, simply reactivate +the environment by running -Now the Python module can be installed in the virtual environment using the -``pip`` command without requiring *root* permissions. When returning to the -project, simply reactivate the environment by running the second command. +.. code-block:: shell + + $ source env/bin/activate The *psamm-import* tool is developed in a separate Git repository. After installing PSAMM, the *psamm-import* tool can be installed using: .. code-block:: shell - $ pip install git+https://github.com/zhanglab/psamm-import.git + (env) $ pip install git+https://github.com/zhanglab/psamm-import.git Dependencies ------------ @@ -42,6 +45,8 @@ one to work. The LP solver *Cplex* is the preferred solver. The rational solver *QSopt_ex* does not support MILP problems which means that some analyses require *Cplex*. +.. _install-cplex: + Cplex ----- From 7598ac0a089c2a26130d3bfb941ed75f2354306f Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Wed, 19 Aug 2015 13:14:10 -0400 Subject: [PATCH 15/46] docs: Add FAQ entry about common Cplex install error --- docs/faq.rst | 26 ++++++++++++++++++++++++++ docs/index.rst | 1 + 2 files changed, 27 insertions(+) create mode 100644 docs/faq.rst diff --git a/docs/faq.rst b/docs/faq.rst new file mode 100644 index 00000000..7366bb74 --- /dev/null +++ b/docs/faq.rst @@ -0,0 +1,26 @@ + +FAQ +=== + +**When I run PSAMM it exits with the error "No solvers available". How can I +fix this?** + +This means that PSAMM is searching for a linear programming solver but was not +able to find one. This can occur even when the Cplex solver is installed +because the Cplex Python-bindings have to be installed separately from the +main Cplex package (see :ref:`install-cplex`). Also, if using a `Virtualenv` +with Python, the Cplex Python-bindings must be installed `in the Virtualenv`. +The bindings will `not` be available in the Virtualenv if they are installed +globally. + +To check whether the Cplex Python-bindings are correctly installed use the +following command: + +.. code-block:: shell + + (env) $ pip list + +This will show a list of Python packages that are available. The package +``cplex`` will appear in this list if the Cplex Python-bindings are correctly +installed. If the package does `not` appear in the list, follow the instuctions +at :ref:`install-cplex` to install the package. diff --git a/docs/index.rst b/docs/index.rst index a5635bf7..ab3de6c5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -12,6 +12,7 @@ Contents: file_format commands development + faq api references From 0d2ff29385ca17a7b4676575d63edeba12e50d75 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 20 Aug 2015 10:49:21 -0400 Subject: [PATCH 16/46] docs: Show (env) prefix on install commands of solvers Show (env) prefix on commands to install Cplex and QSopt_ex solvers to make it clear that the packages need to be installed in the Virtualenv. --- docs/install.rst | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/docs/install.rst b/docs/install.rst index 01628003..59d2bd6e 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -57,7 +57,13 @@ the virtual environment: 1. Locate the directory where Cplex was installed (e.g. ``/path/to/IBM/ILOG/CPLEX_StudioXXX``). 2. Locate the appropriate subdirectory based on your platform: ``cplex/python/`` (e.g. ``cplex/python/x86-64_osx``). -3. Use ``pip`` to install the package from this directory: ``pip install /path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/`` +3. Use ``pip`` to install the package from this directory using the following + command. + +.. code-block:: shell + + (env) $ pip install \ + /path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/ QSopt_ex -------- @@ -68,7 +74,7 @@ can be installed using ``pip``: .. code-block:: shell - $ pip install python-qsoptex + (env) $ pip install python-qsoptex .. _Virtualenv: https://virtualenv.pypa.io/ .. _python-qsoptex: https://pypi.python.org/pypi/python-qsoptex From 77013f823916a19c4ed9cab376dcf995312abba0 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 20 Aug 2015 00:08:10 -0400 Subject: [PATCH 17/46] lpsolver: Add Gurobi LP solver interface --- psamm/lpsolver/generic.py | 13 ++ psamm/lpsolver/gurobi.py | 290 ++++++++++++++++++++++++++++++++++++++ tox.ini | 4 +- 3 files changed, 306 insertions(+), 1 deletion(-) create mode 100644 psamm/lpsolver/gurobi.py diff --git a/psamm/lpsolver/generic.py b/psamm/lpsolver/generic.py index 0e69c945..80823059 100644 --- a/psamm/lpsolver/generic.py +++ b/psamm/lpsolver/generic.py @@ -56,6 +56,19 @@ except ImportError: pass +# Try to load Gurobi solver +try: + from . import gurobi + _solvers.append({ + 'class': gurobi.Solver, + 'name': 'gurobi', + 'integer': True, + 'rational': False, + 'priority': 9 + }) +except ImportError: + pass + class RequirementsError(Exception): """Error resolving solver requirements""" diff --git a/psamm/lpsolver/gurobi.py b/psamm/lpsolver/gurobi.py new file mode 100644 index 00000000..9f85ccdc --- /dev/null +++ b/psamm/lpsolver/gurobi.py @@ -0,0 +1,290 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2015 Jon Lund Steffensen + +"""Linear programming solver using Gurobi.""" + +from __future__ import absolute_import + +import logging +from itertools import repeat, count +import numbers + +from six.moves import zip + +import gurobipy + +from .lp import Solver as BaseSolver +from .lp import Constraint as BaseConstraint +from .lp import Problem as BaseProblem +from .lp import Result as BaseResult +from .lp import (VariableSet, Expression, Relation, + ObjectiveSense, VariableType, + InvalidResultError) + +# Module-level logging +logger = logging.getLogger(__name__) + + +class Solver(BaseSolver): + """Represents an LP-solver using Gurobi.""" + + def create_problem(self, **kwargs): + """Create a new LP-problem using the solver.""" + return Problem(**kwargs) + + +class Problem(BaseProblem): + """Represents an LP-problem of a gurobi.Solver.""" + + VARTYPE_MAP = { + VariableType.Continuous: gurobipy.GRB.CONTINUOUS, + VariableType.Binary: gurobipy.GRB.BINARY, + VariableType.Integer: gurobipy.GRB.INTEGER + } + + CONSTR_SENSE_MAP = { + Relation.Equals: gurobipy.GRB.EQUAL, + Relation.Greater: gurobipy.GRB.GREATER_EQUAL, + Relation.Less: gurobipy.GRB.LESS_EQUAL + } + + OBJ_SENSE_MAP = { + ObjectiveSense.Maximize: gurobipy.GRB.MAXIMIZE, + ObjectiveSense.Minimize: gurobipy.GRB.MINIMIZE + } + + def __init__(self, **kwargs): + self._p = gurobipy.Model() + + # TODO Logging should use the Python logging system like every other + # part of PSAMM. Unfortunately, Gurobi does not make this easy as the + # logging output is directed to stdout/(stderr?) and/or a file. For + # now, let's just disable logging to console so we don't mess up any + # of our output. This can be done with the LogToConsole parameter which + # disables output to console EXCEPT that when this parameter is set + # a log message is generated explaining that the parameter was + # changed... We are fortunate enough that this does not happen when + # using the OutputFlag parameter but instead we lose the log file. + self._p.params.OutputFlag = 0 + + self._variables = {} + self._var_names = ('x'+str(i) for i in count(1)) + self._constr_names = ('c'+str(i) for i in count(1)) + + self._result = None + + @property + def gurobi(self): + """The underlying Gurobi Model object.""" + return self._p + + def define(self, *names, **kwargs): + """Define variable in the problem. + + Variables must be defined before they can be accessed by var() or + set(). This function takes keyword arguments lower and upper to define + the bounds of the variable (default: -inf to inf). The keyword argument + types can be used to select the type of the variable (Continuous + (default), Binary or Integer). Setting any variables different than + Continuous will turn the problem into an MILP problem. + """ + names = tuple(names) + lower = kwargs.get('lower', None) + upper = kwargs.get('upper', None) + vartype = kwargs.get('types', None) + + # Repeat values if a scalar is given + if lower is None or isinstance(lower, numbers.Number): + lower = repeat(lower, len(names)) + if upper is None or isinstance(upper, numbers.Number): + upper = repeat(upper, len(names)) + if vartype is None or vartype in ( + VariableType.Continuous, VariableType.Binary, + VariableType.Integer): + vartype = repeat(vartype, len(names)) + + lp_names = tuple(next(self._var_names) for name in names) + + # Assign default values + lower = (-gurobipy.GRB.INFINITY if value is None else value + for value in lower) + upper = (gurobipy.GRB.INFINITY if value is None else value + for value in upper) + vartype = tuple(VariableType.Continuous if value is None else value + for value in vartype) + + self._variables.update(zip(names, lp_names)) + for name, lb, ub, vt in zip(lp_names, lower, upper, vartype): + self._p.addVar(name=name, lb=lb, ub=ub, vtype=self.VARTYPE_MAP[vt]) + + self._p.update() + + def var(self, name): + """Return the variable as an expression.""" + if name not in self._variables: + raise ValueError('Undefined variable: {}'.format(name)) + return Expression({name: 1}) + + def set(self, names): + """Return the set of variables as an expression.""" + names = tuple(names) + if any(name not in self._variables for name in names): + raise ValueError('Undefined variables: {}'.format( + set(names) - set(self._variables))) + return Expression({VariableSet(names): 1}) + + def _grb_expr_from_value_set(self, value_set): + return gurobipy.LinExpr( + [(val, self._p.getVarByName(self._variables[var])) + for var, val in value_set]) + + def _add_constraints(self, relation): + """Add the given relation as one or more constraints. + + Return a list of the names of the constraints added. + """ + if relation.sense in ( + Relation.StrictlyGreater, Relation.StrictlyLess): + raise ValueError( + 'Strict relations are invalid in LP-problems:' + ' {}'.format(relation)) + + names = [] + expression = relation.expression + for value_set in expression.value_sets(): + name = next(self._constr_names) + self._p.addConstr( + self._grb_expr_from_value_set(value_set), + self.CONSTR_SENSE_MAP[relation.sense], -expression.offset, + name) + names.append(name) + + self._p.update() + + return names + + def add_linear_constraints(self, *relations): + """Add constraints to the problem. + + Each constraint is represented by a Relation, and the + expression in that relation can be a set expression. + """ + constraints = [] + + for relation in relations: + if isinstance(relation, bool): + # A bool in place of a relation is accepted to mean + # a relation that does not involve any variables and + # has therefore been evaluated to a truth-value (e.g + # '0 == 0' or '2 >= 3'). + if not relation: + raise ValueError('Unsatisfiable relation added') + constraints.append(Constraint(self, None)) + else: + for name in self._add_constraints(relation): + constraints.append(Constraint(self, name)) + + return constraints + + def set_linear_objective(self, expression): + """Set linear objective of problem.""" + + if isinstance(expression, numbers.Number): + # Allow expressions with no variables as objective, + # represented as a number + expression = Expression() + + self._p.setObjective( + self._grb_expr_from_value_set(expression.values())) + + def set_objective_sense(self, sense): + """Set type of problem (maximize or minimize).""" + + if sense not in (ObjectiveSense.Minimize, ObjectiveSense.Maximize): + raise ValueError('Invalid objective sense') + + self._p.ModelSense = self.OBJ_SENSE_MAP[sense] + + def solve(self, sense=None): + """Solve problem.""" + if sense is not None: + self.set_objective_sense(sense) + self._p.optimize() + + self._result = Result(self) + return self._result + + @property + def result(self): + return self._result + + +class Constraint(BaseConstraint): + """Represents a constraint in a gurobi.Problem.""" + + def __init__(self, prob, name): + self._prob = prob + self._name = name + + def delete(self): + if self._name is not None: + self._prob._p.remove(self._prob._p.getConstrByName(self._name)) + + +class Result(BaseResult): + """Represents the solution to a gurobi.Problem. + + This object will be returned from the gurobi.Problem.solve() method or by + accessing the gurobi.Problem.result property after solving a problem. This + class should not be instantiated manually. + + Result will evaluate to a boolean according to the success of the + solution, so checking the truth value of the result will immediately + indicate whether solving was successful. + """ + + def __init__(self, prob): + self._problem = prob + + def _check_valid(self): + if self._problem.result != self: + raise InvalidResultError() + + @property + def success(self): + """Return boolean indicating whether a solution was found.""" + self._check_valid() + return self._problem._p.Status == gurobipy.GRB.OPTIMAL + + @property + def status(self): + """Return string indicating the error encountered on failure.""" + self._check_valid() + return str(self._problem._p.Status) + + def get_value(self, expression): + """Return value of expression.""" + + self._check_valid() + if isinstance(expression, Expression): + return sum(self._problem._p.getVarByName( + self._problem._variables[var]).x * value + for var, value in expression.values()) + elif expression not in self._problem._variables: + raise ValueError('Unknown expression: {}'.format(expression)) + return self._problem._p.getVarByName( + self._problem._variables[expression]).x diff --git a/tox.ini b/tox.ini index 96a6063c..0b21c85e 100644 --- a/tox.ini +++ b/tox.ini @@ -1,6 +1,6 @@ [tox] envlist = - py27-{nosolver,cplex,qsoptex}, + py27-{nosolver,cplex,qsoptex,gurobi}, py34-{nosolver,qsoptex}, coverage, flake, @@ -15,10 +15,12 @@ setenv = nosolver: PSAMM_SOLVER=nosolver cplex: PSAMM_SOLVER=cplex qsoptex: PSAMM_SOLVER=qsoptex + gurobi: PSAMM_SOLVER=gurobi deps = coverage cplex: {env:CPLEX_PYTHON_PACKAGE} qsoptex: python-qsoptex>=0.4 + gurobi: {env:GUROBI_PYTHON_PACKAGE} commands = coverage run -p --branch --omit={envdir}/*,env/*,psamm/tests/*,setup.py \ ./setup.py test From cbadf9ef8807e322db925d39b6c28b2a0596dcb2 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Fri, 21 Aug 2015 13:19:28 -0400 Subject: [PATCH 18/46] tests: Change fluxanalysis tests to use assertAlmostEqual The output from the LP solvers are floating point numbers so should not be compared using assertEqual. This fixes test failures occuring for the Gurobi solver. --- psamm/tests/test_fluxanalysis.py | 60 ++++++++++++++++---------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/psamm/tests/test_fluxanalysis.py b/psamm/tests/test_fluxanalysis.py index c68358d0..e1ad0805 100755 --- a/psamm/tests/test_fluxanalysis.py +++ b/psamm/tests/test_fluxanalysis.py @@ -45,29 +45,29 @@ def setUp(self): def test_flux_balance_rxn_1(self): fluxes = dict(fluxanalysis.flux_balance( self.model, 'rxn_1', tfba=False, solver=self.solver)) - self.assertEqual(fluxes['rxn_1'], 500) - self.assertEqual(fluxes['rxn_2'], 0) - self.assertEqual(fluxes['rxn_6'], 1000) + self.assertAlmostEqual(fluxes['rxn_1'], 500) + self.assertAlmostEqual(fluxes['rxn_2'], 0) + self.assertAlmostEqual(fluxes['rxn_6'], 1000) def test_flux_balance_rxn_2(self): fluxes = dict(fluxanalysis.flux_balance( self.model, 'rxn_2', tfba=False, solver=self.solver)) - self.assertEqual(fluxes['rxn_2'], 0) + self.assertAlmostEqual(fluxes['rxn_2'], 0) def test_flux_balance_rxn_3(self): fluxes = dict(fluxanalysis.flux_balance( self.model, 'rxn_3', tfba=False, solver=self.solver)) - self.assertEqual(fluxes['rxn_1'], 500) - self.assertEqual(fluxes['rxn_2'], 0) - self.assertEqual(fluxes['rxn_3'], 1000) - self.assertEqual(fluxes['rxn_6'], 1000) + self.assertAlmostEqual(fluxes['rxn_1'], 500) + self.assertAlmostEqual(fluxes['rxn_2'], 0) + self.assertAlmostEqual(fluxes['rxn_3'], 1000) + self.assertAlmostEqual(fluxes['rxn_6'], 1000) def test_flux_balance_rxn_6(self): fluxes = dict(fluxanalysis.flux_balance( self.model, 'rxn_6', tfba=False, solver=self.solver)) - self.assertEqual(fluxes['rxn_1'], 500) - self.assertEqual(fluxes['rxn_2'], 0) - self.assertEqual(fluxes['rxn_6'], 1000) + self.assertAlmostEqual(fluxes['rxn_1'], 500) + self.assertAlmostEqual(fluxes['rxn_2'], 0) + self.assertAlmostEqual(fluxes['rxn_6'], 1000) class TestFluxBalanceThermodynamic(unittest.TestCase): @@ -94,11 +94,11 @@ def setUp(self): def test_flux_balance_tfba_exchange_d(self): fluxes = dict(fluxanalysis.flux_balance( self.model, 'ex_D', tfba=True, solver=self.solver)) - self.assertEquals(fluxes['ex_A'], -10) - self.assertEquals(fluxes['ex_D'], 10) - self.assertEquals(fluxes['rxn_2'], 10) - self.assertEquals(fluxes['rxn_4'], 0) - self.assertEquals(fluxes['rxn_5'], 0) + self.assertAlmostEqual(fluxes['ex_A'], -10) + self.assertAlmostEqual(fluxes['ex_D'], 10) + self.assertAlmostEqual(fluxes['rxn_2'], 10) + self.assertAlmostEqual(fluxes['rxn_4'], 0) + self.assertAlmostEqual(fluxes['rxn_5'], 0) class TestFluxVariability(unittest.TestCase): @@ -127,15 +127,15 @@ def test_flux_variability(self): self.model, self.model.reactions, {'rxn_6': 200}, tfba=False, solver=self.solver)) - self.assertEqual(fluxes['rxn_1'][0], 100) + self.assertAlmostEqual(fluxes['rxn_1'][0], 100) - self.assertEqual(fluxes['rxn_2'][0], 0) - self.assertEqual(fluxes['rxn_2'][1], 0) + self.assertAlmostEqual(fluxes['rxn_2'][0], 0) + self.assertAlmostEqual(fluxes['rxn_2'][1], 0) - self.assertEqual(fluxes['rxn_5'][0], 0) - self.assertEqual(fluxes['rxn_5'][1], 100) + self.assertAlmostEqual(fluxes['rxn_5'][0], 0) + self.assertAlmostEqual(fluxes['rxn_5'][1], 100) - self.assertEqual(fluxes['rxn_6'][0], 200) + self.assertAlmostEqual(fluxes['rxn_6'][0], 200) self.assertGreater(fluxes['rxn_7'][1], 0) self.assertGreater(fluxes['rxn_8'][1], 0) @@ -145,18 +145,18 @@ def test_flux_variability_thermodynamic(self): self.model, self.model.reactions, {'rxn_6': 200}, tfba=True, solver=self.solver)) - self.assertEqual(fluxes['rxn_1'][0], 100) + self.assertAlmostEqual(fluxes['rxn_1'][0], 100) - self.assertEqual(fluxes['rxn_2'][0], 0) - self.assertEqual(fluxes['rxn_2'][1], 0) + self.assertAlmostEqual(fluxes['rxn_2'][0], 0) + self.assertAlmostEqual(fluxes['rxn_2'][1], 0) - self.assertEqual(fluxes['rxn_5'][0], 0) - self.assertEqual(fluxes['rxn_5'][1], 100) + self.assertAlmostEqual(fluxes['rxn_5'][0], 0) + self.assertAlmostEqual(fluxes['rxn_5'][1], 100) - self.assertEqual(fluxes['rxn_6'][0], 200) + self.assertAlmostEqual(fluxes['rxn_6'][0], 200) - self.assertEqual(fluxes['rxn_7'][1], 0) - self.assertEqual(fluxes['rxn_8'][1], 0) + self.assertAlmostEqual(fluxes['rxn_7'][1], 0) + self.assertAlmostEqual(fluxes['rxn_8'][1], 0) class TestFluxConsistency(unittest.TestCase): From 89c7751fbe7787dce0072e9d98ba423465953889 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Fri, 21 Aug 2015 13:11:39 -0400 Subject: [PATCH 19/46] docs: Document install procedure for Gurobi solver --- docs/install.rst | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/docs/install.rst b/docs/install.rst index 306429b3..0c53f368 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -32,15 +32,17 @@ installing PSAMM, the *psamm-import* tool can be installed using: Dependencies ------------ -- Linear programming solver (*Cplex*, *QSopt_ex*) +- Linear programming solver (*Cplex*, *Gurobi*, or *QSopt_ex*) - PyYAML (for reading the native model format) - NumPy (optional; model matrix can be exported to NumPy matrix if available) PyYAML is installed automatically when PSAMM is installed through ``pip``. The linear programming solver is not strictly required but most analyses require -one to work. The LP solver *Cplex* is the preferred solver. The rational solver -*QSopt_ex* does not support MILP problems which means that some analyses -require *Cplex*. +one to work. The LP solver *Cplex* is the preferred solver. We recently added +support for the LP solver *Gurobi*. + +The rational solver *QSopt_ex* does not support MILP problems which means that +some analyses require one of the other solvers. Cplex ----- @@ -54,6 +56,21 @@ the virtual environment: ``cplex/python/`` (e.g. ``cplex/python/x86-64_osx``). 3. Use ``pip`` to install the package from this directory: ``pip install /path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/`` +Gurobi +------ + +The Gurobi Python bindings will have to be installed into the virtualenv. After +activating the virtualenv: + +1. Locate the directory where the Guropi python bindings were installed. For + example, on OSX this directory is ``/Library/gurobiXXX/mac64`` where ``XXX`` + is a version code. +2. Use ``pip`` to install the package from this directory. For example: + +.. code-block:: shell + + (env) $ pip install /Library/gurobi604/mac64 + QSopt_ex -------- From 00a1f0acf929dbc2a8004f1ba7e169aedc7b23a7 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Fri, 21 Aug 2015 13:15:39 -0400 Subject: [PATCH 20/46] tox: Pass DYLD_LIBRARY_PATH to test environment This makes it easier to set up tox on OSX where the DYLD_LIBRARY_PATH is needed when using a non-Apple Python. --- tox.ini | 1 + 1 file changed, 1 insertion(+) diff --git a/tox.ini b/tox.ini index 0b21c85e..844a6c1c 100644 --- a/tox.ini +++ b/tox.ini @@ -21,6 +21,7 @@ deps = cplex: {env:CPLEX_PYTHON_PACKAGE} qsoptex: python-qsoptex>=0.4 gurobi: {env:GUROBI_PYTHON_PACKAGE} +passenv = DYLD_LIBRARY_PATH commands = coverage run -p --branch --omit={envdir}/*,env/*,psamm/tests/*,setup.py \ ./setup.py test From 0e1521b630f22075c5e28a8e9fcc0184d2a3683b Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 24 Aug 2015 09:57:08 -0400 Subject: [PATCH 21/46] Add more coverage files to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index db4561ea..353a76db 100644 --- a/.gitignore +++ b/.gitignore @@ -36,6 +36,7 @@ pip-delete-this-directory.txt htmlcov/ .tox/ .coverage +.coverage.* .cache nosetests.xml coverage.xml From eeda383afb8805584d66ca04a0816f9d7806160e Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 25 Aug 2015 10:50:32 -0400 Subject: [PATCH 22/46] docs: Add documentation for gurobi module --- docs/api.rst | 1 + docs/conf.py | 2 +- docs/lpsolver_gurobi.rst | 6 ++++++ 3 files changed, 8 insertions(+), 1 deletion(-) create mode 100644 docs/lpsolver_gurobi.rst diff --git a/docs/api.rst b/docs/api.rst index 131495bb..8f0395d2 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -22,6 +22,7 @@ PSAMM API lpsolver_lp lpsolver_generic lpsolver_cplex + lpsolver_gurobi lpsolver_qsoptex massconsistency metabolicmodel diff --git a/docs/conf.py b/docs/conf.py index 2caba659..5b89ad95 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -29,7 +29,7 @@ # Mock optional modules to allow autodoc to run even in the absense of these # modules. -MOCK_MODULES = ['cplex', 'qsoptex'] +MOCK_MODULES = ['cplex', 'qsoptex', 'gurobipy'] for module in MOCK_MODULES: sys.modules[module] = mock.Mock() diff --git a/docs/lpsolver_gurobi.rst b/docs/lpsolver_gurobi.rst new file mode 100644 index 00000000..e55c912a --- /dev/null +++ b/docs/lpsolver_gurobi.rst @@ -0,0 +1,6 @@ + +``psamm.lpsolver.gurobi`` -- Gurobi LP solver +================================================= + +.. automodule:: psamm.lpsolver.gurobi + :members: From 963641e2fbb49b84ca4e73319500eb07800fa84e Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 25 Aug 2015 10:54:09 -0400 Subject: [PATCH 23/46] docs: Move API documentation to subdirectory --- docs/api.rst | 24 ++---------------------- docs/{ => api}/command.rst | 0 docs/{ => api}/database.rst | 0 docs/{ => api}/datasource_kegg.rst | 0 docs/{ => api}/datasource_misc.rst | 0 docs/{ => api}/datasource_modelseed.rst | 0 docs/{ => api}/datasource_native.rst | 0 docs/{ => api}/datasource_sbml.rst | 0 docs/{ => api}/expression_affine.rst | 0 docs/{ => api}/expression_boolean.rst | 0 docs/{ => api}/fastcore.rst | 0 docs/{ => api}/fluxanalysis.rst | 0 docs/{ => api}/fluxcoupling.rst | 0 docs/{ => api}/formula.rst | 0 docs/{ => api}/gapfill.rst | 0 docs/{ => api}/lpsolver_cplex.rst | 0 docs/{ => api}/lpsolver_generic.rst | 0 docs/{ => api}/lpsolver_gurobi.rst | 0 docs/{ => api}/lpsolver_lp.rst | 0 docs/{ => api}/lpsolver_qsoptex.rst | 0 docs/{ => api}/massconsistency.rst | 0 docs/{ => api}/metabolicmodel.rst | 0 docs/{ => api}/reaction.rst | 0 23 files changed, 2 insertions(+), 22 deletions(-) rename docs/{ => api}/command.rst (100%) rename docs/{ => api}/database.rst (100%) rename docs/{ => api}/datasource_kegg.rst (100%) rename docs/{ => api}/datasource_misc.rst (100%) rename docs/{ => api}/datasource_modelseed.rst (100%) rename docs/{ => api}/datasource_native.rst (100%) rename docs/{ => api}/datasource_sbml.rst (100%) rename docs/{ => api}/expression_affine.rst (100%) rename docs/{ => api}/expression_boolean.rst (100%) rename docs/{ => api}/fastcore.rst (100%) rename docs/{ => api}/fluxanalysis.rst (100%) rename docs/{ => api}/fluxcoupling.rst (100%) rename docs/{ => api}/formula.rst (100%) rename docs/{ => api}/gapfill.rst (100%) rename docs/{ => api}/lpsolver_cplex.rst (100%) rename docs/{ => api}/lpsolver_generic.rst (100%) rename docs/{ => api}/lpsolver_gurobi.rst (100%) rename docs/{ => api}/lpsolver_lp.rst (100%) rename docs/{ => api}/lpsolver_qsoptex.rst (100%) rename docs/{ => api}/massconsistency.rst (100%) rename docs/{ => api}/metabolicmodel.rst (100%) rename docs/{ => api}/reaction.rst (100%) diff --git a/docs/api.rst b/docs/api.rst index 8f0395d2..1f6c5bee 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -4,26 +4,6 @@ PSAMM API .. toctree:: :maxdepth: 2 + :glob: - command - database - datasource_kegg - datasource_misc - datasource_modelseed - datasource_native - datasource_sbml - expression_affine - expression_boolean - fastcore - fluxanalysis - fluxcoupling - formula - gapfill - lpsolver_lp - lpsolver_generic - lpsolver_cplex - lpsolver_gurobi - lpsolver_qsoptex - massconsistency - metabolicmodel - reaction + api/* diff --git a/docs/command.rst b/docs/api/command.rst similarity index 100% rename from docs/command.rst rename to docs/api/command.rst diff --git a/docs/database.rst b/docs/api/database.rst similarity index 100% rename from docs/database.rst rename to docs/api/database.rst diff --git a/docs/datasource_kegg.rst b/docs/api/datasource_kegg.rst similarity index 100% rename from docs/datasource_kegg.rst rename to docs/api/datasource_kegg.rst diff --git a/docs/datasource_misc.rst b/docs/api/datasource_misc.rst similarity index 100% rename from docs/datasource_misc.rst rename to docs/api/datasource_misc.rst diff --git a/docs/datasource_modelseed.rst b/docs/api/datasource_modelseed.rst similarity index 100% rename from docs/datasource_modelseed.rst rename to docs/api/datasource_modelseed.rst diff --git a/docs/datasource_native.rst b/docs/api/datasource_native.rst similarity index 100% rename from docs/datasource_native.rst rename to docs/api/datasource_native.rst diff --git a/docs/datasource_sbml.rst b/docs/api/datasource_sbml.rst similarity index 100% rename from docs/datasource_sbml.rst rename to docs/api/datasource_sbml.rst diff --git a/docs/expression_affine.rst b/docs/api/expression_affine.rst similarity index 100% rename from docs/expression_affine.rst rename to docs/api/expression_affine.rst diff --git a/docs/expression_boolean.rst b/docs/api/expression_boolean.rst similarity index 100% rename from docs/expression_boolean.rst rename to docs/api/expression_boolean.rst diff --git a/docs/fastcore.rst b/docs/api/fastcore.rst similarity index 100% rename from docs/fastcore.rst rename to docs/api/fastcore.rst diff --git a/docs/fluxanalysis.rst b/docs/api/fluxanalysis.rst similarity index 100% rename from docs/fluxanalysis.rst rename to docs/api/fluxanalysis.rst diff --git a/docs/fluxcoupling.rst b/docs/api/fluxcoupling.rst similarity index 100% rename from docs/fluxcoupling.rst rename to docs/api/fluxcoupling.rst diff --git a/docs/formula.rst b/docs/api/formula.rst similarity index 100% rename from docs/formula.rst rename to docs/api/formula.rst diff --git a/docs/gapfill.rst b/docs/api/gapfill.rst similarity index 100% rename from docs/gapfill.rst rename to docs/api/gapfill.rst diff --git a/docs/lpsolver_cplex.rst b/docs/api/lpsolver_cplex.rst similarity index 100% rename from docs/lpsolver_cplex.rst rename to docs/api/lpsolver_cplex.rst diff --git a/docs/lpsolver_generic.rst b/docs/api/lpsolver_generic.rst similarity index 100% rename from docs/lpsolver_generic.rst rename to docs/api/lpsolver_generic.rst diff --git a/docs/lpsolver_gurobi.rst b/docs/api/lpsolver_gurobi.rst similarity index 100% rename from docs/lpsolver_gurobi.rst rename to docs/api/lpsolver_gurobi.rst diff --git a/docs/lpsolver_lp.rst b/docs/api/lpsolver_lp.rst similarity index 100% rename from docs/lpsolver_lp.rst rename to docs/api/lpsolver_lp.rst diff --git a/docs/lpsolver_qsoptex.rst b/docs/api/lpsolver_qsoptex.rst similarity index 100% rename from docs/lpsolver_qsoptex.rst rename to docs/api/lpsolver_qsoptex.rst diff --git a/docs/massconsistency.rst b/docs/api/massconsistency.rst similarity index 100% rename from docs/massconsistency.rst rename to docs/api/massconsistency.rst diff --git a/docs/metabolicmodel.rst b/docs/api/metabolicmodel.rst similarity index 100% rename from docs/metabolicmodel.rst rename to docs/api/metabolicmodel.rst diff --git a/docs/reaction.rst b/docs/api/reaction.rst similarity index 100% rename from docs/reaction.rst rename to docs/api/reaction.rst From 2bba241f1a563c9d49b665677cbd431bb61817ae Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 24 Aug 2015 16:09:12 -0400 Subject: [PATCH 24/46] tox: Add Python 3.3 to tox environment --- tox.ini | 1 + 1 file changed, 1 insertion(+) diff --git a/tox.ini b/tox.ini index 844a6c1c..fa0527a6 100644 --- a/tox.ini +++ b/tox.ini @@ -1,6 +1,7 @@ [tox] envlist = py27-{nosolver,cplex,qsoptex,gurobi}, + py33-{nosolver,qsoptex}, py34-{nosolver,qsoptex}, coverage, flake, From 766709bf0933ed328a4e5bbd0dba7350172d94f4 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 24 Aug 2015 16:10:37 -0400 Subject: [PATCH 25/46] setup.py: Add Python 3.3 classification --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index 8494ef75..29b7b6cb 100755 --- a/setup.py +++ b/setup.py @@ -37,6 +37,7 @@ 'Development Status :: 4 - Beta', 'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)', 'Programming Language :: Python :: 2.7', + 'Programming Language :: Python :: 3.3' 'Programming Language :: Python :: 3.4' ], From 6a1be0c7036bbac70663268f4189263860884778 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 10:15:33 -0400 Subject: [PATCH 26/46] util: Suppress error output when calling git describe --- psamm/util.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/psamm/util.py b/psamm/util.py index 2074dc4b..83ad67c3 100644 --- a/psamm/util.py +++ b/psamm/util.py @@ -121,7 +121,8 @@ def git_try_describe(repo_path): """ try: p = subprocess.Popen(['git', 'describe', '--always', '--dirty'], - cwd=repo_path, stdout=subprocess.PIPE) + cwd=repo_path, stdout=subprocess.PIPE, + stderr=subprocess.PIPE) output, _ = p.communicate() except: return None From 27cd4e91f4f8d91ce7c873f7a11304fdbdb6b3a3 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 10:30:08 -0400 Subject: [PATCH 27/46] Add command line tool to list available LP solvers --- psamm/lpsolver/generic.py | 36 +++++++++++++++++++++++++++++------- setup.py | 3 ++- 2 files changed, 31 insertions(+), 8 deletions(-) diff --git a/psamm/lpsolver/generic.py b/psamm/lpsolver/generic.py index 80823059..a310de57 100644 --- a/psamm/lpsolver/generic.py +++ b/psamm/lpsolver/generic.py @@ -17,7 +17,7 @@ """Generic interface to LP solver instantiation.""" -from __future__ import absolute_import +from __future__ import absolute_import, print_function import os import logging @@ -28,6 +28,7 @@ logger = logging.getLogger(__name__) _solvers = [] +_solver_import_errors = {} # Try to load Cplex solver @@ -40,8 +41,8 @@ 'rational': False, 'priority': 10 }) -except ImportError: - pass +except ImportError as e: + _solver_import_errors['cplex'] = str(e) # Try to load QSopt_ex solver try: @@ -53,8 +54,8 @@ 'rational': True, 'priority': 5 }) -except ImportError: - pass +except ImportError as e: + _solver_import_errors['qsoptex'] = str(e) # Try to load Gurobi solver try: @@ -66,8 +67,8 @@ 'rational': False, 'priority': 9 }) -except ImportError: - pass +except ImportError as e: + _solver_import_errors['gurobi'] = str(e) class RequirementsError(Exception): @@ -139,3 +140,24 @@ def parse_solver_setting(s): value = float(value) return key, value + + +def list_solvers(): + """Print list of solvers.""" + if len(_solvers) > 0: + print('Available solvers:') + for solver in _solvers: + print('Name: {}'.format(solver['name'])) + print('Priority: {}'.format(solver['priority'])) + print('MILP (integer) problem support: {}'.format( + solver['integer'])) + print('Rational solution: {}'.format(solver['rational'])) + print('Class: {}'.format(solver['class'])) + print() + else: + print('No solvers available!') + print() + + print('Unavailable solvers:') + for solver, error in iteritems(_solver_import_errors): + print('{}: {}'.format(solver, error)) diff --git a/setup.py b/setup.py index 29b7b6cb..90dbf8c2 100755 --- a/setup.py +++ b/setup.py @@ -44,7 +44,8 @@ packages=find_packages(), entry_points = { 'console_scripts': [ - 'psamm-model = psamm.command:main' + 'psamm-model = psamm.command:main', + 'psamm-list-lpsolvers = psamm.lpsolver.generic:list_solvers' ] }, From cc7e0d167e4733743c7d96151129b0c34ac3bb7c Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 12:10:05 -0400 Subject: [PATCH 28/46] lpsolver: Improve output of psamm-list-lpsolvers Changed command to use argparse. This improves the help messages and allows additional requirements to be specified on the command line. The PSAMM_SOLVER environment variable is also used to decide priority of solvers. --- psamm/lpsolver/generic.py | 84 +++++++++++++++++++++++++++++++++------ 1 file changed, 71 insertions(+), 13 deletions(-) diff --git a/psamm/lpsolver/generic.py b/psamm/lpsolver/generic.py index a310de57..b3f91f04 100644 --- a/psamm/lpsolver/generic.py +++ b/psamm/lpsolver/generic.py @@ -19,6 +19,7 @@ from __future__ import absolute_import, print_function +import argparse import os import logging @@ -75,6 +76,17 @@ class RequirementsError(Exception): """Error resolving solver requirements""" +def filter_solvers(solvers, requirements): + """Yield solvers that fullfil the requirements.""" + for solver in solvers: + for req, value in iteritems(requirements): + if (req in ('integer', 'rational', 'name') and + (req not in solver or solver[req] != value)): + break + else: + yield solver + + class Solver(BaseSolver): """Generic solver interface based on requirements @@ -87,15 +99,12 @@ class Solver(BaseSolver): """ def __init__(self, **kwargs): - solvers = _solvers - if len(solvers) == 0: + if len(_solvers) == 0: raise RequirementsError('No solvers available') self._requirements = {key: value for key, value in iteritems(kwargs) if value is not None} - for req, value in iteritems(self._requirements): - if req in ('integer', 'rational', 'name'): - solvers = [s for s in solvers if req in s and s[req] == value] + solvers = filter_solvers(_solvers, self._requirements) # Obtain solver priority from environment variable, if specified. priority = {} @@ -143,10 +152,52 @@ def parse_solver_setting(s): def list_solvers(): - """Print list of solvers.""" - if len(_solvers) > 0: - print('Available solvers:') - for solver in _solvers: + """Entry point for listing available solvers.""" + parser = argparse.ArgumentParser( + description='''List LP solver available in PSAMM. This will produce a + list of all of the available LP solvers in prioritized + order. Addtional requirements can be imposed with the + arguments (e.g. integer=yes to select only solvers that + support MILP problems). The list will also be influenced + by the PSAMM_SOLVER environment variable which can be + used to only allow specific solvers (e.g. + PSAMM_SOLVER=cplex).''') + parser.add_argument( + 'requirement', nargs='*', type=str, + help='Additional requirements on the selected solvers') + args = parser.parse_args() + + requirements = {} + for arg in args.requirement: + try: + key, value = parse_solver_setting(arg) + except ValueError as e: + parser.error(str(e)) + else: + requirements[key] = value + + solvers = list(filter_solvers(_solvers, requirements)) + solver_names = set(solver['name'] for solver in solvers) + + # Obtain solver priority from environment variable, if specified. + priority = {} + if 'PSAMM_SOLVER' in os.environ: + names = os.environ['PSAMM_SOLVER'].split(',') + for i, solver_name in enumerate(names): + priority[solver_name] = len(names) - i + solvers = [s for s in solvers if s['name'] in priority] + solver_names = set(priority) + else: + # Use built-in priorities + for solver in solvers: + priority[solver['name']] = solver['priority'] + + solvers = sorted(solvers, key=lambda s: priority.get(s['name'], 0), + reverse=True) + + if len(solvers) > 0: + print('Prioritized solvers:') + for solver in solvers: print('Name: {}'.format(solver['name'])) print('Priority: {}'.format(solver['priority'])) print('MILP (integer) problem support: {}'.format( @@ -155,9 +206,16 @@ def list_solvers(): print('Class: {}'.format(solver['class'])) print() else: - print('No solvers available!') + print('No solvers fullfil the requirements!') print() - print('Unavailable solvers:') - for solver, error in iteritems(_solver_import_errors): - print('{}: {}'.format(solver, error)) + filtered_solvers_count = len(_solvers) - len(solvers) + if filtered_solvers_count > 0 or len(_solver_import_errors) > 0: + print('Unavailable solvers:') + for solver in _solvers: + if solver['name'] not in solver_names: + print('{}: Does not fullfil the specified requirements'.format( + solver['name'])) + + for solver, error in iteritems(_solver_import_errors): + print('{}: Error loading solver: {}'.format(solver, error)) From 7041aa57070dc282645114b685c92b8a3b22c996 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 14:01:34 -0400 Subject: [PATCH 29/46] Fix exception when PSAMM_SOLVER is not set --- psamm/lpsolver/generic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psamm/lpsolver/generic.py b/psamm/lpsolver/generic.py index b3f91f04..1cf4257d 100644 --- a/psamm/lpsolver/generic.py +++ b/psamm/lpsolver/generic.py @@ -104,7 +104,7 @@ def __init__(self, **kwargs): self._requirements = {key: value for key, value in iteritems(kwargs) if value is not None} - solvers = filter_solvers(_solvers, self._requirements) + solvers = list(filter_solvers(_solvers, self._requirements)) # Obtain solver priority from environment variable, if specified. priority = {} From f181a005efe76c9af3e32c0749a3fb5ae7d183f4 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Wed, 26 Aug 2015 13:43:56 -0400 Subject: [PATCH 30/46] command: Separate commands into psamm.commands package Splits the command module into one module for each command. These modules exist as submodules of psamm.commands. The Command subclasses in these modules are references using the entry point mechanism in setup.py thereby making the commands discoverable. --- psamm/command.py | 1230 +------------------------------- psamm/commands/__init__.py | 0 psamm/commands/chargecheck.py | 82 +++ psamm/commands/console.py | 69 ++ psamm/commands/fastgapfill.py | 150 ++++ psamm/commands/fba.py | 114 +++ psamm/commands/fluxcheck.py | 126 ++++ psamm/commands/fluxcoupling.py | 138 ++++ psamm/commands/formulacheck.py | 98 +++ psamm/commands/fva.py | 75 ++ psamm/commands/gapfill.py | 80 +++ psamm/commands/masscheck.py | 115 +++ psamm/commands/randomsparse.py | 136 ++++ psamm/commands/robustness.py | 125 ++++ psamm/commands/sbmlexport.py | 33 + psamm/commands/search.py | 145 ++++ setup.py | 16 + 17 files changed, 1521 insertions(+), 1211 deletions(-) create mode 100644 psamm/commands/__init__.py create mode 100644 psamm/commands/chargecheck.py create mode 100644 psamm/commands/console.py create mode 100644 psamm/commands/fastgapfill.py create mode 100644 psamm/commands/fba.py create mode 100644 psamm/commands/fluxcheck.py create mode 100644 psamm/commands/fluxcoupling.py create mode 100644 psamm/commands/formulacheck.py create mode 100644 psamm/commands/fva.py create mode 100644 psamm/commands/gapfill.py create mode 100644 psamm/commands/masscheck.py create mode 100644 psamm/commands/randomsparse.py create mode 100644 psamm/commands/robustness.py create mode 100644 psamm/commands/sbmlexport.py create mode 100644 psamm/commands/search.py diff --git a/psamm/command.py b/psamm/command.py index 4ec23110..86a94b09 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -14,44 +14,32 @@ # along with PSAMM. If not, see . # # Copyright 2014-2015 Jon Lund Steffensen -# Copyright 2015 Keith Dufault-Thompson """Command line interface. Each command in the command line interface is implemented as a subclass of -:class:`Command`. Commands are automatically discovered based on this -inheritance so new commands are added by subclassing :class:`Command`. +:class:`Command`. Commands are also referenced from ``setup.py`` using the +entry point mechanism which allows the commands to be automatically +discovered. The :func:`.main` function is the entry point of command line interface. """ -from __future__ import print_function - -import sys import os import argparse -import operator -import re import logging -import random -import math import abc +import pkg_resources +from six import add_metaclass, iteritems + from . import __version__ as package_version -from .formula import Formula -from .gapfill import gapfind, gapfill -from .database import DictDatabase -from .metabolicmodel import MetabolicModel -from .reaction import Compound from .datasource.native import NativeModel -from .datasource import sbml -from . import fluxanalysis, fluxcoupling, massconsistency, fastcore +from .metabolicmodel import MetabolicModel +from .database import DictDatabase from .lpsolver import generic from . import util -from six import add_metaclass, iteritems - -# Module-level logging logger = logging.getLogger(__name__) @@ -146,1181 +134,6 @@ def _get_solver(self, **kwargs): return generic.Solver(**solver_args) -class ChargeBalanceCommand(Command): - """Check whether compound charge in a given database or model is balanced. - - Balanced reactions are those reactions where the total charge - is consistent on the left and right side of the reaction equation. - Reactions that are not balanced will be printed out. - """ - - name = 'chargecheck' - title = 'Check charge balance on a model or database' - - def run(self): - """Run charge balance command""" - - # Load compound information - compound_name = {} - compound_charge = {} - for compound in self._model.parse_compounds(): - compound_name[compound.id] = ( - compound.name if compound.name is not None else compound.id) - if hasattr(compound, 'charge') and compound.charge is not None: - compound_charge[compound.id] = compound.charge - - # Create a set of known charge-inconsistent reactions - exchange = set() - for reaction_id in self._mm.reactions: - if self._mm.is_exchange(reaction_id): - exchange.add(reaction_id) - - def reaction_charges(reaction_id): - for compound, value in self._mm.get_reaction_values(reaction_id): - charge = compound_charge.get(compound.name, float('nan')) - yield charge * float(value) - - count = 0 - unbalanced = 0 - unchecked = 0 - for reaction in sorted(self._mm.reactions): - if reaction in exchange: - continue - - count += 1 - charge = sum(reaction_charges(reaction)) - if math.isnan(charge): - logger.debug('Not checking reaction {};' - ' missing charge'.format(reaction)) - unchecked += 1 - elif charge != 0: - unbalanced += 1 - rx = self._mm.get_reaction(reaction) - rxt = rx.translated_compounds( - lambda x: compound_name.get(x, x)) - print('{}\t{}\t{}'.format(reaction, charge, rxt)) - - logger.info('Unbalanced reactions: {}/{}'.format(unbalanced, count)) - logger.info('Unchecked reactions due to missing charge: {}/{}'.format( - unchecked, count)) - - -class ConsoleCommand(Command): - """Start an interactive Python console with the given model loaded.""" - - name = 'console' - title = 'Start Python console with metabolic model loaded' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--type', choices=('python', 'ipython', 'ipython-kernel'), - default='python', help='type of console to open') - - def open_python(self, message, namespace): - """Open interactive python console""" - - # Importing readline will in some cases print weird escape - # characters to stdout. To avoid this we only import readline - # and related packages at this point when we are certain - # they are needed. - from code import InteractiveConsole - import readline - import rlcompleter - - readline.set_completer(rlcompleter.Completer(namespace).complete) - readline.parse_and_bind('tab: complete') - console = InteractiveConsole(namespace) - console.interact(message) - - def open_ipython(self, message, namespace): - from IPython.terminal.embed import InteractiveShellEmbed - console = InteractiveShellEmbed(user_ns=namespace, banner2=message) - console() - - def open_ipython_kernel(self, message, namespace): - from IPython import embed_kernel - embed_kernel(local_ns=namespace) - - def run(self): - message = ('Native model has been loaded into: "model"\n' + - 'Metabolic model has been loaded into: "mm"') - namespace = {'model': self._model, 'mm': self._mm} - console_type = self._args.type - - if console_type == 'python': - self.open_python(message, namespace) - elif console_type == 'ipython': - self.open_ipython(message, namespace) - elif console_type == 'ipython-kernel': - self.open_ipython_kernel(message, namespace) - - -class FastGapFillCommand(SolverCommandMixin, Command): - """Run FastGapFill algorithm on a metabolic model.""" - - name = 'fastgapfill' - title = 'Run FastGapFill on a metabolic model' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--penalty', metavar='file', type=argparse.FileType('r'), - help='List of penalty scores for database reactions') - parser.add_argument( - '--db-weight', metavar='weight', type=float, - help='Default weight for database reactions') - parser.add_argument( - '--tp-weight', metavar='weight', type=float, - help='Default weight for transport reactions') - parser.add_argument( - '--ex-weight', metavar='weight', type=float, - help='Default weight for exchange reactions') - parser.add_argument( - '--epsilon', type=float, help='Threshold for Fastcore', - default=1e-5) - parser.add_argument( - '--no-tfba', help='Disable thermodynamic constraints on FBA', - action='store_true') - parser.add_argument( - 'reaction', help='Reaction to maximize', nargs='?') - super(FastGapFillCommand, cls).init_parser(parser) - - def run(self): - """Run FastGapFill command""" - - # Create solver - enable_tfba = not self._args.no_tfba - if enable_tfba: - solver = self._get_solver(integer=True) - else: - solver = self._get_solver() - - # Load compound information - compound_name = {} - for compound in self._model.parse_compounds(): - compound_name[compound.id] = ( - compound.name if compound.name is not None else compound.id) - - epsilon = self._args.epsilon - model_compartments = set(self._mm.compartments) - - # Add exchange and transport reactions to database - model_complete = self._mm.copy() - logger.info('Adding database, exchange and transport reactions') - db_added = model_complete.add_all_database_reactions( - model_compartments) - ex_added = model_complete.add_all_exchange_reactions() - tp_added = model_complete.add_all_transport_reactions() - - # Add penalty weights on reactions - weights = {} - if self._args.db_weight is not None: - weights.update((rxnid, self._args.db_weight) for rxnid in db_added) - if self._args.tp_weight is not None: - weights.update((rxnid, self._args.tp_weight) for rxnid in tp_added) - if self._args.ex_weight is not None: - weights.update((rxnid, self._args.ex_weight) for rxnid in ex_added) - - if self._args.penalty is not None: - for line in self._args.penalty: - line, _, comment = line.partition('#') - line = line.strip() - if line == '': - continue - rxnid, weight = line.split(None, 1) - weights[rxnid] = float(weight) - - # Run Fastcore and print the induced reaction set - logger.info('Calculating Fastcore induced set on model') - core = set(self._mm.reactions) - - induced = fastcore.fastcore(model_complete, core, epsilon, - weights=weights, solver=solver) - logger.info('Result: |A| = {}, A = {}'.format(len(induced), induced)) - added_reactions = induced - core - logger.info('Extended: |E| = {}, E = {}'.format( - len(added_reactions), added_reactions)) - - if self._args.reaction is not None: - maximized_reaction = self._args.reaction - else: - maximized_reaction = self._model.get_biomass_reaction() - if maximized_reaction is None: - raise CommandError('The maximized reaction was not specified') - - if not not self._mm.has_reaction(maximized_reaction): - raise CommandError( - 'The biomass reaction is not a valid model' - ' reaction: {}'.format(maximized_reaction)) - - logger.info('Flux balance on induced model maximizing {}'.format( - maximized_reaction)) - model_induced = self._mm.copy() - for rxnid in induced: - model_induced.add_reaction(rxnid) - for rxnid, flux in sorted(fluxanalysis.flux_balance( - model_induced, maximized_reaction, tfba=enable_tfba, - solver=solver)): - reaction_class = 'Dbase' - weight = weights.get(rxnid, 1) - if self._mm.has_reaction(rxnid): - reaction_class = 'Model' - weight = 0 - rx = model_complete.get_reaction(rxnid) - rxt = rx.translated_compounds(lambda x: compound_name.get(x, x)) - print('{}\t{}\t{}\t{}\t{}'.format( - rxnid, reaction_class, weight, flux, rxt)) - - logger.info('Calculating Fastcc consistent subset of induced model') - consistent_core = fastcore.fastcc_consistent_subset( - model_induced, epsilon, solver=solver) - logger.info('Result: |A| = {}, A = {}'.format( - len(consistent_core), consistent_core)) - removed_reactions = set(model_induced.reactions) - consistent_core - logger.info('Removed: |R| = {}, R = {}'.format( - len(removed_reactions), removed_reactions)) - - -class FluxBalanceCommand(SolverCommandMixin, Command): - """Run flux balance analysis on a metabolic model.""" - - name = 'fba' - title = 'Run flux balance analysis on a metabolic model' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--no-tfba', help='Disable thermodynamic constraints on FBA', - action='store_true') - parser.add_argument( - '--epsilon', type=float, help='Threshold for flux minimization', - default=1e-5) - parser.add_argument('reaction', help='Reaction to maximize', nargs='?') - super(FluxBalanceCommand, cls).init_parser(parser) - - def run(self): - """Run flux analysis command""" - - # Load compound information - compound_name = {} - for compound in self._model.parse_compounds(): - if 'name' in compound.properties: - compound_name[compound.id] = compound.properties['name'] - elif compound.id not in compound_name: - compound_name[compound.id] = compound.id - - if self._args.reaction is not None: - reaction = self._args.reaction - else: - reaction = self._model.get_biomass_reaction() - if reaction is None: - raise CommandError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise CommandError('Specified reaction is not in model: {}'.format( - reaction)) - - if self._args.no_tfba: - result = self.run_fba_minimized(reaction) - else: - result = self.run_tfba(reaction) - - optimum = None - for reaction_id, fba_flux, flux in sorted(result): - rx = self._mm.get_reaction(reaction_id) - print('{}\t{}\t{}\t{}'.format( - reaction_id, fba_flux, flux, - rx.translated_compounds(lambda x: compound_name.get(x, x)))) - # Remember flux of requested reaction - if reaction_id == reaction: - optimum = flux - - logger.info('Maximum flux: {}'.format(optimum)) - - def run_fba_minimized(self, reaction): - """Run normal FBA and flux minimization on model, then print output""" - - solver = self._get_solver() - fba_fluxes = dict(fluxanalysis.flux_balance( - self._mm, reaction, tfba=False, solver=solver)) - optimum = fba_fluxes[reaction] - epsilon = self._args.epsilon - - # Run flux minimization - fmin_fluxes = dict(fluxanalysis.flux_minimization( - self._mm, {reaction: optimum}, solver=solver)) - count = 0 - for reaction_id, flux in iteritems(fmin_fluxes): - if fba_fluxes[reaction_id] - epsilon > flux: - count += 1 - yield reaction_id, fba_fluxes[reaction_id], flux - logger.info('Minimized reactions: {}'.format(count)) - - def run_tfba(self, reaction): - """Run FBA and tFBA on model""" - - solver = self._get_solver(integer=True) - - fba_fluxes = dict(fluxanalysis.flux_balance( - self._mm, reaction, tfba=False, solver=solver)) - fluxes = dict(fluxanalysis.flux_balance( - self._mm, reaction, tfba=True, solver=solver)) - - for reaction_id, flux in iteritems(fluxes): - yield reaction_id, fba_fluxes[reaction_id], flux - - -class FluxConsistencyCommand(SolverCommandMixin, Command): - """Check that reactions are flux consistent in a model. - - A reaction is flux consistent if there exists any steady-state flux - solution where the flux of the given reaction is non-zero. The - bounds on the exchange reactions can be removed when performing the - consistency check by providing the ``--unrestricted`` option. - """ - - name = 'fluxcheck' - title = 'Check that the model is flux consistent' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--fastcore', help='Enable use of Fastcore algorithm', - action='store_true') - parser.add_argument( - '--no-tfba', - help='Disable thermodynamic constraints on flux check', - action='store_true') - parser.add_argument( - '--epsilon', type=float, help='Flux threshold', - default=1e-5) - parser.add_argument( - '--unrestricted', action='store_true', - help='Remove limits on exchange reactions before checking' - ) - super(FluxConsistencyCommand, cls).init_parser(parser) - - def run(self): - """Run flux consistency check command""" - - # Load compound information - compound_name = {} - for compound in self._model.parse_compounds(): - if 'name' in compound.properties: - compound_name[compound.id] = compound.properties['name'] - elif compound.id not in compound_name: - compound_name[compound.id] = compound.id - - epsilon = self._args.epsilon - - if self._args.unrestricted: - # 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 - - if self._args.fastcore: - solver = self._get_solver() - inconsistent = set(fastcore.fastcc( - self._mm, epsilon, solver=solver)) - else: - enable_tfba = not self._args.no_tfba - if enable_tfba: - solver = self._get_solver(integer=True) - else: - solver = self._get_solver() - inconsistent = set( - fluxanalysis.consistency_check( - self._mm, self._mm.reactions, epsilon, - tfba=enable_tfba, solver=solver)) - - # Count the number of reactions that are fixed at zero. While these - # reactions are still inconsistent, they are inconsistent because they - # have been explicitly disabled. - disabled_exchange = 0 - disabled_internal = 0 - - count_exchange = 0 - total_exchange = 0 - - count_internal = 0 - total_internal = 0 - - # Print result - for reaction in sorted(self._mm.reactions): - disabled = self._mm.limits[reaction].bounds == (0, 0) - - if self._mm.is_exchange(reaction): - total_exchange += 1 - count_exchange += int(reaction in inconsistent) - disabled_exchange += int(disabled) - else: - total_internal += 1 - count_internal += int(reaction in inconsistent) - disabled_internal += int(disabled) - - if reaction in inconsistent: - rx = self._mm.get_reaction(reaction) - rxt = rx.translated_compounds( - lambda x: compound_name.get(x, x)) - print('{}\t{}'.format(reaction, rxt)) - - logger.info('Model has {}/{} inconsistent internal reactions' - ' ({} disabled by user)'.format( - count_internal, total_internal, disabled_internal)) - logger.info('Model has {}/{} inconsistent exchange reactions' - ' ({} disabled by user)'.format( - count_exchange, total_exchange, disabled_exchange)) - - -class FluxCouplingCommand(SolverCommandMixin, Command): - """Find flux coupled reactions in the model. - - This identifies any reaction pairs where the flux of one reaction - constrains the flux of another reaction. The reactions can be coupled in - three distinct ways depending on the ratio between the reaction fluxes. - The reactions can be fully coupled (the ratio is static and non-zero); - partially coupled (the ratio is bounded and non-zero); or directionally - coupled (the ratio is non-zero). - """ - - name = 'fluxcoupling' - title = 'Find flux coupled reactions in model' - - def run(self): - solver = self._get_solver() - - max_reaction = self._model.get_biomass_reaction() - if max_reaction is None: - raise CommandError('The biomass reaction was not specified') - - fba_fluxes = dict(fluxanalysis.flux_balance( - self._mm, max_reaction, tfba=False, solver=solver)) - optimum = fba_fluxes[max_reaction] - - self._fcp = fluxcoupling.FluxCouplingProblem( - self._mm, {max_reaction: 0.999 * optimum}, solver) - - self._coupled = {} - self._groups = [] - - reactions = sorted(self._mm.reactions) - for i, reaction1 in enumerate(reactions): - if reaction1 in self._coupled: - continue - - for reaction2 in reactions[i+1:]: - if (reaction2 in self._coupled and - (self._coupled[reaction2] == - self._coupled.get(reaction1))): - continue - - self._check_reactions(reaction1, reaction2) - - logger.info('Coupled groups:') - for i, group in enumerate(self._groups): - if group is not None: - logger.info('{}: {}'.format(i, ', '.join(sorted(group)))) - - def _check_reactions(self, reaction1, reaction2): - logger.debug('Solving for {}, {}'.format(reaction1, reaction2)) - lower, upper = self._fcp.solve(reaction1, reaction2) - - logger.debug('Result: {}, {}'.format(lower, upper)) - - coupling = fluxcoupling.classify_coupling((lower, upper)) - if coupling in (fluxcoupling.CouplingClass.DirectionalForward, - fluxcoupling.CouplingClass.DirectionalReverse): - text = 'Directional, v1 / v2 in [{}, {}]'.format(lower, upper) - if (coupling == fluxcoupling.CouplingClass.DirectionalReverse and - not self._mm.is_reversible(reaction1) and lower == 0.0): - return - elif coupling == fluxcoupling.CouplingClass.Full: - text = 'Full: v1 / v2 = {}'.format(lower) - self._couple_reactions(reaction1, reaction2) - elif coupling == fluxcoupling.CouplingClass.Partial: - text = 'Partial: v1 / v2 in [{}; {}]'.format(lower, upper) - self._couple_reactions(reaction1, reaction2) - else: - return - - print('{}\t{}\t{}\t{}\t{}'.format( - reaction1, reaction2, lower, upper, text)) - - def _couple_reactions(self, reaction1, reaction2): - logger.debug('Couple {} and {}'.format(reaction1, reaction2)) - - if reaction1 in self._coupled and reaction2 in self._coupled: - if self._coupled[reaction1] == self._coupled[reaction2]: - return - logger.debug('Merge groups {}, {}'.format( - self._coupled[reaction1], self._coupled[reaction2])) - group_index = len(self._groups) - group = (self._groups[self._coupled[reaction1]] | - self._groups[self._coupled[reaction2]]) - logger.debug('New group is {}: {}'.format( - group_index, sorted(group))) - self._groups.append(group) - self._groups[self._coupled[reaction1]] = None - self._groups[self._coupled[reaction2]] = None - for reaction in group: - self._coupled[reaction] = group_index - elif reaction1 in self._coupled: - group_index = self._coupled[reaction1] - group = self._groups[group_index] - logger.debug('Put {} into existing group {}: {}'.format( - reaction2, self._coupled[reaction1], group)) - self._coupled[reaction2] = group_index - group.add(reaction2) - elif reaction2 in self._coupled: - group_index = self._coupled[reaction2] - group = self._groups[group_index] - logger.debug('Put {} into existing group {}: {}'.format( - reaction1, self._coupled[reaction2], group)) - self._coupled[reaction1] = group_index - group.add(reaction1) - else: - group = set([reaction1, reaction2]) - group_index = len(self._groups) - logger.debug('Creating new group {}'.format(group_index)) - self._coupled[reaction1] = group_index - self._coupled[reaction2] = group_index - self._groups.append(group) - - -class FluxVariabilityCommand(SolverCommandMixin, Command): - """Run flux variablity analysis on a metabolic model.""" - - name = 'fva' - title = 'Run flux variability analysis on a metabolic model' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--no-tfba', help='Disable thermodynamic constraints on FVA', - action='store_true') - parser.add_argument('reaction', help='Reaction to maximize', nargs='?') - super(FluxVariabilityCommand, cls).init_parser(parser) - - def run(self): - """Run flux variability command""" - - # Load compound information - compound_name = {} - for compound in self._model.parse_compounds(): - if 'name' in compound.properties: - compound_name[compound.id] = compound.properties['name'] - elif compound.id not in compound_name: - compound_name[compound.id] = compound.id - - if self._args.reaction is not None: - reaction = self._args.reaction - else: - reaction = self._model.get_biomass_reaction() - if reaction is None: - raise CommandError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise CommandError('Specified reaction is not in model: {}'.format( - reaction)) - - enable_tfba = not self._args.no_tfba - if enable_tfba: - solver = self._get_solver(integer=True) - else: - solver = self._get_solver() - - fba_fluxes = dict(fluxanalysis.flux_balance( - self._mm, reaction, tfba=False, solver=solver)) - optimum = fba_fluxes[reaction] - - flux_bounds = fluxanalysis.flux_variability( - self._mm, sorted(self._mm.reactions), {reaction: optimum}, - tfba=enable_tfba, solver=solver) - for reaction_id, bounds in flux_bounds: - rx = self._mm.get_reaction(reaction_id) - rxt = rx.translated_compounds(lambda x: compound_name.get(x, x)) - print('{}\t{}\t{}\t{}'.format( - reaction_id, bounds[0], bounds[1], rxt)) - - -class FormulaBalanceCommand(Command): - """Check whether reactions in a model are elementally balanced. - - Balanced reactions are those reactions where the number of elements - (atoms) is consistent on the left and right side of the reaction equation. - Reactions that are not balanced will be printed out. - """ - - name = 'formulacheck' - title = 'Check formula balance on a model or database' - - def run(self): - """Run formula balance command""" - - # Mapping from compound id to formula - compound_formula = {} - for compound in self._model.parse_compounds(): - # Only cpd11632 (Photon) is allowed to have an empty formula. - if compound.formula is None: - if compound.id == 'cpd11632': - compound_formula[compound.id] = Formula() - else: - try: - f = Formula.parse(compound.formula).flattened() - compound_formula[compound.id] = f - except ValueError as e: - logger.warning( - 'Error parsing {}: {}'.format(compound.formula, e)) - - # Create a set of known mass-inconsistent reactions - exchange = set() - for reaction_id in self._mm.reactions: - if self._mm.is_exchange(reaction_id): - exchange.add(reaction_id) - - def multiply_formula(compound_list): - for compound, count in compound_list: - yield count * compound_formula.get(compound.name, Formula()) - - count = 0 - unbalanced = 0 - unchecked = 0 - for reaction in self._mm.reactions: - if reaction in exchange: - continue - - count += 1 - rx = self._mm.get_reaction(reaction) - - # Skip reaction if any compounds have undefined formula - for compound, _ in rx.compounds: - if compound.name not in compound_formula: - unchecked += 1 - break - else: - left_form = reduce( - operator.or_, multiply_formula(rx.left), Formula()) - right_form = reduce( - operator.or_, multiply_formula(rx.right), Formula()) - - if right_form != left_form: - unbalanced += 1 - right_missing, left_missing = Formula.balance( - right_form, left_form) - - print('{}\t{}\t{}\t{}\t{}'.format( - reaction, left_form, right_form, - left_missing, right_missing)) - - logger.info('Unbalanced reactions: {}/{}'.format(unbalanced, count)) - logger.info('Unchecked reactions due to missing formula: {}/{}'.format( - unchecked, count)) - - -class GapFillCommand(SolverCommandMixin, Command): - """Run GapFind and GapFill on a metabolic model.""" - - name = 'gapfill' - title = 'Run GapFind and GapFill on a metabolic model' - - def run(self): - """Run GapFill command""" - - # Load compound information - compound_name = {} - for compound in self._model.parse_compounds(): - if 'name' in compound.properties: - compound_name[compound.id] = compound.properties['name'] - elif compound.id not in compound_name: - compound_name[compound.id] = compound.id - - model_compartments = set(self._mm.compartments) - - solver = self._get_solver(integer=True) - - # Run GapFind on model - logger.info('Searching for blocked compounds') - blocked = set(compound for compound in gapfind(self._mm, solver=solver) - if compound.compartment is not 'e') - if len(blocked) > 0: - logger.info('Blocked compounds') - for compound in blocked: - print(compound.translate(lambda x: compound_name.get(x, x))) - - if len(blocked) > 0: - # Add exchange and transport reactions to database - model_complete = self._mm.copy() - logger.info('Adding database, exchange and transport reactions') - model_complete.add_all_database_reactions(model_compartments) - model_complete.add_all_exchange_reactions() - model_complete.add_all_transport_reactions() - - logger.info('Searching for reactions to fill gaps') - added_reactions, reversed_reactions = gapfill( - model_complete, self._mm.reactions, blocked, solver=solver) - - for rxnid in added_reactions: - rx = model_complete.get_reaction(rxnid) - rxt = rx.translated_compounds( - lambda x: compound_name.get(x, x)) - print('{}\t{}\t{}'.format(rxnid, 'Add', rxt)) - - for rxnid in reversed_reactions: - rx = model_complete.get_reaction(rxnid) - rxt = rx.translated_compounds( - lambda x: compound_name.get(x, x)) - print('{}\t{}\t{}'.format(rxnid, 'Reverse', rxt)) - else: - logger.info('No blocked compounds found') - - -class MassConsistencyCommand(SolverCommandMixin, Command): - """Check whether a model is mass consistent.""" - - name = 'masscheck' - title = 'Run mass consistency check on a database' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--exclude', metavar='reaction', action='append', type=str, - default=[], help='Exclude reaction from mass consistency') - parser.add_argument( - '--epsilon', type=float, help='Mass threshold', - default=1e-5) - parser.add_argument( - '--checked', metavar='reaction', action='append', type=str, - default=[], help='Mark reaction as already checked (no residual)') - super(MassConsistencyCommand, cls).init_parser(parser) - - def run(self): - # Load compound information - compound_name = {} - zeromass = set() - for compound in self._model.parse_compounds(): - if 'name' in compound.properties: - compound_name[compound.id] = compound.properties['name'] - elif compound.id not in compound_name: - compound_name[compound.id] = compound.id - - if compound.properties.get('zeromass', False): - zeromass.add(compound.id) - - # Create a set of known mass-inconsistent reactions - exchange = set() - for reaction_id in self._mm.reactions: - if self._mm.is_exchange(reaction_id): - exchange.add(reaction_id) - - # Other reactions to exclude from consistency check - exclude = set(self._args.exclude) - - # Add biomass reaction to be excluded - biomass_reaction = self._model.get_biomass_reaction() - if biomass_reaction is not None: - exclude.add(biomass_reaction) - - # Create set of checked reactions - checked = set(self._args.checked) - - solver = self._get_solver() - - known_inconsistent = exclude | exchange - - logger.info('Mass consistency on database') - epsilon = self._args.epsilon - compound_iter = massconsistency.check_compound_consistency( - self._mm, solver, known_inconsistent, zeromass) - - good = 0 - total = 0 - for compound, mass in sorted(compound_iter, key=lambda x: (x[1], x[0]), - reverse=True): - if mass >= 1-epsilon or compound.name in zeromass: - good += 1 - total += 1 - print('{}: {}'.format(compound.translate( - lambda x: compound_name.get(x, x)), mass)) - logger.info('Consistent compounds: {}/{}'.format(good, total)) - - logger.info('Is consistent? {}'.format( - massconsistency.is_consistent( - self._mm, solver, known_inconsistent, zeromass))) - - good = 0 - total = 0 - reaction_iter, compound_iter = ( - massconsistency.check_reaction_consistency( - self._mm, solver=solver, exchange=known_inconsistent, - checked=checked, zeromass=zeromass)) - for reaction_id, residual in sorted( - reaction_iter, key=lambda x: abs(x[1]), reverse=True): - total += 1 - if abs(residual) >= epsilon: - reaction = self._mm.get_reaction(reaction_id) - rxt = reaction.translated_compounds( - lambda x: compound_name.get(x, x)) - print('{}\t{}\t{}'.format(reaction_id, residual, rxt)) - else: - good += 1 - logger.info('Consistent reactions: {}/{}'.format(good, total)) - - -class RandomSparseNetworkCommand(SolverCommandMixin, Command): - """Find random minimal network of model. - - Given a reaction to optimize and a threshold, delete reactions randomly - until the flux of the reaction to optimize falls under the threshold. - Keep deleting reactions until no more reactions can be deleted. By default - this uses standard FBA (not tFBA). Since the internal fluxes are irrelevant - the FBA and tFBA are equivalent for this purpose. - - The threshold can be specified as an absolute flux (e.g. '1.23') or a - relative flux of the full model flux (e.g. '40.5%'). - """ - - name = 'randomsparse' - title = 'Generate a random sparse network of model reactions' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--tfba', help='Enable thermodynamic constraints on FBA', - action='store_true', default=False) - parser.add_argument( - '--reaction', help='Reaction to maximize') - parser.add_argument( - 'threshold', help='Threshold of max reaction flux', - type=util.MaybeRelative) - parser.add_argument( - '--exchange', help='Only analyze the exchange reaction in model', - action='store_true') - super(RandomSparseNetworkCommand, cls).init_parser(parser) - - def run(self): - if self._args.reaction is not None: - reaction = self._args.reaction - else: - reaction = self._model.get_biomass_reaction() - if reaction is None: - raise CommandError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise CommandError('Specified reaction is not in model: {}'.format( - reaction)) - - if not self._args.tfba: - fb_problem = fluxanalysis.FluxBalanceProblem - solver = self._get_solver() - else: - fb_problem = fluxanalysis.FluxBalanceTDProblem - solver = self._get_solver(integer=True) - - p = fb_problem(self._mm, solver) - - threshold = self._args.threshold - if threshold.relative: - p.solve(reaction) - threshold.reference = p.get_flux(reaction) - - flux_threshold = float(threshold) - - logger.info('Flux threshold for {} is {}'.format( - reaction, flux_threshold)) - - essential = {reaction} - deleted = set() - if self._args.exchange: - reactions = set() - for reaction_id in self._mm.reactions: - if self._mm.is_exchange(reaction_id): - reactions.add(reaction_id) - else: - reactions = set(self._mm.reactions) - - test_set = reactions - essential - - while len(test_set) > 0: - testing_reaction = random.sample(test_set, 1)[0] - test_set.remove(testing_reaction) - - flux_var = p.get_flux_var(testing_reaction) - c, = p.prob.add_linear_constraints(flux_var == 0) - - logger.info('Trying FBA without reaction {}...'.format( - testing_reaction)) - - try: - p.solve(reaction) - except fluxanalysis.FluxBalanceError: - logger.info( - 'FBA is infeasible, marking {} as essential'.format( - testing_reaction)) - c.delete() - essential.add(testing_reaction) - continue - - logger.debug('Reaction {} has flux {}'.format( - reaction, p.get_flux(reaction))) - - if p.get_flux(reaction) < flux_threshold: - c.delete() - essential.add(testing_reaction) - logger.info('Reaction {} was essential'.format( - testing_reaction)) - else: - deleted.add(testing_reaction) - logger.info('Reaction {} was deleted'.format(testing_reaction)) - - for reaction_id in sorted(reactions): - value = 0 if reaction_id in deleted else 1 - print('{}\t{}'.format(reaction_id, value)) - - -class RobustnessCommand(SolverCommandMixin, Command): - """Run robustness analysis on metabolic model. - - Given a reaction to maximize and a reaction to vary, - the robustness analysis will run FBA while fixing the - reaction to vary at each iteration. The reaction will - be fixed at the specified number of steps between the - minimum and maximum flux value specified in the model. - """ - - name = 'robustness' - title = 'Run robustness analysis on a metabolic model' - - @classmethod - def init_parser(cls, parser): - parser.add_argument( - '--steps', metavar='N', type=int, default=10, - help='Number of flux value steps for varying reaction') - parser.add_argument( - '--minimum', metavar='V', type=float, - help='Minumum flux value of varying reacton') - parser.add_argument( - '--maximum', metavar='V', type=float, - help='Maximum flux value of varying reacton') - parser.add_argument( - '--no-tfba', help='Disable thermodynamic constraints on FBA', - action='store_true') - parser.add_argument( - '--reaction', help='Reaction to maximize', nargs='?') - parser.add_argument('varying', help='Reaction to vary') - super(RobustnessCommand, cls).init_parser(parser) - - def run(self): - """Run flux analysis command""" - - def run_fba_fmin(model, reaction): - fba = fluxanalysis.FluxBalanceProblem(model, solver=solver) - fba.solve(reaction) - optimum = fba.get_flux(reaction) - return fluxanalysis.flux_minimization( - model, {reaction: optimum}, solver=solver) - - def run_tfba(model, reaction): - return fluxanalysis.flux_balance( - model, reaction, tfba=True, solver=solver) - - if self._args.no_tfba: - run_fba = run_fba_fmin - solver = self._get_solver() - else: - run_fba = run_tfba - solver = self._get_solver(integer=True) - - # Load compound information - compound_name = {} - for compound in self._model.parse_compounds(): - if 'name' in compound.properties: - compound_name[compound.id] = compound.properties['name'] - elif compound.id not in compound_name: - compound_name[compound.id] = compound.id - - if self._args.reaction is not None: - reaction = self._args.reaction - else: - reaction = self._model.get_biomass_reaction() - if reaction is None: - raise CommandError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise CommandError('Specified reaction is not in model: {}'.format( - reaction)) - - varying_reaction = self._args.varying - if not self._mm.has_reaction(varying_reaction): - raise CommandError('Specified reaction is not in model: {}'.format( - varying_reaction)) - - steps = self._args.steps - if steps <= 0: - raise CommandError('Invalid number of steps: {}\n'.format(steps)) - - # Run FBA on model at different fixed flux values - flux_min = self._mm.limits[varying_reaction].lower - flux_max = self._mm.limits[varying_reaction].upper - if self._args.minimum is not None: - flux_min = self._args.minimum - if self._args.maximum is not None: - flux_max = self._args.maximum - - if flux_min > flux_max: - raise CommandError('Invalid flux range: {}, {}\n'.format( - flux_min, flux_max)) - - for i in xrange(steps): - fixed_flux = flux_min + i*(flux_max - flux_min)/float(steps-1) - test_model = self._mm.copy() - test_model.limits[varying_reaction].bounds = fixed_flux, fixed_flux - - try: - for other_reaction, flux in run_fba(test_model, reaction): - print('{}\t{}\t{}'.format( - other_reaction, fixed_flux, flux)) - except fluxanalysis.FluxBalanceError: - pass - - -class SBMLExport(Command): - """Export model as SBML file.""" - - name = 'sbmlexport' - title = 'Export model as SBML file' - - def run(self): - writer = sbml.SBMLWriter() - writer.write_model( - sys.stdout, self._mm, self._model.parse_compounds()) - - -class SearchCommand(Command): - """Search for reactions and compounds in a model.""" - - name = 'search' - title = 'Search the database of reactions or compounds' - - @classmethod - def init_parser(cls, parser): - """Initialize argument parser""" - subparsers = parser.add_subparsers(title='Search domain') - - # Compound subcommand - parser_compound = subparsers.add_parser( - 'compound', help='Search in compounds') - parser_compound.set_defaults(which='compound') - parser_compound.add_argument( - '--id', '-i', dest='id', metavar='id', type=str, - default=[], action='append', help='Compound ID') - parser_compound.add_argument( - '--name', '-n', dest='name', metavar='name', type=str, - default=[], action='append', help='Name of compound') - - # Reaction subcommand - parser_reaction = subparsers.add_parser( - 'reaction', help='Search in reactions') - parser_reaction.set_defaults(which='reaction') - parser_reaction.add_argument( - '--id', '-i', dest='id', metavar='id', type=str, - default=[], action='append', help='Reaction ID') - parser_reaction.add_argument( - '--compound', '-c', dest='compound', metavar='compound', type=str, - default=[], action='append', - help='Comma-separated list of compound IDs') - - def run(self): - """Run search command.""" - - def parse_compound(s): - """Parse a compound specification with optional compartment""" - m = re.match(r'([^\[]+)\[(\w+)\]', s) - if m is not None: - return Compound(m.group(1), compartment=m.group(2)) - return Compound(s) - - def filter_search_term(s): - return re.sub(r'[^a-z0-9]+', '', s.lower()) - - which_command = self._args.which - if which_command == 'compound': - selected_compounds = set() - - for compound in self._model.parse_compounds(): - if len(self._args.id) > 0: - if any(c == compound.id for c in self._args.id): - selected_compounds.add(compound) - continue - - if len(self._args.name) > 0: - names = set() - if 'name' in compound.properties: - names.add(compound.properties['name']) - names.update(compound.properties.get('names', [])) - names = set(filter_search_term(n) for n in names) - if any(filter_search_term(n) in names - for n in self._args.name): - selected_compounds.add(compound) - continue - - # Show results - for compound in selected_compounds: - print('ID: {}'.format(compound.id)) - if 'name' in compound.properties: - print('Name: {}'.format(compound.properties['name'])) - if 'names' in compound.properties: - print('Additional names: {}'.format( - ', '.join(compound.properties['names']))) - if 'formula' in compound.properties: - print('Formula: {}'.format(compound.properties['formula'])) - if compound.filemark is not None: - print('Parsed from: {}'.format(compound.filemark)) - print() - elif which_command == 'reaction': - selected_reactions = set() - - # Prepare translation table from compound id to name - compound_name = {} - for compound in self._model.parse_compounds(): - if 'name' in compound.properties: - compound_name[compound.id] = compound.properties['name'] - - # Prepare sets of matched compounds - search_compounds = [] - for compound_list in self._args.compound: - compound_set = set() - for compound_spec in compound_list.split(','): - compound_set.add(parse_compound(compound_spec.strip())) - search_compounds.append(compound_set) - - for reaction in self._model.parse_reactions(): - if len(self._args.id) > 0: - if any(r == reaction.id for r in self._args.id): - selected_reactions.add(reaction) - continue - - if len(search_compounds) > 0: - compounds = set(c for c, _ in reaction.equation.compounds) - if any(c.issubset(compounds) for c in search_compounds): - selected_reactions.add(reaction) - continue - - # Show results - for reaction in selected_reactions: - print('ID: {}'.format(reaction.id)) - print('Reaction (IDs): {}'.format( - reaction.equation)) - translated_equation = reaction.equation.translated_compounds( - lambda x: compound_name.get(x, x)) - if reaction.equation != translated_equation: - print('Reaction (names): {}'.format(translated_equation)) - if reaction.filemark is not None: - print('Parsed from: {}'.format(reaction.filemark)) - print() - - def main(command_class=None): """Run the command line interface with the given :class:`Command`. @@ -1354,22 +167,21 @@ def main(command_class=None): else: # Discover all available commands commands = {} - for command_class in Command.__subclasses__(): - command_name = getattr(command_class, 'name', None) - command_title = getattr(command_class, 'title', None) - if (command_name is not None and - command_title is not None and - command_name not in commands): - commands[command_name] = (command_class.title, - command_class) + for entry in pkg_resources.iter_entry_points('psamm.commands'): + canonical = entry.name.lower() + if canonical not in commands: + command_class = entry.load() + commands[canonical] = command_class + else: + logger.warning('Command {} was found more than once!'.format( + canonical)) # Create parsers for subcommands subparsers = parser.add_subparsers(title='Command') - for name, values in sorted(iteritems(commands)): - title, command_class = values - description = command_class.__doc__ + for name, command_class in sorted(iteritems(commands)): subparser = subparsers.add_parser( - name, help=title, description=description) + name, help=command_class.title, + description=command_class.__doc__) subparser.set_defaults(command=command_class) command_class.init_parser(subparser) @@ -1384,7 +196,3 @@ def main(command_class=None): command.run() except CommandError as e: parser.error(str(e)) - - -if __name__ == '__main__': - main() diff --git a/psamm/commands/__init__.py b/psamm/commands/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/psamm/commands/chargecheck.py b/psamm/commands/chargecheck.py new file mode 100644 index 00000000..6ddbe24a --- /dev/null +++ b/psamm/commands/chargecheck.py @@ -0,0 +1,82 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import math +import logging + +from ..command import Command + +logger = logging.getLogger(__name__) + + +class ChargeBalanceCommand(Command): + """Check whether compound charge in a given database or model is balanced. + + Balanced reactions are those reactions where the total charge + is consistent on the left and right side of the reaction equation. + Reactions that are not balanced will be printed out. + """ + + name = 'chargecheck' + title = 'Check charge balance on a model or database' + + def run(self): + """Run charge balance command""" + + # Load compound information + compound_name = {} + compound_charge = {} + for compound in self._model.parse_compounds(): + compound_name[compound.id] = ( + compound.name if compound.name is not None else compound.id) + if hasattr(compound, 'charge') and compound.charge is not None: + compound_charge[compound.id] = compound.charge + + # Create a set of known charge-inconsistent reactions + exchange = set() + for reaction_id in self._mm.reactions: + if self._mm.is_exchange(reaction_id): + exchange.add(reaction_id) + + def reaction_charges(reaction_id): + for compound, value in self._mm.get_reaction_values(reaction_id): + charge = compound_charge.get(compound.name, float('nan')) + yield charge * float(value) + + count = 0 + unbalanced = 0 + unchecked = 0 + for reaction in sorted(self._mm.reactions): + if reaction in exchange: + continue + + count += 1 + charge = sum(reaction_charges(reaction)) + if math.isnan(charge): + logger.debug('Not checking reaction {};' + ' missing charge'.format(reaction)) + unchecked += 1 + elif charge != 0: + unbalanced += 1 + rx = self._mm.get_reaction(reaction) + rxt = rx.translated_compounds( + lambda x: compound_name.get(x, x)) + print('{}\t{}\t{}'.format(reaction, charge, rxt)) + + logger.info('Unbalanced reactions: {}/{}'.format(unbalanced, count)) + logger.info('Unchecked reactions due to missing charge: {}/{}'.format( + unchecked, count)) diff --git a/psamm/commands/console.py b/psamm/commands/console.py new file mode 100644 index 00000000..7a04e3e1 --- /dev/null +++ b/psamm/commands/console.py @@ -0,0 +1,69 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +from ..command import Command + + +class ConsoleCommand(Command): + """Start an interactive Python console with the given model loaded.""" + + name = 'console' + title = 'Start Python console with metabolic model loaded' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--type', choices=('python', 'ipython', 'ipython-kernel'), + default='python', help='type of console to open') + + def open_python(self, message, namespace): + """Open interactive python console""" + + # Importing readline will in some cases print weird escape + # characters to stdout. To avoid this we only import readline + # and related packages at this point when we are certain + # they are needed. + from code import InteractiveConsole + import readline + import rlcompleter + + readline.set_completer(rlcompleter.Completer(namespace).complete) + readline.parse_and_bind('tab: complete') + console = InteractiveConsole(namespace) + console.interact(message) + + def open_ipython(self, message, namespace): + from IPython.terminal.embed import InteractiveShellEmbed + console = InteractiveShellEmbed(user_ns=namespace, banner2=message) + console() + + def open_ipython_kernel(self, message, namespace): + from IPython import embed_kernel + embed_kernel(local_ns=namespace) + + def run(self): + message = ('Native model has been loaded into: "model"\n' + + 'Metabolic model has been loaded into: "mm"') + namespace = {'model': self._model, 'mm': self._mm} + console_type = self._args.type + + if console_type == 'python': + self.open_python(message, namespace) + elif console_type == 'ipython': + self.open_ipython(message, namespace) + elif console_type == 'ipython-kernel': + self.open_ipython_kernel(message, namespace) diff --git a/psamm/commands/fastgapfill.py b/psamm/commands/fastgapfill.py new file mode 100644 index 00000000..d33055dd --- /dev/null +++ b/psamm/commands/fastgapfill.py @@ -0,0 +1,150 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import argparse +import logging + +from ..command import Command, SolverCommandMixin, CommandError +from .. import fastcore, fluxanalysis + +logger = logging.getLogger(__name__) + + +class FastGapFillCommand(SolverCommandMixin, Command): + """Run FastGapFill algorithm on a metabolic model.""" + + name = 'fastgapfill' + title = 'Run FastGapFill on a metabolic model' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--penalty', metavar='file', type=argparse.FileType('r'), + help='List of penalty scores for database reactions') + parser.add_argument( + '--db-weight', metavar='weight', type=float, + help='Default weight for database reactions') + parser.add_argument( + '--tp-weight', metavar='weight', type=float, + help='Default weight for transport reactions') + parser.add_argument( + '--ex-weight', metavar='weight', type=float, + help='Default weight for exchange reactions') + parser.add_argument( + '--epsilon', type=float, help='Threshold for Fastcore', + default=1e-5) + parser.add_argument( + '--no-tfba', help='Disable thermodynamic constraints on FBA', + action='store_true') + parser.add_argument( + 'reaction', help='Reaction to maximize', nargs='?') + super(FastGapFillCommand, cls).init_parser(parser) + + def run(self): + """Run FastGapFill command""" + + # Create solver + enable_tfba = not self._args.no_tfba + if enable_tfba: + solver = self._get_solver(integer=True) + else: + solver = self._get_solver() + + # Load compound information + compound_name = {} + for compound in self._model.parse_compounds(): + compound_name[compound.id] = ( + compound.name if compound.name is not None else compound.id) + + epsilon = self._args.epsilon + model_compartments = set(self._mm.compartments) + + # Add exchange and transport reactions to database + model_complete = self._mm.copy() + logger.info('Adding database, exchange and transport reactions') + db_added = model_complete.add_all_database_reactions( + model_compartments) + ex_added = model_complete.add_all_exchange_reactions() + tp_added = model_complete.add_all_transport_reactions() + + # Add penalty weights on reactions + weights = {} + if self._args.db_weight is not None: + weights.update((rxnid, self._args.db_weight) for rxnid in db_added) + if self._args.tp_weight is not None: + weights.update((rxnid, self._args.tp_weight) for rxnid in tp_added) + if self._args.ex_weight is not None: + weights.update((rxnid, self._args.ex_weight) for rxnid in ex_added) + + if self._args.penalty is not None: + for line in self._args.penalty: + line, _, comment = line.partition('#') + line = line.strip() + if line == '': + continue + rxnid, weight = line.split(None, 1) + weights[rxnid] = float(weight) + + # Run Fastcore and print the induced reaction set + logger.info('Calculating Fastcore induced set on model') + core = set(self._mm.reactions) + + induced = fastcore.fastcore(model_complete, core, epsilon, + weights=weights, solver=solver) + logger.info('Result: |A| = {}, A = {}'.format(len(induced), induced)) + added_reactions = induced - core + logger.info('Extended: |E| = {}, E = {}'.format( + len(added_reactions), added_reactions)) + + if self._args.reaction is not None: + maximized_reaction = self._args.reaction + else: + maximized_reaction = self._model.get_biomass_reaction() + if maximized_reaction is None: + raise CommandError('The maximized reaction was not specified') + + if not not self._mm.has_reaction(maximized_reaction): + raise CommandError( + 'The biomass reaction is not a valid model' + ' reaction: {}'.format(maximized_reaction)) + + logger.info('Flux balance on induced model maximizing {}'.format( + maximized_reaction)) + model_induced = self._mm.copy() + for rxnid in induced: + model_induced.add_reaction(rxnid) + for rxnid, flux in sorted(fluxanalysis.flux_balance( + model_induced, maximized_reaction, tfba=enable_tfba, + solver=solver)): + reaction_class = 'Dbase' + weight = weights.get(rxnid, 1) + if self._mm.has_reaction(rxnid): + reaction_class = 'Model' + weight = 0 + rx = model_complete.get_reaction(rxnid) + rxt = rx.translated_compounds(lambda x: compound_name.get(x, x)) + print('{}\t{}\t{}\t{}\t{}'.format( + rxnid, reaction_class, weight, flux, rxt)) + + logger.info('Calculating Fastcc consistent subset of induced model') + consistent_core = fastcore.fastcc_consistent_subset( + model_induced, epsilon, solver=solver) + logger.info('Result: |A| = {}, A = {}'.format( + len(consistent_core), consistent_core)) + removed_reactions = set(model_induced.reactions) - consistent_core + logger.info('Removed: |R| = {}, R = {}'.format( + len(removed_reactions), removed_reactions)) diff --git a/psamm/commands/fba.py b/psamm/commands/fba.py new file mode 100644 index 00000000..2cdacf07 --- /dev/null +++ b/psamm/commands/fba.py @@ -0,0 +1,114 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import logging + +from six import iteritems + +from ..command import SolverCommandMixin, Command, CommandError +from .. import fluxanalysis + +logger = logging.getLogger(__name__) + + +class FluxBalanceCommand(SolverCommandMixin, Command): + """Run flux balance analysis on a metabolic model.""" + + name = 'fba' + title = 'Run flux balance analysis on a metabolic model' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--no-tfba', help='Disable thermodynamic constraints on FBA', + action='store_true') + parser.add_argument( + '--epsilon', type=float, help='Threshold for flux minimization', + default=1e-5) + parser.add_argument('reaction', help='Reaction to maximize', nargs='?') + super(FluxBalanceCommand, cls).init_parser(parser) + + def run(self): + """Run flux analysis command""" + + # Load compound information + compound_name = {} + for compound in self._model.parse_compounds(): + if 'name' in compound.properties: + compound_name[compound.id] = compound.properties['name'] + elif compound.id not in compound_name: + compound_name[compound.id] = compound.id + + if self._args.reaction is not None: + reaction = self._args.reaction + else: + reaction = self._model.get_biomass_reaction() + if reaction is None: + raise CommandError('The biomass reaction was not specified') + + if not self._mm.has_reaction(reaction): + raise CommandError('Specified reaction is not in model: {}'.format( + reaction)) + + if self._args.no_tfba: + result = self.run_fba_minimized(reaction) + else: + result = self.run_tfba(reaction) + + optimum = None + for reaction_id, fba_flux, flux in sorted(result): + rx = self._mm.get_reaction(reaction_id) + print('{}\t{}\t{}\t{}'.format( + reaction_id, fba_flux, flux, + rx.translated_compounds(lambda x: compound_name.get(x, x)))) + # Remember flux of requested reaction + if reaction_id == reaction: + optimum = flux + + logger.info('Maximum flux: {}'.format(optimum)) + + def run_fba_minimized(self, reaction): + """Run normal FBA and flux minimization on model, then print output""" + + solver = self._get_solver() + fba_fluxes = dict(fluxanalysis.flux_balance( + self._mm, reaction, tfba=False, solver=solver)) + optimum = fba_fluxes[reaction] + epsilon = self._args.epsilon + + # Run flux minimization + fmin_fluxes = dict(fluxanalysis.flux_minimization( + self._mm, {reaction: optimum}, solver=solver)) + count = 0 + for reaction_id, flux in iteritems(fmin_fluxes): + if fba_fluxes[reaction_id] - epsilon > flux: + count += 1 + yield reaction_id, fba_fluxes[reaction_id], flux + logger.info('Minimized reactions: {}'.format(count)) + + def run_tfba(self, reaction): + """Run FBA and tFBA on model""" + + solver = self._get_solver(integer=True) + + fba_fluxes = dict(fluxanalysis.flux_balance( + self._mm, reaction, tfba=False, solver=solver)) + fluxes = dict(fluxanalysis.flux_balance( + self._mm, reaction, tfba=True, solver=solver)) + + for reaction_id, flux in iteritems(fluxes): + yield reaction_id, fba_fluxes[reaction_id], flux diff --git a/psamm/commands/fluxcheck.py b/psamm/commands/fluxcheck.py new file mode 100644 index 00000000..ac3b99c1 --- /dev/null +++ b/psamm/commands/fluxcheck.py @@ -0,0 +1,126 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import logging + +from ..command import Command, SolverCommandMixin +from .. import fluxanalysis, fastcore + +logger = logging.getLogger(__name__) + + +class FluxConsistencyCommand(SolverCommandMixin, Command): + """Check that reactions are flux consistent in a model. + + A reaction is flux consistent if there exists any steady-state flux + solution where the flux of the given reaction is non-zero. The + bounds on the exchange reactions can be removed when performing the + consistency check by providing the ``--unrestricted`` option. + """ + + name = 'fluxcheck' + title = 'Check that the model is flux consistent' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--fastcore', help='Enable use of Fastcore algorithm', + action='store_true') + parser.add_argument( + '--no-tfba', + help='Disable thermodynamic constraints on flux check', + action='store_true') + parser.add_argument( + '--epsilon', type=float, help='Flux threshold', + default=1e-5) + parser.add_argument( + '--unrestricted', action='store_true', + help='Remove limits on exchange reactions before checking' + ) + super(FluxConsistencyCommand, cls).init_parser(parser) + + def run(self): + """Run flux consistency check command""" + + # Load compound information + compound_name = {} + for compound in self._model.parse_compounds(): + if 'name' in compound.properties: + compound_name[compound.id] = compound.properties['name'] + elif compound.id not in compound_name: + compound_name[compound.id] = compound.id + + epsilon = self._args.epsilon + + if self._args.unrestricted: + # 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 + + if self._args.fastcore: + solver = self._get_solver() + inconsistent = set(fastcore.fastcc( + self._mm, epsilon, solver=solver)) + else: + enable_tfba = not self._args.no_tfba + if enable_tfba: + solver = self._get_solver(integer=True) + else: + solver = self._get_solver() + inconsistent = set( + fluxanalysis.consistency_check( + self._mm, self._mm.reactions, epsilon, + tfba=enable_tfba, solver=solver)) + + # Count the number of reactions that are fixed at zero. While these + # reactions are still inconsistent, they are inconsistent because they + # have been explicitly disabled. + disabled_exchange = 0 + disabled_internal = 0 + + count_exchange = 0 + total_exchange = 0 + + count_internal = 0 + total_internal = 0 + + # Print result + for reaction in sorted(self._mm.reactions): + disabled = self._mm.limits[reaction].bounds == (0, 0) + + if self._mm.is_exchange(reaction): + total_exchange += 1 + count_exchange += int(reaction in inconsistent) + disabled_exchange += int(disabled) + else: + total_internal += 1 + count_internal += int(reaction in inconsistent) + disabled_internal += int(disabled) + + if reaction in inconsistent: + rx = self._mm.get_reaction(reaction) + rxt = rx.translated_compounds( + lambda x: compound_name.get(x, x)) + print('{}\t{}'.format(reaction, rxt)) + + logger.info('Model has {}/{} inconsistent internal reactions' + ' ({} disabled by user)'.format( + count_internal, total_internal, disabled_internal)) + logger.info('Model has {}/{} inconsistent exchange reactions' + ' ({} disabled by user)'.format( + count_exchange, total_exchange, disabled_exchange)) diff --git a/psamm/commands/fluxcoupling.py b/psamm/commands/fluxcoupling.py new file mode 100644 index 00000000..6a18c2a9 --- /dev/null +++ b/psamm/commands/fluxcoupling.py @@ -0,0 +1,138 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import logging + +from ..command import Command, SolverCommandMixin, CommandError +from .. import fluxanalysis, fluxcoupling + +logger = logging.getLogger(__name__) + + +class FluxCouplingCommand(SolverCommandMixin, Command): + """Find flux coupled reactions in the model. + + This identifies any reaction pairs where the flux of one reaction + constrains the flux of another reaction. The reactions can be coupled in + three distinct ways depending on the ratio between the reaction fluxes. + The reactions can be fully coupled (the ratio is static and non-zero); + partially coupled (the ratio is bounded and non-zero); or directionally + coupled (the ratio is non-zero). + """ + + name = 'fluxcoupling' + title = 'Find flux coupled reactions in model' + + def run(self): + solver = self._get_solver() + + max_reaction = self._model.get_biomass_reaction() + if max_reaction is None: + raise CommandError('The biomass reaction was not specified') + + fba_fluxes = dict(fluxanalysis.flux_balance( + self._mm, max_reaction, tfba=False, solver=solver)) + optimum = fba_fluxes[max_reaction] + + self._fcp = fluxcoupling.FluxCouplingProblem( + self._mm, {max_reaction: 0.999 * optimum}, solver) + + self._coupled = {} + self._groups = [] + + reactions = sorted(self._mm.reactions) + for i, reaction1 in enumerate(reactions): + if reaction1 in self._coupled: + continue + + for reaction2 in reactions[i+1:]: + if (reaction2 in self._coupled and + (self._coupled[reaction2] == + self._coupled.get(reaction1))): + continue + + self._check_reactions(reaction1, reaction2) + + logger.info('Coupled groups:') + for i, group in enumerate(self._groups): + if group is not None: + logger.info('{}: {}'.format(i, ', '.join(sorted(group)))) + + def _check_reactions(self, reaction1, reaction2): + logger.debug('Solving for {}, {}'.format(reaction1, reaction2)) + lower, upper = self._fcp.solve(reaction1, reaction2) + + logger.debug('Result: {}, {}'.format(lower, upper)) + + coupling = fluxcoupling.classify_coupling((lower, upper)) + if coupling in (fluxcoupling.CouplingClass.DirectionalForward, + fluxcoupling.CouplingClass.DirectionalReverse): + text = 'Directional, v1 / v2 in [{}, {}]'.format(lower, upper) + if (coupling == fluxcoupling.CouplingClass.DirectionalReverse and + not self._mm.is_reversible(reaction1) and lower == 0.0): + return + elif coupling == fluxcoupling.CouplingClass.Full: + text = 'Full: v1 / v2 = {}'.format(lower) + self._couple_reactions(reaction1, reaction2) + elif coupling == fluxcoupling.CouplingClass.Partial: + text = 'Partial: v1 / v2 in [{}; {}]'.format(lower, upper) + self._couple_reactions(reaction1, reaction2) + else: + return + + print('{}\t{}\t{}\t{}\t{}'.format( + reaction1, reaction2, lower, upper, text)) + + def _couple_reactions(self, reaction1, reaction2): + logger.debug('Couple {} and {}'.format(reaction1, reaction2)) + + if reaction1 in self._coupled and reaction2 in self._coupled: + if self._coupled[reaction1] == self._coupled[reaction2]: + return + logger.debug('Merge groups {}, {}'.format( + self._coupled[reaction1], self._coupled[reaction2])) + group_index = len(self._groups) + group = (self._groups[self._coupled[reaction1]] | + self._groups[self._coupled[reaction2]]) + logger.debug('New group is {}: {}'.format( + group_index, sorted(group))) + self._groups.append(group) + self._groups[self._coupled[reaction1]] = None + self._groups[self._coupled[reaction2]] = None + for reaction in group: + self._coupled[reaction] = group_index + elif reaction1 in self._coupled: + group_index = self._coupled[reaction1] + group = self._groups[group_index] + logger.debug('Put {} into existing group {}: {}'.format( + reaction2, self._coupled[reaction1], group)) + self._coupled[reaction2] = group_index + group.add(reaction2) + elif reaction2 in self._coupled: + group_index = self._coupled[reaction2] + group = self._groups[group_index] + logger.debug('Put {} into existing group {}: {}'.format( + reaction1, self._coupled[reaction2], group)) + self._coupled[reaction1] = group_index + group.add(reaction1) + else: + group = set([reaction1, reaction2]) + group_index = len(self._groups) + logger.debug('Creating new group {}'.format(group_index)) + self._coupled[reaction1] = group_index + self._coupled[reaction2] = group_index + self._groups.append(group) diff --git a/psamm/commands/formulacheck.py b/psamm/commands/formulacheck.py new file mode 100644 index 00000000..876efde3 --- /dev/null +++ b/psamm/commands/formulacheck.py @@ -0,0 +1,98 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import operator +import logging + +from ..command import Command +from ..formula import Formula + +logger = logging.getLogger(__name__) + + +class FormulaBalanceCommand(Command): + """Check whether reactions in a model are elementally balanced. + + Balanced reactions are those reactions where the number of elements + (atoms) is consistent on the left and right side of the reaction equation. + Reactions that are not balanced will be printed out. + """ + + name = 'formulacheck' + title = 'Check formula balance on a model or database' + + def run(self): + """Run formula balance command""" + + # Mapping from compound id to formula + compound_formula = {} + for compound in self._model.parse_compounds(): + # Only cpd11632 (Photon) is allowed to have an empty formula. + if compound.formula is None: + if compound.id == 'cpd11632': + compound_formula[compound.id] = Formula() + else: + try: + f = Formula.parse(compound.formula).flattened() + compound_formula[compound.id] = f + except ValueError as e: + logger.warning( + 'Error parsing {}: {}'.format(compound.formula, e)) + + # Create a set of known mass-inconsistent reactions + exchange = set() + for reaction_id in self._mm.reactions: + if self._mm.is_exchange(reaction_id): + exchange.add(reaction_id) + + def multiply_formula(compound_list): + for compound, count in compound_list: + yield count * compound_formula.get(compound.name, Formula()) + + count = 0 + unbalanced = 0 + unchecked = 0 + for reaction in self._mm.reactions: + if reaction in exchange: + continue + + count += 1 + rx = self._mm.get_reaction(reaction) + + # Skip reaction if any compounds have undefined formula + for compound, _ in rx.compounds: + if compound.name not in compound_formula: + unchecked += 1 + break + else: + left_form = reduce( + operator.or_, multiply_formula(rx.left), Formula()) + right_form = reduce( + operator.or_, multiply_formula(rx.right), Formula()) + + if right_form != left_form: + unbalanced += 1 + right_missing, left_missing = Formula.balance( + right_form, left_form) + + print('{}\t{}\t{}\t{}\t{}'.format( + reaction, left_form, right_form, + left_missing, right_missing)) + + logger.info('Unbalanced reactions: {}/{}'.format(unbalanced, count)) + logger.info('Unchecked reactions due to missing formula: {}/{}'.format( + unchecked, count)) diff --git a/psamm/commands/fva.py b/psamm/commands/fva.py new file mode 100644 index 00000000..d4080e04 --- /dev/null +++ b/psamm/commands/fva.py @@ -0,0 +1,75 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +from ..command import Command, SolverCommandMixin, CommandError +from .. import fluxanalysis + + +class FluxVariabilityCommand(SolverCommandMixin, Command): + """Run flux variablity analysis on a metabolic model.""" + + name = 'fva' + title = 'Run flux variability analysis on a metabolic model' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--no-tfba', help='Disable thermodynamic constraints on FVA', + action='store_true') + parser.add_argument('reaction', help='Reaction to maximize', nargs='?') + super(FluxVariabilityCommand, cls).init_parser(parser) + + def run(self): + """Run flux variability command""" + + # Load compound information + compound_name = {} + for compound in self._model.parse_compounds(): + if 'name' in compound.properties: + compound_name[compound.id] = compound.properties['name'] + elif compound.id not in compound_name: + compound_name[compound.id] = compound.id + + if self._args.reaction is not None: + reaction = self._args.reaction + else: + reaction = self._model.get_biomass_reaction() + if reaction is None: + raise CommandError('The biomass reaction was not specified') + + if not self._mm.has_reaction(reaction): + raise CommandError('Specified reaction is not in model: {}'.format( + reaction)) + + enable_tfba = not self._args.no_tfba + if enable_tfba: + solver = self._get_solver(integer=True) + else: + solver = self._get_solver() + + fba_fluxes = dict(fluxanalysis.flux_balance( + self._mm, reaction, tfba=False, solver=solver)) + optimum = fba_fluxes[reaction] + + flux_bounds = fluxanalysis.flux_variability( + self._mm, sorted(self._mm.reactions), {reaction: optimum}, + tfba=enable_tfba, solver=solver) + for reaction_id, bounds in flux_bounds: + rx = self._mm.get_reaction(reaction_id) + rxt = rx.translated_compounds(lambda x: compound_name.get(x, x)) + print('{}\t{}\t{}\t{}'.format( + reaction_id, bounds[0], bounds[1], rxt)) diff --git a/psamm/commands/gapfill.py b/psamm/commands/gapfill.py new file mode 100644 index 00000000..e27278e4 --- /dev/null +++ b/psamm/commands/gapfill.py @@ -0,0 +1,80 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import logging + +from ..command import Command, SolverCommandMixin +from ..gapfill import gapfind, gapfill + +logger = logging.getLogger(__name__) + + +class GapFillCommand(SolverCommandMixin, Command): + """Run GapFind and GapFill on a metabolic model.""" + + name = 'gapfill' + title = 'Run GapFind and GapFill on a metabolic model' + + def run(self): + """Run GapFill command""" + + # Load compound information + compound_name = {} + for compound in self._model.parse_compounds(): + if 'name' in compound.properties: + compound_name[compound.id] = compound.properties['name'] + elif compound.id not in compound_name: + compound_name[compound.id] = compound.id + + model_compartments = set(self._mm.compartments) + + solver = self._get_solver(integer=True) + + # Run GapFind on model + logger.info('Searching for blocked compounds') + blocked = set(compound for compound in gapfind(self._mm, solver=solver) + if compound.compartment is not 'e') + if len(blocked) > 0: + logger.info('Blocked compounds') + for compound in blocked: + print(compound.translate(lambda x: compound_name.get(x, x))) + + if len(blocked) > 0: + # Add exchange and transport reactions to database + model_complete = self._mm.copy() + logger.info('Adding database, exchange and transport reactions') + model_complete.add_all_database_reactions(model_compartments) + model_complete.add_all_exchange_reactions() + model_complete.add_all_transport_reactions() + + logger.info('Searching for reactions to fill gaps') + added_reactions, reversed_reactions = gapfill( + model_complete, self._mm.reactions, blocked, solver=solver) + + for rxnid in added_reactions: + rx = model_complete.get_reaction(rxnid) + rxt = rx.translated_compounds( + lambda x: compound_name.get(x, x)) + print('{}\t{}\t{}'.format(rxnid, 'Add', rxt)) + + for rxnid in reversed_reactions: + rx = model_complete.get_reaction(rxnid) + rxt = rx.translated_compounds( + lambda x: compound_name.get(x, x)) + print('{}\t{}\t{}'.format(rxnid, 'Reverse', rxt)) + else: + logger.info('No blocked compounds found') diff --git a/psamm/commands/masscheck.py b/psamm/commands/masscheck.py new file mode 100644 index 00000000..7c8eeec4 --- /dev/null +++ b/psamm/commands/masscheck.py @@ -0,0 +1,115 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import logging + +from ..command import Command, SolverCommandMixin +from .. import massconsistency + +logger = logging.getLogger(__name__) + + +class MassConsistencyCommand(SolverCommandMixin, Command): + """Check whether a model is mass consistent.""" + + name = 'masscheck' + title = 'Run mass consistency check on a database' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--exclude', metavar='reaction', action='append', type=str, + default=[], help='Exclude reaction from mass consistency') + parser.add_argument( + '--epsilon', type=float, help='Mass threshold', + default=1e-5) + parser.add_argument( + '--checked', metavar='reaction', action='append', type=str, + default=[], help='Mark reaction as already checked (no residual)') + super(MassConsistencyCommand, cls).init_parser(parser) + + def run(self): + # Load compound information + compound_name = {} + zeromass = set() + for compound in self._model.parse_compounds(): + if 'name' in compound.properties: + compound_name[compound.id] = compound.properties['name'] + elif compound.id not in compound_name: + compound_name[compound.id] = compound.id + + if compound.properties.get('zeromass', False): + zeromass.add(compound.id) + + # Create a set of known mass-inconsistent reactions + exchange = set() + for reaction_id in self._mm.reactions: + if self._mm.is_exchange(reaction_id): + exchange.add(reaction_id) + + # Other reactions to exclude from consistency check + exclude = set(self._args.exclude) + + # Add biomass reaction to be excluded + biomass_reaction = self._model.get_biomass_reaction() + if biomass_reaction is not None: + exclude.add(biomass_reaction) + + # Create set of checked reactions + checked = set(self._args.checked) + + solver = self._get_solver() + + known_inconsistent = exclude | exchange + + logger.info('Mass consistency on database') + epsilon = self._args.epsilon + compound_iter = massconsistency.check_compound_consistency( + self._mm, solver, known_inconsistent, zeromass) + + good = 0 + total = 0 + for compound, mass in sorted(compound_iter, key=lambda x: (x[1], x[0]), + reverse=True): + if mass >= 1-epsilon or compound.name in zeromass: + good += 1 + total += 1 + print('{}: {}'.format(compound.translate( + lambda x: compound_name.get(x, x)), mass)) + logger.info('Consistent compounds: {}/{}'.format(good, total)) + + logger.info('Is consistent? {}'.format( + massconsistency.is_consistent( + self._mm, solver, known_inconsistent, zeromass))) + + good = 0 + total = 0 + reaction_iter, compound_iter = ( + massconsistency.check_reaction_consistency( + self._mm, solver=solver, exchange=known_inconsistent, + checked=checked, zeromass=zeromass)) + for reaction_id, residual in sorted( + reaction_iter, key=lambda x: abs(x[1]), reverse=True): + total += 1 + if abs(residual) >= epsilon: + reaction = self._mm.get_reaction(reaction_id) + rxt = reaction.translated_compounds( + lambda x: compound_name.get(x, x)) + print('{}\t{}\t{}'.format(reaction_id, residual, rxt)) + else: + good += 1 + logger.info('Consistent reactions: {}/{}'.format(good, total)) diff --git a/psamm/commands/randomsparse.py b/psamm/commands/randomsparse.py new file mode 100644 index 00000000..55d06aba --- /dev/null +++ b/psamm/commands/randomsparse.py @@ -0,0 +1,136 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen +# Copyright 2015 Keith Dufault-Thompson + +import random +import logging + +from ..command import Command, SolverCommandMixin, CommandError +from .. import fluxanalysis, util + +logger = logging.getLogger(__name__) + + +class RandomSparseNetworkCommand(SolverCommandMixin, Command): + """Find random minimal network of model. + + Given a reaction to optimize and a threshold, delete reactions randomly + until the flux of the reaction to optimize falls under the threshold. + Keep deleting reactions until no more reactions can be deleted. By default + this uses standard FBA (not tFBA). Since the internal fluxes are irrelevant + the FBA and tFBA are equivalent for this purpose. + + The threshold can be specified as an absolute flux (e.g. '1.23') or a + relative flux of the full model flux (e.g. '40.5%'). + """ + + name = 'randomsparse' + title = 'Generate a random sparse network of model reactions' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--tfba', help='Enable thermodynamic constraints on FBA', + action='store_true', default=False) + parser.add_argument( + '--reaction', help='Reaction to maximize') + parser.add_argument( + 'threshold', help='Threshold of max reaction flux', + type=util.MaybeRelative) + parser.add_argument( + '--exchange', help='Only analyze the exchange reaction in model', + action='store_true') + super(RandomSparseNetworkCommand, cls).init_parser(parser) + + def run(self): + if self._args.reaction is not None: + reaction = self._args.reaction + else: + reaction = self._model.get_biomass_reaction() + if reaction is None: + raise CommandError('The biomass reaction was not specified') + + if not self._mm.has_reaction(reaction): + raise CommandError('Specified reaction is not in model: {}'.format( + reaction)) + + if not self._args.tfba: + fb_problem = fluxanalysis.FluxBalanceProblem + solver = self._get_solver() + else: + fb_problem = fluxanalysis.FluxBalanceTDProblem + solver = self._get_solver(integer=True) + + p = fb_problem(self._mm, solver) + + threshold = self._args.threshold + if threshold.relative: + p.solve(reaction) + threshold.reference = p.get_flux(reaction) + + flux_threshold = float(threshold) + + logger.info('Flux threshold for {} is {}'.format( + reaction, flux_threshold)) + + essential = {reaction} + deleted = set() + if self._args.exchange: + reactions = set() + for reaction_id in self._mm.reactions: + if self._mm.is_exchange(reaction_id): + reactions.add(reaction_id) + else: + reactions = set(self._mm.reactions) + + test_set = reactions - essential + + while len(test_set) > 0: + testing_reaction = random.sample(test_set, 1)[0] + test_set.remove(testing_reaction) + + flux_var = p.get_flux_var(testing_reaction) + c, = p.prob.add_linear_constraints(flux_var == 0) + + logger.info('Trying FBA without reaction {}...'.format( + testing_reaction)) + + try: + p.solve(reaction) + except fluxanalysis.FluxBalanceError: + logger.info( + 'FBA is infeasible, marking {} as essential'.format( + testing_reaction)) + c.delete() + essential.add(testing_reaction) + continue + + logger.debug('Reaction {} has flux {}'.format( + reaction, p.get_flux(reaction))) + + if p.get_flux(reaction) < flux_threshold: + c.delete() + essential.add(testing_reaction) + logger.info('Reaction {} was essential'.format( + testing_reaction)) + else: + deleted.add(testing_reaction) + logger.info('Reaction {} was deleted'.format(testing_reaction)) + + for reaction_id in sorted(reactions): + value = 0 if reaction_id in deleted else 1 + print('{}\t{}'.format(reaction_id, value)) diff --git a/psamm/commands/robustness.py b/psamm/commands/robustness.py new file mode 100644 index 00000000..2f957605 --- /dev/null +++ b/psamm/commands/robustness.py @@ -0,0 +1,125 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +from ..command import Command, SolverCommandMixin, CommandError +from .. import fluxanalysis + + +class RobustnessCommand(SolverCommandMixin, Command): + """Run robustness analysis on metabolic model. + + Given a reaction to maximize and a reaction to vary, + the robustness analysis will run FBA while fixing the + reaction to vary at each iteration. The reaction will + be fixed at the specified number of steps between the + minimum and maximum flux value specified in the model. + """ + + name = 'robustness' + title = 'Run robustness analysis on a metabolic model' + + @classmethod + def init_parser(cls, parser): + parser.add_argument( + '--steps', metavar='N', type=int, default=10, + help='Number of flux value steps for varying reaction') + parser.add_argument( + '--minimum', metavar='V', type=float, + help='Minumum flux value of varying reacton') + parser.add_argument( + '--maximum', metavar='V', type=float, + help='Maximum flux value of varying reacton') + parser.add_argument( + '--no-tfba', help='Disable thermodynamic constraints on FBA', + action='store_true') + parser.add_argument( + '--reaction', help='Reaction to maximize', nargs='?') + parser.add_argument('varying', help='Reaction to vary') + super(RobustnessCommand, cls).init_parser(parser) + + def run(self): + """Run flux analysis command""" + + def run_fba_fmin(model, reaction): + fba = fluxanalysis.FluxBalanceProblem(model, solver=solver) + fba.solve(reaction) + optimum = fba.get_flux(reaction) + return fluxanalysis.flux_minimization( + model, {reaction: optimum}, solver=solver) + + def run_tfba(model, reaction): + return fluxanalysis.flux_balance( + model, reaction, tfba=True, solver=solver) + + if self._args.no_tfba: + run_fba = run_fba_fmin + solver = self._get_solver() + else: + run_fba = run_tfba + solver = self._get_solver(integer=True) + + # Load compound information + compound_name = {} + for compound in self._model.parse_compounds(): + if 'name' in compound.properties: + compound_name[compound.id] = compound.properties['name'] + elif compound.id not in compound_name: + compound_name[compound.id] = compound.id + + if self._args.reaction is not None: + reaction = self._args.reaction + else: + reaction = self._model.get_biomass_reaction() + if reaction is None: + raise CommandError('The biomass reaction was not specified') + + if not self._mm.has_reaction(reaction): + raise CommandError('Specified reaction is not in model: {}'.format( + reaction)) + + varying_reaction = self._args.varying + if not self._mm.has_reaction(varying_reaction): + raise CommandError('Specified reaction is not in model: {}'.format( + varying_reaction)) + + steps = self._args.steps + if steps <= 0: + raise CommandError('Invalid number of steps: {}\n'.format(steps)) + + # Run FBA on model at different fixed flux values + flux_min = self._mm.limits[varying_reaction].lower + flux_max = self._mm.limits[varying_reaction].upper + if self._args.minimum is not None: + flux_min = self._args.minimum + if self._args.maximum is not None: + flux_max = self._args.maximum + + if flux_min > flux_max: + raise CommandError('Invalid flux range: {}, {}\n'.format( + flux_min, flux_max)) + + for i in xrange(steps): + fixed_flux = flux_min + i*(flux_max - flux_min)/float(steps-1) + test_model = self._mm.copy() + test_model.limits[varying_reaction].bounds = fixed_flux, fixed_flux + + try: + for other_reaction, flux in run_fba(test_model, reaction): + print('{}\t{}\t{}'.format( + other_reaction, fixed_flux, flux)) + except fluxanalysis.FluxBalanceError: + pass diff --git a/psamm/commands/sbmlexport.py b/psamm/commands/sbmlexport.py new file mode 100644 index 00000000..a72df78a --- /dev/null +++ b/psamm/commands/sbmlexport.py @@ -0,0 +1,33 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import sys + +from ..datasource import sbml +from ..command import Command + + +class SBMLExport(Command): + """Export model as SBML file.""" + + name = 'sbmlexport' + title = 'Export model as SBML file' + + def run(self): + writer = sbml.SBMLWriter() + writer.write_model( + sys.stdout, self._mm, self._model.parse_compounds()) diff --git a/psamm/commands/search.py b/psamm/commands/search.py new file mode 100644 index 00000000..6d95d820 --- /dev/null +++ b/psamm/commands/search.py @@ -0,0 +1,145 @@ +# This file is part of PSAMM. +# +# PSAMM is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# PSAMM is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with PSAMM. If not, see . +# +# Copyright 2014-2015 Jon Lund Steffensen + +import re + +from ..command import Command +from ..reaction import Compound + + +class SearchCommand(Command): + """Search for reactions and compounds in a model.""" + + name = 'search' + title = 'Search the database of reactions or compounds' + + @classmethod + def init_parser(cls, parser): + """Initialize argument parser""" + subparsers = parser.add_subparsers(title='Search domain') + + # Compound subcommand + parser_compound = subparsers.add_parser( + 'compound', help='Search in compounds') + parser_compound.set_defaults(which='compound') + parser_compound.add_argument( + '--id', '-i', dest='id', metavar='id', type=str, + default=[], action='append', help='Compound ID') + parser_compound.add_argument( + '--name', '-n', dest='name', metavar='name', type=str, + default=[], action='append', help='Name of compound') + + # Reaction subcommand + parser_reaction = subparsers.add_parser( + 'reaction', help='Search in reactions') + parser_reaction.set_defaults(which='reaction') + parser_reaction.add_argument( + '--id', '-i', dest='id', metavar='id', type=str, + default=[], action='append', help='Reaction ID') + parser_reaction.add_argument( + '--compound', '-c', dest='compound', metavar='compound', type=str, + default=[], action='append', + help='Comma-separated list of compound IDs') + + def run(self): + """Run search command.""" + + def parse_compound(s): + """Parse a compound specification with optional compartment""" + m = re.match(r'([^\[]+)\[(\w+)\]', s) + if m is not None: + return Compound(m.group(1), compartment=m.group(2)) + return Compound(s) + + def filter_search_term(s): + return re.sub(r'[^a-z0-9]+', '', s.lower()) + + which_command = self._args.which + if which_command == 'compound': + selected_compounds = set() + + for compound in self._model.parse_compounds(): + if len(self._args.id) > 0: + if any(c == compound.id for c in self._args.id): + selected_compounds.add(compound) + continue + + if len(self._args.name) > 0: + names = set() + if 'name' in compound.properties: + names.add(compound.properties['name']) + names.update(compound.properties.get('names', [])) + names = set(filter_search_term(n) for n in names) + if any(filter_search_term(n) in names + for n in self._args.name): + selected_compounds.add(compound) + continue + + # Show results + for compound in selected_compounds: + print('ID: {}'.format(compound.id)) + if 'name' in compound.properties: + print('Name: {}'.format(compound.properties['name'])) + if 'names' in compound.properties: + print('Additional names: {}'.format( + ', '.join(compound.properties['names']))) + if 'formula' in compound.properties: + print('Formula: {}'.format(compound.properties['formula'])) + if compound.filemark is not None: + print('Parsed from: {}'.format(compound.filemark)) + print() + elif which_command == 'reaction': + selected_reactions = set() + + # Prepare translation table from compound id to name + compound_name = {} + for compound in self._model.parse_compounds(): + if 'name' in compound.properties: + compound_name[compound.id] = compound.properties['name'] + + # Prepare sets of matched compounds + search_compounds = [] + for compound_list in self._args.compound: + compound_set = set() + for compound_spec in compound_list.split(','): + compound_set.add(parse_compound(compound_spec.strip())) + search_compounds.append(compound_set) + + for reaction in self._model.parse_reactions(): + if len(self._args.id) > 0: + if any(r == reaction.id for r in self._args.id): + selected_reactions.add(reaction) + continue + + if len(search_compounds) > 0: + compounds = set(c for c, _ in reaction.equation.compounds) + if any(c.issubset(compounds) for c in search_compounds): + selected_reactions.add(reaction) + continue + + # Show results + for reaction in selected_reactions: + print('ID: {}'.format(reaction.id)) + print('Reaction (IDs): {}'.format( + reaction.equation)) + translated_equation = reaction.equation.translated_compounds( + lambda x: compound_name.get(x, x)) + if reaction.equation != translated_equation: + print('Reaction (names): {}'.format(translated_equation)) + if reaction.filemark is not None: + print('Parsed from: {}'.format(reaction.filemark)) + print() diff --git a/setup.py b/setup.py index 90dbf8c2..bebff9f8 100755 --- a/setup.py +++ b/setup.py @@ -46,6 +46,22 @@ 'console_scripts': [ 'psamm-model = psamm.command:main', 'psamm-list-lpsolvers = psamm.lpsolver.generic:list_solvers' + ], + 'psamm.commands': [ + 'chargecheck = psamm.commands.chargecheck:ChargeBalanceCommand', + 'console = psamm.commands.console:ConsoleCommand', + 'fastgapfill = psamm.commands.fastgapfill:FastGapFillCommand', + 'fba = psamm.commands.fba:FluxBalanceCommand', + 'fluxcheck = psamm.commands.fluxcheck:FluxConsistencyCommand', + 'fluxcoupling = psamm.commands.fluxcoupling:FluxCouplingCommand', + 'formulacheck = psamm.commands.formulacheck:FormulaBalanceCommand', + 'fva = psamm.commands.fva:FluxVariabilityCommand', + 'gapfill = psamm.commands.gapfill:GapFillCommand', + 'masscheck = psamm.commands.masscheck:MassConsistencyCommand', + 'randomsparse = psamm.commands.randomsparse:RandomSparseNetworkCommand', + 'robustness = psamm.commands.robustness:RobustnessCommand', + 'sbmlexport = psamm.commands.sbmlexport:SBMLExport', + 'search = psamm.commands.search:SearchCommand' ] }, From 93c04ec45d27405288958721bd3bea5374073755 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 15:05:30 -0400 Subject: [PATCH 31/46] command: Obtain command name and help from docstring Parse the docstring on the Command subclasses to obtain the name and help parameters provided to add_parser(). --- psamm/command.py | 5 +++-- psamm/commands/chargecheck.py | 5 +---- psamm/commands/console.py | 5 +---- psamm/commands/fastgapfill.py | 5 +---- psamm/commands/fba.py | 5 +---- psamm/commands/fluxcheck.py | 3 --- psamm/commands/fluxcoupling.py | 3 --- psamm/commands/formulacheck.py | 5 +---- psamm/commands/fva.py | 5 +---- psamm/commands/gapfill.py | 5 +---- psamm/commands/masscheck.py | 5 +---- psamm/commands/randomsparse.py | 5 +---- psamm/commands/robustness.py | 5 +---- psamm/commands/sbmlexport.py | 3 --- psamm/commands/search.py | 5 +---- 15 files changed, 14 insertions(+), 55 deletions(-) diff --git a/psamm/command.py b/psamm/command.py index 86a94b09..0c37dc9c 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -177,10 +177,11 @@ def main(command_class=None): canonical)) # Create parsers for subcommands - subparsers = parser.add_subparsers(title='Command') + subparsers = parser.add_subparsers(title='Commands', metavar='command') for name, command_class in sorted(iteritems(commands)): + title, _, _ = command_class.__doc__.partition('\n\n') subparser = subparsers.add_parser( - name, help=command_class.title, + name, help=title.rstrip('.'), description=command_class.__doc__) subparser.set_defaults(command=command_class) command_class.init_parser(subparser) diff --git a/psamm/commands/chargecheck.py b/psamm/commands/chargecheck.py index 6ddbe24a..7255cb68 100644 --- a/psamm/commands/chargecheck.py +++ b/psamm/commands/chargecheck.py @@ -24,16 +24,13 @@ class ChargeBalanceCommand(Command): - """Check whether compound charge in a given database or model is balanced. + """Check whether compound charge is balanced. Balanced reactions are those reactions where the total charge is consistent on the left and right side of the reaction equation. Reactions that are not balanced will be printed out. """ - name = 'chargecheck' - title = 'Check charge balance on a model or database' - def run(self): """Run charge balance command""" diff --git a/psamm/commands/console.py b/psamm/commands/console.py index 7a04e3e1..37f40c94 100644 --- a/psamm/commands/console.py +++ b/psamm/commands/console.py @@ -19,10 +19,7 @@ class ConsoleCommand(Command): - """Start an interactive Python console with the given model loaded.""" - - name = 'console' - title = 'Start Python console with metabolic model loaded' + """Start an interactive Python console with the model loaded.""" @classmethod def init_parser(cls, parser): diff --git a/psamm/commands/fastgapfill.py b/psamm/commands/fastgapfill.py index d33055dd..6c5f2ca0 100644 --- a/psamm/commands/fastgapfill.py +++ b/psamm/commands/fastgapfill.py @@ -25,10 +25,7 @@ class FastGapFillCommand(SolverCommandMixin, Command): - """Run FastGapFill algorithm on a metabolic model.""" - - name = 'fastgapfill' - title = 'Run FastGapFill on a metabolic model' + """Run the FastGapFill gap-filling algorithm on model.""" @classmethod def init_parser(cls, parser): diff --git a/psamm/commands/fba.py b/psamm/commands/fba.py index 2cdacf07..d9df1f7c 100644 --- a/psamm/commands/fba.py +++ b/psamm/commands/fba.py @@ -26,10 +26,7 @@ class FluxBalanceCommand(SolverCommandMixin, Command): - """Run flux balance analysis on a metabolic model.""" - - name = 'fba' - title = 'Run flux balance analysis on a metabolic model' + """Run flux balance analysis on the model.""" @classmethod def init_parser(cls, parser): diff --git a/psamm/commands/fluxcheck.py b/psamm/commands/fluxcheck.py index ac3b99c1..08ca6a08 100644 --- a/psamm/commands/fluxcheck.py +++ b/psamm/commands/fluxcheck.py @@ -32,9 +32,6 @@ class FluxConsistencyCommand(SolverCommandMixin, Command): consistency check by providing the ``--unrestricted`` option. """ - name = 'fluxcheck' - title = 'Check that the model is flux consistent' - @classmethod def init_parser(cls, parser): parser.add_argument( diff --git a/psamm/commands/fluxcoupling.py b/psamm/commands/fluxcoupling.py index 6a18c2a9..24922a2c 100644 --- a/psamm/commands/fluxcoupling.py +++ b/psamm/commands/fluxcoupling.py @@ -34,9 +34,6 @@ class FluxCouplingCommand(SolverCommandMixin, Command): coupled (the ratio is non-zero). """ - name = 'fluxcoupling' - title = 'Find flux coupled reactions in model' - def run(self): solver = self._get_solver() diff --git a/psamm/commands/formulacheck.py b/psamm/commands/formulacheck.py index 876efde3..cde99a73 100644 --- a/psamm/commands/formulacheck.py +++ b/psamm/commands/formulacheck.py @@ -25,16 +25,13 @@ class FormulaBalanceCommand(Command): - """Check whether reactions in a model are elementally balanced. + """Check whether reactions in the model are elementally balanced. Balanced reactions are those reactions where the number of elements (atoms) is consistent on the left and right side of the reaction equation. Reactions that are not balanced will be printed out. """ - name = 'formulacheck' - title = 'Check formula balance on a model or database' - def run(self): """Run formula balance command""" diff --git a/psamm/commands/fva.py b/psamm/commands/fva.py index d4080e04..d938e11e 100644 --- a/psamm/commands/fva.py +++ b/psamm/commands/fva.py @@ -20,10 +20,7 @@ class FluxVariabilityCommand(SolverCommandMixin, Command): - """Run flux variablity analysis on a metabolic model.""" - - name = 'fva' - title = 'Run flux variability analysis on a metabolic model' + """Run flux variablity analysis on the model.""" @classmethod def init_parser(cls, parser): diff --git a/psamm/commands/gapfill.py b/psamm/commands/gapfill.py index e27278e4..6450e75d 100644 --- a/psamm/commands/gapfill.py +++ b/psamm/commands/gapfill.py @@ -24,10 +24,7 @@ class GapFillCommand(SolverCommandMixin, Command): - """Run GapFind and GapFill on a metabolic model.""" - - name = 'gapfill' - title = 'Run GapFind and GapFill on a metabolic model' + """Run the GapFind and GapFill algorithms on the model.""" def run(self): """Run GapFill command""" diff --git a/psamm/commands/masscheck.py b/psamm/commands/masscheck.py index 7c8eeec4..f14c7f0f 100644 --- a/psamm/commands/masscheck.py +++ b/psamm/commands/masscheck.py @@ -24,10 +24,7 @@ class MassConsistencyCommand(SolverCommandMixin, Command): - """Check whether a model is mass consistent.""" - - name = 'masscheck' - title = 'Run mass consistency check on a database' + """Check whether the model is mass consistent.""" @classmethod def init_parser(cls, parser): diff --git a/psamm/commands/randomsparse.py b/psamm/commands/randomsparse.py index 55d06aba..4888d639 100644 --- a/psamm/commands/randomsparse.py +++ b/psamm/commands/randomsparse.py @@ -26,7 +26,7 @@ class RandomSparseNetworkCommand(SolverCommandMixin, Command): - """Find random minimal network of model. + """Find a random minimal network of model reactions. Given a reaction to optimize and a threshold, delete reactions randomly until the flux of the reaction to optimize falls under the threshold. @@ -38,9 +38,6 @@ class RandomSparseNetworkCommand(SolverCommandMixin, Command): relative flux of the full model flux (e.g. '40.5%'). """ - name = 'randomsparse' - title = 'Generate a random sparse network of model reactions' - @classmethod def init_parser(cls, parser): parser.add_argument( diff --git a/psamm/commands/robustness.py b/psamm/commands/robustness.py index 2f957605..768450b7 100644 --- a/psamm/commands/robustness.py +++ b/psamm/commands/robustness.py @@ -20,7 +20,7 @@ class RobustnessCommand(SolverCommandMixin, Command): - """Run robustness analysis on metabolic model. + """Run robustness analysis on the model. Given a reaction to maximize and a reaction to vary, the robustness analysis will run FBA while fixing the @@ -29,9 +29,6 @@ class RobustnessCommand(SolverCommandMixin, Command): minimum and maximum flux value specified in the model. """ - name = 'robustness' - title = 'Run robustness analysis on a metabolic model' - @classmethod def init_parser(cls, parser): parser.add_argument( diff --git a/psamm/commands/sbmlexport.py b/psamm/commands/sbmlexport.py index a72df78a..f7f7f697 100644 --- a/psamm/commands/sbmlexport.py +++ b/psamm/commands/sbmlexport.py @@ -24,9 +24,6 @@ class SBMLExport(Command): """Export model as SBML file.""" - name = 'sbmlexport' - title = 'Export model as SBML file' - def run(self): writer = sbml.SBMLWriter() writer.write_model( diff --git a/psamm/commands/search.py b/psamm/commands/search.py index 6d95d820..9a531bcd 100644 --- a/psamm/commands/search.py +++ b/psamm/commands/search.py @@ -22,10 +22,7 @@ class SearchCommand(Command): - """Search for reactions and compounds in a model.""" - - name = 'search' - title = 'Search the database of reactions or compounds' + """Search for reactions and compounds in the model.""" @classmethod def init_parser(cls, parser): From 2f6123bfeb78591d0d43ecc4323846cc328ad8a9 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 15:19:21 -0400 Subject: [PATCH 32/46] search: Import print_function from __future__ to fix output --- psamm/commands/search.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/psamm/commands/search.py b/psamm/commands/search.py index 9a531bcd..7fdc06cb 100644 --- a/psamm/commands/search.py +++ b/psamm/commands/search.py @@ -15,6 +15,8 @@ # # Copyright 2014-2015 Jon Lund Steffensen +from __future__ import print_function + import re from ..command import Command From 58dd4f957920cde70b25ad2c3ff3dcb1cef3ee7b Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 16:42:39 -0400 Subject: [PATCH 33/46] command: Change name of debug environment variable to PSAMM_DEBUG We use a PSAMM_ prefix to not interfere with other software using the DEBUG variable. --- psamm/command.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psamm/command.py b/psamm/command.py index 0c37dc9c..ab9c9bc9 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -142,8 +142,8 @@ def main(command_class=None): """ # Set up logging for the command line interface - if 'DEBUG' in os.environ: - level = getattr(logging, os.environ['DEBUG'].upper(), None) + if 'PSAMM_DEBUG' in os.environ: + level = getattr(logging, os.environ['PSAMM_DEBUG'].upper(), None) if level is not None: logging.basicConfig(level=level) else: From a0b6cde0bc930d4da7ac2a2ecdfa45e05c0464de Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Thu, 27 Aug 2015 16:43:48 -0400 Subject: [PATCH 34/46] command: Use a simpler logging format when PSAMM_DEBUG is not set --- psamm/command.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/psamm/command.py b/psamm/command.py index ab9c9bc9..84b3e21e 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -147,7 +147,8 @@ def main(command_class=None): if level is not None: logging.basicConfig(level=level) else: - logging.basicConfig(level=logging.INFO) + logging.basicConfig( + level=logging.INFO, format='%(levelname)s: %(message)s') title = 'Metabolic modeling tools' if command_class is not None: From c9fa2e5980d4a36a25cb3e777131422d95d53fe1 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 31 Aug 2015 13:47:54 -0400 Subject: [PATCH 35/46] tox: Use Python3 Cplex package --- tox.ini | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tox.ini b/tox.ini index fa0527a6..18ddddc7 100644 --- a/tox.ini +++ b/tox.ini @@ -2,7 +2,7 @@ envlist = py27-{nosolver,cplex,qsoptex,gurobi}, py33-{nosolver,qsoptex}, - py34-{nosolver,qsoptex}, + py34-{nosolver,cplex,qsoptex}, coverage, flake, docs @@ -19,7 +19,8 @@ setenv = gurobi: PSAMM_SOLVER=gurobi deps = coverage - cplex: {env:CPLEX_PYTHON_PACKAGE} + py27-cplex: {env:CPLEX_PYTHON2_PACKAGE} + py34-cplex: {env:CPLEX_PYTHON3_PACKAGE} qsoptex: python-qsoptex>=0.4 gurobi: {env:GUROBI_PYTHON_PACKAGE} passenv = DYLD_LIBRARY_PATH From 694293ff58e694e9d04b61b9a825ef7bafd0c780 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 31 Aug 2015 13:48:29 -0400 Subject: [PATCH 36/46] docs: Update install instructions based on the latest Cplex release The location of the Cplex python package changed slightly in the new release because Python 3.4 is supported now with a separate package. Updates the install instructions for this change. Also, updates the instructions for running tox. --- docs/development.rst | 24 +++++++++++++++++++----- docs/install.rst | 13 ++++++++++--- 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/docs/development.rst b/docs/development.rst index 3ce94b75..7f0ac578 100644 --- a/docs/development.rst +++ b/docs/development.rst @@ -21,14 +21,28 @@ for coding style (PEP8) and building documentation, use tox_: $ tox -When testing with tox, the URL for the Cplex python module must be provided in -the environment variable ``CPLEX_PYTHON_PACKAGE`` for the Cplex tests to -succeed. For example: +When testing with tox, the local path for the Cplex python module must be +provided in the environment variables ``CPLEX_PYTHON2_PACKAGE`` and +``CPLEX_PYTHON3_PACKAGE`` for Python 2 and Python 3, respectively. For example: .. code-block:: shell - $ export CPLEX_PYTHON_PACKAGE=file:///path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/x86-64_osx - $ tox + $ export CPLEX_PYTHON2_PACKAGE=/path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/2.7/x86-64_osx + $ export CPLEX_PYTHON3_PACKAGE=/path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/3.4/x86-64_osx + $ tox -e py27-cplex,py34-cplex + +.. note:: + + Python 3 support was added in a recent release of Cplex. Older versions + only support Python 2. + +Similarly, the local path to the Gurobi package must be specified in the +environment variable ``GUROBI_PYTHON_PACKAGE``: + +.. code-block:: shell + + $ export GUROBI_PYTHON_PACKAGE=/Library/gurobi604/mac64 + $ tox -e py27-gurobi Adding new tests ---------------- diff --git a/docs/install.rst b/docs/install.rst index 8c135b23..6c363472 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -57,15 +57,22 @@ a virtual environment (as described above) this should be done after activating the virtual environment: 1. Locate the directory where Cplex was installed (e.g. ``/path/to/IBM/ILOG/CPLEX_StudioXXX``). -2. Locate the appropriate subdirectory based on your platform: - ``cplex/python/`` (e.g. ``cplex/python/x86-64_osx``). +2. Locate the appropriate subdirectory based on your platform and Python + version: ``cplex/python//`` + (e.g. ``cplex/python/2.7/x86-64_osx``). 3. Use ``pip`` to install the package from this directory using the following command. .. code-block:: shell (env) $ pip install \ - /path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/ + /path/to/IBM/ILOG/CPLEX_Studio1262/cplex/python// + +.. note:: + + Python 3 support was added in a recent release of Cplex. Older versions + only support Python 2. If you are using Python 3 make sure that you have + the latest version of Cplex installed. Gurobi ------ From 9e5517cbdc3ec12f246fa9efb7d174654f2c1396 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 31 Aug 2015 14:21:37 -0400 Subject: [PATCH 37/46] README: Update PSAMM description --- README.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/README.rst b/README.rst index 554abdbb..4d460cc4 100644 --- a/README.rst +++ b/README.rst @@ -11,8 +11,10 @@ PSAMM metabolic modeling tools :alt: Package Status :target: https://pypi.python.org/pypi/psamm -Tools related to metabolic modeling, reconstruction, data parsing and -formatting, consistency checking, automatic gap filling, and model simulations. +PSAMM is an open source software that is designed for the curation and analysis +of metabolic models. It supports model version tracking, model annotation, data +integration, data parsing and formatting, consistency checking, automatic gap +filling, and model simulations. See NEWS_ for information on recent changes. The ``master`` branch tracks the latest release while the ``develop`` branch is the latest version in From fa96c49dce198ba635f31d2f38281d2ef0e7f4a1 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Mon, 31 Aug 2015 14:21:54 -0400 Subject: [PATCH 38/46] docs: Update PSAMM description --- docs/overview.rst | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/docs/overview.rst b/docs/overview.rst index 49dafdcd..4ad26072 100644 --- a/docs/overview.rst +++ b/docs/overview.rst @@ -2,14 +2,15 @@ Overview ======== -PSAMM is a collection of tools related to metabolic modeling, model -reconstruction, data parsing and formatting, consistency checking, -automatic gap filling, and model simulations. +PSAMM is an open source software that is designed for the curation and analysis +of metabolic models. It supports model version tracking, model annotation, data +integration, data parsing and formatting, consistency checking, automatic gap +filling, and model simulations. PSAMM is developed as an open source project, coordinated through `Github`_. -The PSAMM software is being developed in the `Zhang Laboratory`_ at the University of Rhode Island and -is the subject of a research article that has been submitted for -review. [Steffensen15]_ +The PSAMM software is being developed in the `Zhang Laboratory`_ at the +University of Rhode Island and is the subject of a research article that has +been submitted for review. [Steffensen15]_ .. _Github: https://github.com/zhanglab/psamm .. _Zhang Laboratory: http://zhanglab.uri.edu/ From 57ddedd1abb527df8f120c8ceaca53d0db1d6145 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 12:01:57 -0400 Subject: [PATCH 39/46] sbml: Log warning on level 1 non-integer stoichiometry --- psamm/datasource/sbml.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/psamm/datasource/sbml.py b/psamm/datasource/sbml.py index b60afe21..30b25d9c 100644 --- a/psamm/datasource/sbml.py +++ b/psamm/datasource/sbml.py @@ -198,8 +198,12 @@ def _parse_species_references(self, name): denom = int(species.get('denominator', 1)) species_value = Fraction(value, denom) except ValueError: + message = ('Non-integer stoichiometry is not allowed in' + ' SBML level 1 (reaction {})'.format(self.id)) if self._reader._strict: - raise + raise ParseError(message) + else: + logger.warning(message) species_value = Decimal(species.get('stoichiometry', 1)) elif self._reader._level == 2: # Stoichiometric value is a double but can alternatively be From 60e1dbda05fb9857820539a01d890080aa3dbe5a Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 12:57:21 -0400 Subject: [PATCH 40/46] fluxanalysis: Merge into one FluxBalanceProblem object Remove the FluxBalanceTDProblem and move the contents into the FluxBalanceProblem as a method add_thermodynamic() that adds the necessary variables and thermodynamic contraints. The solve() method is renamed to maximize() and additional methods are added for minimize_l1() and the convenience method max_min_l1(). The method minimize_l1() applies the minimization of the L1 norm of fluxes previously found in flux_minimization() and flux_minimization() is now just a wrapper. The max_min_l1() method first uses maximize(), then contrains the maximized reaction to the optimal value, and finally call minimize_l1(). --- psamm/commands/randomsparse.py | 11 +- psamm/commands/robustness.py | 2 +- psamm/fluxanalysis.py | 258 +++++++++++++++++++------------ psamm/tests/test_fluxanalysis.py | 33 ++++ 4 files changed, 196 insertions(+), 108 deletions(-) diff --git a/psamm/commands/randomsparse.py b/psamm/commands/randomsparse.py index 4888d639..66a1a8d5 100644 --- a/psamm/commands/randomsparse.py +++ b/psamm/commands/randomsparse.py @@ -66,17 +66,18 @@ def run(self): reaction)) if not self._args.tfba: - fb_problem = fluxanalysis.FluxBalanceProblem solver = self._get_solver() else: - fb_problem = fluxanalysis.FluxBalanceTDProblem solver = self._get_solver(integer=True) - p = fb_problem(self._mm, solver) + p = fluxanalysis.FluxBalanceProblem(self._mm, solver) + + if self._args.tfba: + p.add_thermodynamic() threshold = self._args.threshold if threshold.relative: - p.solve(reaction) + p.maximize(reaction) threshold.reference = p.get_flux(reaction) flux_threshold = float(threshold) @@ -107,7 +108,7 @@ def run(self): testing_reaction)) try: - p.solve(reaction) + p.maximize(reaction) except fluxanalysis.FluxBalanceError: logger.info( 'FBA is infeasible, marking {} as essential'.format( diff --git a/psamm/commands/robustness.py b/psamm/commands/robustness.py index de130971..95678090 100644 --- a/psamm/commands/robustness.py +++ b/psamm/commands/robustness.py @@ -55,7 +55,7 @@ def run(self): def run_fba_fmin(model, reaction): fba = fluxanalysis.FluxBalanceProblem(model, solver=solver) - fba.solve(reaction) + fba.maximize(reaction) optimum = fba.get_flux(reaction) return fluxanalysis.flux_minimization( model, {reaction: optimum}, solver=solver) diff --git a/psamm/fluxanalysis.py b/psamm/fluxanalysis.py index 90c0c522..96609701 100644 --- a/psamm/fluxanalysis.py +++ b/psamm/fluxanalysis.py @@ -30,10 +30,10 @@ def _get_fba_problem(model, tfba, solver): """Convenience function for returning the right FBA problem instance""" + p = FluxBalanceProblem(model, solver) if tfba: - return FluxBalanceTDProblem(model, solver=solver) - else: - return FluxBalanceProblem(model, solver=solver) + p.add_thermodynamic() + return p class FluxBalanceError(Exception): @@ -41,24 +41,41 @@ class FluxBalanceError(Exception): class FluxBalanceProblem(object): - """Maximize the flux of a specific reaction + """Model as a flux optimization problem with steady state assumption. - The solve method solves the problem for a specific parameter. After - solving, the flux can be obtained using the get_flux method. The problem - can be solved again with a new objective as many times as needed. + Create a representation of the model as an LP optimization problem with + steady state assumption, i.e. the concentrations of compounds are always + zero. + + The problem can be modified and solved as many times as needed. The flux + of a reaction can be obtained after solving using :meth:`.get_flux`. + + Args: + model: :class:`MetabolicModel` to solve. + solver: LP solver instance to use. """ def __init__(self, model, solver): self._prob = solver.create_problem() + self._model = model + + self._has_minimization_vars = False + + # We keep track of temporary constraints from the last optimization so + # that we can remove them during the next call to solve(). This is + # necessary because removing the constraint immediately after solving + # a model can invalidate the solution in some solvers. + self._temp_constr = [] + self._remove_constr = [] # Define flux variables - for reaction_id in model.reactions: - lower, upper = model.limits[reaction_id] + for reaction_id in self._model.reactions: + lower, upper = self._model.limits[reaction_id] self._prob.define(('v', reaction_id), lower=lower, upper=upper) # Define constraints massbalance_lhs = {compound: 0 for compound in model.compounds} - for spec, value in iteritems(model.matrix): + for spec, value in iteritems(self._model.matrix): compound, reaction_id = spec massbalance_lhs[compound] += self.get_flux_var(reaction_id) * value for compound, lhs in iteritems(massbalance_lhs): @@ -66,11 +83,57 @@ def __init__(self, model, solver): @property def prob(self): - """Return the underlying LP problem""" + """Return the underlying LP problem. + + This can be used to add additional constraints on the problem. Calling + solve on the underlying problem is not guaranteed to work correctly, + instead use the methods on this object that solves the problem or + make a subclass with a method that calls :meth:`_solve`. + """ return self._prob - def solve(self, reaction): - """Solve problem maximizing the given reaction + def add_thermodynamic(self, epsilon=1e-5, em=1e5): + """Apply thermodynamic constraints to the model. + + Adding these constraints restricts the solution space to only + contain solutions that have no internal loops. This is solved as a + MILP problem as described in [Muller13]_. + """ + + p = self._prob + for reaction_id in self._model.reactions: + # Constrain internal reactions to a direction determined + # by alpha. + if not self._model.is_exchange(reaction_id): + p.define(('alpha', reaction_id), types=lp.VariableType.Binary) + p.define(('dmu', reaction_id)) # Delta mu + + flux = self.get_flux_var(reaction_id) + alpha = p.var(('alpha', reaction_id)) + dmu = p.var(('dmu', reaction_id)) + + lower, upper = self._model.limits[reaction_id] + p.add_linear_constraints(flux >= lower*(1 - alpha), + flux <= upper*alpha, + dmu >= -em*alpha + epsilon, + dmu <= em*(1 - alpha) - epsilon) + + # Define mu variables + p = self._prob + p.define(*(('mu', compound) for compound in self._model.compounds)) + + tdbalance_lhs = {reaction_id: 0 + for reaction_id in self._model.reactions} + for spec, value in iteritems(self._model.matrix): + compound, reaction_id = spec + if not self._model.is_exchange(reaction_id): + tdbalance_lhs[reaction_id] += p.var(('mu', compound)) * value + for reaction_id, lhs in iteritems(tdbalance_lhs): + if not self._model.is_exchange(reaction_id): + p.add_linear_constraints(lhs == p.var(('dmu', reaction_id))) + + def maximize(self, reaction): + """Solve the model by maximizing the given reaction. If reaction is a dictionary object, each entry is interpreted as a weight on the objective for that reaction (non-existent reaction will @@ -83,68 +146,85 @@ def solve(self, reaction): else: objective = self.get_flux_var(reaction) - # Set objective and solve self._prob.set_linear_objective(objective) - result = self._prob.solve(lp.ObjectiveSense.Maximize) - if not result: - raise FluxBalanceError('Non-optimal solution: {}'.format( - result.status)) - def get_flux_var(self, reaction): - """Get LP variable representing the reaction flux""" - return self._prob.var(('v', reaction)) + self._solve() - def get_flux(self, reaction): - """Get resulting flux value for reaction""" - return self._prob.result.get_value(('v', reaction)) + def _add_minimization_vars(self): + """Add variables and constraints for L1 norm minimization.""" + if self._has_minimization_vars: + return + var_names = (('z', reaction_id) + for reaction_id in self._model.reactions) + self._prob.define(*var_names, lower=0) -class FluxBalanceTDProblem(FluxBalanceProblem): - """Maximize the flux of a specific reaction with thermodynamic constraints + # Define constraints + v = self._prob.set(('v', rxnid) for rxnid in self._model.reactions) + z = self._prob.set(('z', rxnid) for rxnid in self._model.reactions) + self._prob.add_linear_constraints(z >= v, v >= -z) - Like FBA, but with the additional constraint that the flux must be - thermodynamically feasible. This is solved as a MILP problem and the - problem has been shown to be NP-hard in general. + self._has_minimization_vars = True - Described in [Muller13]_. - """ + def minimize_l1(self, weights={}): + """Solve the model by minimizing the L1 norm of the fluxes. - def __init__(self, model, solver, epsilon=1e-5, em=1e5): - super(FluxBalanceTDProblem, self).__init__(model, solver) - p = self.prob + If the weights dictionary is given, the weighted L1 norm if minimized + instead. The dictionary contains the weights of each reaction + (default 1). + """ - for reaction_id in model.reactions: - # Constrain internal reactions to a direction determined - # by alpha. - if not model.is_exchange(reaction_id): - p.define(('alpha', reaction_id), types=lp.VariableType.Binary) - p.define(('dmu', reaction_id)) # Delta mu + self._add_minimization_vars() - flux = self.get_flux_var(reaction_id) - alpha = p.var(('alpha', reaction_id)) - dmu = p.var(('dmu', reaction_id)) + objective = sum( + self._prob.var(('z', reaction_id)) * -weights.get(reaction_id, 1) + for reaction_id in self._model.reactions) + self._prob.set_linear_objective(objective) - lower, upper = model.limits[reaction_id] - p.add_linear_constraints(flux >= lower*(1 - alpha), - flux <= upper*alpha, - dmu >= -em*alpha + epsilon, - dmu <= em*(1 - alpha) - epsilon) + self._solve() - # Define mu variables - p.define(*(('mu', compound) for compound in model.compounds)) + def max_min_l1(self, reaction, weights={}): + """Maximize flux of reaction then minimize the L1 norm. - tdbalance_lhs = {reaction_id: 0 for reaction_id in model.reactions} - for spec, value in iteritems(model.matrix): - compound, reaction_id = spec - if not model.is_exchange(reaction_id): - tdbalance_lhs[reaction_id] += p.var(('mu', compound)) * value - for reaction_id, lhs in iteritems(tdbalance_lhs): - if not model.is_exchange(reaction_id): - p.add_linear_constraints(lhs == p.var(('dmu', reaction_id))) + During minimization the given reaction will be fixed at the maximum + obtained from the first solution. + """ + + self.maximize(reaction) + flux = self.get_flux(reaction) + flux_var = self.get_flux_var(reaction) + c, = self._prob.add_linear_constraints(flux_var == flux) + self._temp_constr.append(c) + self.minimize_l1(weights) + + def _solve(self): + """Solve the problem with the current objective.""" + + # Remove temporary constraints + while len(self._remove_constr) > 0: + self._remove_constr.pop().delete() + + try: + result = self._prob.solve(lp.ObjectiveSense.Maximize) + if not result: + raise FluxBalanceError('Non-optimal solution: {}'.format( + result.status)) + finally: + # Set temporary constraints to be removed on next solve call + self._remove_constr = self._temp_constr + self._temp_constr = [] + + def get_flux_var(self, reaction): + """Get LP variable representing the reaction flux.""" + return self._prob.var(('v', reaction)) + + def get_flux(self, reaction): + """Get resulting flux value for reaction.""" + return self._prob.result.get_value(('v', reaction)) def flux_balance(model, reaction, tfba, solver): - """Run flux balance analysis on the given model + """Run flux balance analysis on the given model. Yields the reaction id and flux value for each reaction in the model. @@ -165,13 +245,13 @@ def flux_balance(model, reaction, tfba, solver): """ fba = _get_fba_problem(model, tfba, solver) - fba.solve(reaction) + fba.maximize(reaction) for reaction in model.reactions: yield reaction, fba.get_flux(reaction) def flux_variability(model, reactions, fixed, tfba, solver): - """Find the variability of each reaction while fixing certain fluxes + """Find the variability of each reaction while fixing certain fluxes. Yields the reaction id, and a tuple of minimum and maximum value for each of the given reactions. The fixed reactions are given in a dictionary as @@ -197,7 +277,7 @@ def flux_variability(model, reactions, fixed, tfba, solver): def min_max_solve(reaction_id): for direction in (-1, 1): - fba.solve({reaction_id: direction}) + fba.maximize({reaction_id: direction}) yield fba.get_flux(reaction_id) # Solve for each reaction @@ -206,7 +286,7 @@ def min_max_solve(reaction_id): def flux_minimization(model, fixed, solver, weights={}): - """Minimize flux of all reactions while keeping certain fluxes fixed + """Minimize flux of all reactions while keeping certain fluxes fixed. The fixed reactions are given in a dictionary as reaction id to value mapping. The weighted L1-norm of the fluxes is minimized. @@ -221,46 +301,20 @@ def flux_minimization(model, fixed, solver, weights={}): An iterator of reaction ID and reaction flux pairs. """ - prob = solver.create_problem() + fba = FluxBalanceProblem(model, solver) - # Define flux variables - for reaction_id in model.reactions: - lower, upper = model.limits[reaction_id] - if reaction_id in fixed: - lower = max(lower, fixed[reaction_id]) - prob.define(('v', reaction_id), lower=lower, upper=upper) - - # Define z variables - prob.define(*(('z', reaction_id) for reaction_id in model.reactions), - lower=0) - objective = sum(prob.var(('z', reaction_id)) * weights.get(reaction_id, 1) - for reaction_id in model.reactions) - prob.set_linear_objective(objective) - - # Define constraints - v = prob.set(('v', rxnid) for rxnid in model.reactions) - z = prob.set(('z', rxnid) for rxnid in model.reactions) - prob.add_linear_constraints(z >= v, v >= -z) - - massbalance_lhs = {compound: 0 for compound in model.compounds} - for spec, value in iteritems(model.matrix): - compound, rxnid = spec - massbalance_lhs[compound] += value * prob.var(('v', rxnid)) - for compound, lhs in iteritems(massbalance_lhs): - prob.add_linear_constraints(lhs == 0) - - # Solve - result = prob.solve(lp.ObjectiveSense.Minimize) - if not result: - raise FluxBalanceError('Non-optimal solution: {}'.format( - result.status)) - - return ((reaction_id, result.get_value(('v', reaction_id))) + for reaction_id, value in iteritems(fixed): + flux = fba.get_flux_var(reaction_id) + fba.prob.add_linear_constraints(flux >= value) + + fba.minimize_l1() + + return ((reaction_id, fba.get_flux(reaction_id)) for reaction_id in model.reactions) def flux_randomization(model, fixed, tfba, solver): - """Find a uniformly random flux mode in the solution space + """Find a uniformly random flux mode in the solution space. This can be used to generate a random sample from the solution space. If many random samples are needed more efficient methods @@ -290,13 +344,13 @@ def flux_randomization(model, fixed, tfba, solver): for reaction_id, value in iteritems(fixed): fba.prob.add_linear_constraints(fba.get_flux_var(reaction_id) >= value) - fba.solve(optimize) + fba.maximize(optimize) for reaction_id in model.reactions: yield reaction_id, fba.get_flux(reaction_id) def consistency_check(model, subset, epsilon, tfba, solver): - """Check that reaction subset of model is consistent using FBA + """Check that reaction subset of model is consistent using FBA. Yields all reactions that are *not* flux consistent. A reaction is consistent if there is at least one flux solution to the model that both @@ -330,14 +384,14 @@ def consistency_check(model, subset, epsilon, tfba, solver): logger.debug('{} left, checking {}...'.format(len(subset), reaction)) - fba.solve(reaction) + fba.maximize(reaction) support = set(rxnid for rxnid in model.reactions if abs(fba.get_flux(rxnid)) >= epsilon) subset -= support if reaction in support: continue elif model.is_reversible(reaction): - fba.solve({reaction: -1}) + fba.maximize({reaction: -1}) support = set(rxnid for rxnid in model.reactions if abs(fba.get_flux(rxnid)) >= epsilon) subset -= support diff --git a/psamm/tests/test_fluxanalysis.py b/psamm/tests/test_fluxanalysis.py index e1ad0805..ed864971 100755 --- a/psamm/tests/test_fluxanalysis.py +++ b/psamm/tests/test_fluxanalysis.py @@ -69,6 +69,39 @@ def test_flux_balance_rxn_6(self): self.assertAlmostEqual(fluxes['rxn_2'], 0) self.assertAlmostEqual(fluxes['rxn_6'], 1000) + def test_flux_balance_object_maximize(self): + p = fluxanalysis.FluxBalanceProblem(self.model, self.solver) + p.maximize('rxn_6') + self.assertAlmostEqual(p.get_flux('rxn_1'), 500) + self.assertAlmostEqual(p.get_flux('rxn_2'), 0) + self.assertAlmostEqual(p.get_flux('rxn_6'), 1000) + + def test_flux_balance_object_minimize_l1(self): + p = fluxanalysis.FluxBalanceProblem(self.model, self.solver) + p.prob.add_linear_constraints(p.get_flux_var('rxn_6') == 1000) + p.minimize_l1() + self.assertAlmostEqual(p.get_flux('rxn_1'), 500) + self.assertAlmostEqual(p.get_flux('rxn_2'), 0) + self.assertAlmostEqual(p.get_flux('rxn_3'), 1000) + self.assertAlmostEqual(p.get_flux('rxn_4'), 0) + self.assertAlmostEqual(p.get_flux('rxn_5'), 0) + self.assertAlmostEqual(p.get_flux('rxn_6'), 1000) + + def test_flux_balance_object_max_min_l1(self): + p = fluxanalysis.FluxBalanceProblem(self.model, self.solver) + p.max_min_l1('rxn_6') + self.assertAlmostEqual(p.get_flux('rxn_1'), 500) + self.assertAlmostEqual(p.get_flux('rxn_2'), 0) + self.assertAlmostEqual(p.get_flux('rxn_3'), 1000) + self.assertAlmostEqual(p.get_flux('rxn_4'), 0) + self.assertAlmostEqual(p.get_flux('rxn_5'), 0) + self.assertAlmostEqual(p.get_flux('rxn_6'), 1000) + + # The temporary constraint on the reaction rxn_6 should go away. If + # not, the next maximize will raise a FluxBalanceError. + p.prob.add_linear_constraints(p.get_flux_var('rxn_1') == 10) + p.maximize('rxn_6') + self.assertAlmostEqual(p.get_flux('rxn_6'), 20) class TestFluxBalanceThermodynamic(unittest.TestCase): def setUp(self): From 20b5e6bc2e0fd3ba136ef452ada2bf1cf616f833 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 13:29:06 -0400 Subject: [PATCH 41/46] fba: Use new interface to speed up fba command slightly --- psamm/commands/fba.py | 40 ++++++++++++++++++++++++---------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/psamm/commands/fba.py b/psamm/commands/fba.py index d9df1f7c..2ea0a9d3 100644 --- a/psamm/commands/fba.py +++ b/psamm/commands/fba.py @@ -40,7 +40,7 @@ def init_parser(cls, parser): super(FluxBalanceCommand, cls).init_parser(parser) def run(self): - """Run flux analysis command""" + """Run flux analysis command.""" # Load compound information compound_name = {} @@ -79,33 +79,41 @@ def run(self): logger.info('Maximum flux: {}'.format(optimum)) def run_fba_minimized(self, reaction): - """Run normal FBA and flux minimization on model, then print output""" + """Run normal FBA and flux minimization on model.""" - solver = self._get_solver() - fba_fluxes = dict(fluxanalysis.flux_balance( - self._mm, reaction, tfba=False, solver=solver)) - optimum = fba_fluxes[reaction] epsilon = self._args.epsilon + solver = self._get_solver() + p = fluxanalysis.FluxBalanceProblem(self._mm, solver) + + # Maximize reaction flux + p.maximize(reaction) + fluxes = {r: p.get_flux(r) for r in self._mm.reactions} + # Run flux minimization - fmin_fluxes = dict(fluxanalysis.flux_minimization( - self._mm, {reaction: optimum}, solver=solver)) + flux_var = p.get_flux_var(reaction) + p.prob.add_linear_constraints(flux_var == p.get_flux(reaction)) + p.minimize_l1() count = 0 - for reaction_id, flux in iteritems(fmin_fluxes): - if fba_fluxes[reaction_id] - epsilon > flux: + for reaction_id in self._mm.reactions: + flux = p.get_flux(reaction_id) + if abs(flux - fluxes[reaction_id]) > epsilon: count += 1 - yield reaction_id, fba_fluxes[reaction_id], flux + yield reaction_id, fluxes[reaction_id], flux logger.info('Minimized reactions: {}'.format(count)) def run_tfba(self, reaction): - """Run FBA and tFBA on model""" + """Run FBA and tFBA on model.""" solver = self._get_solver(integer=True) + p = fluxanalysis.FluxBalanceProblem(self._mm, solver) + + p.maximize(reaction) + fba_fluxes = {r: p.get_flux(r) for r in self._mm.reactions} - fba_fluxes = dict(fluxanalysis.flux_balance( - self._mm, reaction, tfba=False, solver=solver)) - fluxes = dict(fluxanalysis.flux_balance( - self._mm, reaction, tfba=True, solver=solver)) + p.add_thermodynamic() + p.maximize(reaction) + fluxes = {r: p.get_flux(r) for r in self._mm.reactions} for reaction_id, flux in iteritems(fluxes): yield reaction_id, fba_fluxes[reaction_id], flux From 7814bc500ff363944dc927599ba5b308662b7258 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 13:28:06 -0400 Subject: [PATCH 42/46] fastgapfill: Fix error checking whether biomass reaction is in model --- psamm/commands/fastgapfill.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psamm/commands/fastgapfill.py b/psamm/commands/fastgapfill.py index 6c5f2ca0..87ee77ac 100644 --- a/psamm/commands/fastgapfill.py +++ b/psamm/commands/fastgapfill.py @@ -114,7 +114,7 @@ def run(self): if maximized_reaction is None: raise CommandError('The maximized reaction was not specified') - if not not self._mm.has_reaction(maximized_reaction): + if not self._mm.has_reaction(maximized_reaction): raise CommandError( 'The biomass reaction is not a valid model' ' reaction: {}'.format(maximized_reaction)) From 680f760796b57fa4843dd0e2176afbe8c537dfe1 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 13:30:48 -0400 Subject: [PATCH 43/46] robustness: Use new interface to speed up command --- psamm/commands/robustness.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/psamm/commands/robustness.py b/psamm/commands/robustness.py index 95678090..443eb564 100644 --- a/psamm/commands/robustness.py +++ b/psamm/commands/robustness.py @@ -53,22 +53,9 @@ def init_parser(cls, parser): def run(self): """Run flux analysis command""" - def run_fba_fmin(model, reaction): - fba = fluxanalysis.FluxBalanceProblem(model, solver=solver) - fba.maximize(reaction) - optimum = fba.get_flux(reaction) - return fluxanalysis.flux_minimization( - model, {reaction: optimum}, solver=solver) - - def run_tfba(model, reaction): - return fluxanalysis.flux_balance( - model, reaction, tfba=True, solver=solver) - if self._args.no_tfba: - run_fba = run_fba_fmin solver = self._get_solver() else: - run_fba = run_tfba solver = self._get_solver(integer=True) # Load compound information @@ -111,14 +98,27 @@ def run_tfba(model, reaction): raise CommandError('Invalid flux range: {}, {}\n'.format( flux_min, flux_max)) + # Remove the limits on the varying reaction. We will instead fix this + # reaction explicitly in the following loop. + del self._mm.limits[varying_reaction].bounds + + p = fluxanalysis.FluxBalanceProblem(self._mm, solver) + for i in range(steps): fixed_flux = flux_min + i*(flux_max - flux_min)/float(steps-1) - test_model = self._mm.copy() - test_model.limits[varying_reaction].bounds = fixed_flux, fixed_flux + flux_var = p.get_flux_var(varying_reaction) + c, = p.prob.add_linear_constraints(flux_var == fixed_flux) try: - for other_reaction, flux in run_fba(test_model, reaction): + if self._args.no_tfba: + p.max_min_l1(reaction) + else: + p.maximize(reaction) + for other_reaction in self._mm.reactions: print('{}\t{}\t{}'.format( - other_reaction, fixed_flux, flux)) + other_reaction, fixed_flux, + p.get_flux(other_reaction))) except fluxanalysis.FluxBalanceError: pass + finally: + c.delete() From 61aecadda483648492bf9bbaa18083187253e300 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 12:52:51 -0400 Subject: [PATCH 44/46] robustness: Use range from six for Py3 compatability --- psamm/commands/robustness.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/psamm/commands/robustness.py b/psamm/commands/robustness.py index 768450b7..de130971 100644 --- a/psamm/commands/robustness.py +++ b/psamm/commands/robustness.py @@ -18,6 +18,8 @@ from ..command import Command, SolverCommandMixin, CommandError from .. import fluxanalysis +from six.moves import range + class RobustnessCommand(SolverCommandMixin, Command): """Run robustness analysis on the model. @@ -109,7 +111,7 @@ def run_tfba(model, reaction): raise CommandError('Invalid flux range: {}, {}\n'.format( flux_min, flux_max)) - for i in xrange(steps): + for i in range(steps): fixed_flux = flux_min + i*(flux_max - flux_min)/float(steps-1) test_model = self._mm.copy() test_model.limits[varying_reaction].bounds = fixed_flux, fixed_flux From 2a8883742e1edb16797e6c4f4aa7fee96c4844fa Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 15:19:19 -0400 Subject: [PATCH 45/46] setup.py: Bump version to 0.12 --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index bebff9f8..46f16fc9 100755 --- a/setup.py +++ b/setup.py @@ -24,7 +24,7 @@ setup( name='psamm', - version='0.11', + version='0.12', description='PSAMM metabolic modeling tools', maintainer='Jon Lund Steffensen', maintainer_email='jon_steffensen@uri.edu', From 7478ca7e5b7789e273dce25796f6fa3fb13a2dee Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 1 Sep 2015 15:19:31 -0400 Subject: [PATCH 46/46] NEWS: Update for 0.12 --- NEWS.md | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/NEWS.md b/NEWS.md index 97cd6696..8340a084 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,37 @@ +v0.12 (2015-09-01) +------------------ + +- Add support for the Gurobi LP solver (Python 2.7 only). +- Fix a bug in the `fastgapfill` command that caused the biomass reaction to + never be found in the model. +- Add command line tool `psamm-list-lpsolvers` which will list the available + solvers, and report their attributes and whether they can be successfully + loaded. +- Change the way the threshold parameter on the `randomsparse` command is + parsed. A threshold relative to the maximum flux is now given using + percentage notation (e.g. `95%`) while an absolute flux threshold is given as + a number (e.g `1.34`). +- Log the model name at the beginning of every command. If the model exists in + a Git repository, the current Git commit ID is also reported. +- Produce warnings when encountering invalid SBML constructs using non-strict + parsing. Previously, some of these could only be noticed when using the + strict mode. +- Improve error reports from commands. This means that usage information is now + printed when a command fails to parse an argument. +- Improve the API in the `fluxanalysis` module to allow commands to more + easily reuse LP problems. This improves the speed of robustness analysis in + particular. +- Improve and expand the install instructions in the documentation. +- Reorganize implementations of commands into the `psamm.commands` package. The + commands are now discovered through the entry point mechanism from + `setuptools` allowing any package to install additional commands into the + PSAMM user interface. +- Use a simpler output format for log messages unless the environment variable + `PSAMM_DEBUG` is set. If set, this variable determines the logging level + (e.g. `PSAMM_DEBUG=debug` enables more messages). +- Include Python 3.3 in the tox test suite. + v0.11 (2015-08-06) ------------------