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 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) ------------------ 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 diff --git a/docs/api.rst b/docs/api.rst index 131495bb..1f6c5bee 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -4,25 +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_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/api/lpsolver_gurobi.rst b/docs/api/lpsolver_gurobi.rst new file mode 100644 index 00000000..e55c912a --- /dev/null +++ b/docs/api/lpsolver_gurobi.rst @@ -0,0 +1,6 @@ + +``psamm.lpsolver.gurobi`` -- Gurobi LP solver +================================================= + +.. automodule:: psamm.lpsolver.gurobi + :members: 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 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``) ----------------------------------------- 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/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/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 diff --git a/docs/install.rst b/docs/install.rst index 306429b3..6c363472 100644 --- a/docs/install.rst +++ b/docs/install.rst @@ -2,45 +2,52 @@ 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 + +.. code-block:: shell -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. + $ 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 ------------ -- 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. + +.. _install-cplex: Cplex ----- @@ -50,9 +57,37 @@ 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``). -3. Use ``pip`` to install the package from this directory: ``pip install /path/to/IBM/ILOG/CPLEX_StudioXXX/cplex/python/`` +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_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 +------ + +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 -------- @@ -63,7 +98,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 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/ diff --git a/psamm/command.py b/psamm/command.py index 4bd64e92..84b3e21e 100644 --- a/psamm/command.py +++ b/psamm/command.py @@ -14,46 +14,45 @@ # 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__) +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. @@ -72,6 +71,15 @@ 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)) + + 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(): @@ -126,1185 +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 ValueError('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)) - - 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 ValueError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise ValueError('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 ValueError('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 ValueError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise ValueError('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. - """ - - 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', nargs='?') - parser.add_argument( - 'threshold', help='Relative threshold of max reaction flux', - type=float) - 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 ValueError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise ValueError('Specified reaction is not in model: {}'.format( - reaction)) - - threshold = self._args.threshold - if threshold < 0.0 or threshold > 1.0: - raise ValueError( - 'Invalid threshold, must be in [0;1]: {}'.format(threshold)) - - 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.solve(reaction) - flux_threshold = p.get_flux(reaction) * threshold - - logger.info('Flux threshold for {} is {}'.format( - reaction, flux_threshold)) - - if self._args.exchange: - model_test = self._mm.copy() - essential = {reaction} - deleted = set() - exchange = set() - for reaction_id in self._mm.reactions: - if self._mm.is_exchange(reaction_id): - exchange.add(reaction_id) - test_set = set(exchange) - essential - else: - model_test = self._mm.copy() - essential = {reaction} - deleted = set() - test_set = set(self._mm.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 - - 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 - essential.add(testing_reaction) - continue - - logger.debug('Reaction {} has flux {}'.format( - reaction, p.get_flux(reaction))) - - if p.get_flux(reaction) < flux_threshold: - model_test.limits[testing_reaction].bounds = saved_bounds - 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)) - - 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)) - - -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 ValueError('The biomass reaction was not specified') - - if not self._mm.has_reaction(reaction): - raise ValueError('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( - varying_reaction)) - - steps = self._args.steps - if steps <= 0: - raise ValueError('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 ValueError('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`. @@ -1313,12 +142,13 @@ 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: - 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: @@ -1338,22 +168,22 @@ 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__ + 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=title, description=description) + name, help=title.rstrip('.'), + description=command_class.__doc__) subparser.set_defaults(command=command_class) command_class.init_parser(subparser) @@ -1364,8 +194,7 @@ def main(command_class=None): # Instantiate command with model and run command = args.command(model, args) - command.run() - - -if __name__ == '__main__': - main() + try: + command.run() + except CommandError as e: + parser.error(str(e)) 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..7255cb68 --- /dev/null +++ b/psamm/commands/chargecheck.py @@ -0,0 +1,79 @@ +# 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 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. + """ + + 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..37f40c94 --- /dev/null +++ b/psamm/commands/console.py @@ -0,0 +1,66 @@ +# 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 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..87ee77ac --- /dev/null +++ b/psamm/commands/fastgapfill.py @@ -0,0 +1,147 @@ +# 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 the FastGapFill gap-filling algorithm on 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 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..2ea0a9d3 --- /dev/null +++ b/psamm/commands/fba.py @@ -0,0 +1,119 @@ +# 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 the 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.""" + + 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 + 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 in self._mm.reactions: + flux = p.get_flux(reaction_id) + if abs(flux - fluxes[reaction_id]) > epsilon: + count += 1 + yield reaction_id, 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) + p = fluxanalysis.FluxBalanceProblem(self._mm, solver) + + p.maximize(reaction) + fba_fluxes = {r: p.get_flux(r) for r in self._mm.reactions} + + 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 diff --git a/psamm/commands/fluxcheck.py b/psamm/commands/fluxcheck.py new file mode 100644 index 00000000..08ca6a08 --- /dev/null +++ b/psamm/commands/fluxcheck.py @@ -0,0 +1,123 @@ +# 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. + """ + + @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..24922a2c --- /dev/null +++ b/psamm/commands/fluxcoupling.py @@ -0,0 +1,135 @@ +# 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). + """ + + 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..cde99a73 --- /dev/null +++ b/psamm/commands/formulacheck.py @@ -0,0 +1,95 @@ +# 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 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. + """ + + 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..d938e11e --- /dev/null +++ b/psamm/commands/fva.py @@ -0,0 +1,72 @@ +# 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 the 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..6450e75d --- /dev/null +++ b/psamm/commands/gapfill.py @@ -0,0 +1,77 @@ +# 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 the GapFind and GapFill algorithms on the 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..f14c7f0f --- /dev/null +++ b/psamm/commands/masscheck.py @@ -0,0 +1,112 @@ +# 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 the model is mass consistent.""" + + @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..66a1a8d5 --- /dev/null +++ b/psamm/commands/randomsparse.py @@ -0,0 +1,134 @@ +# 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 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. + 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%'). + """ + + @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: + solver = self._get_solver() + else: + solver = self._get_solver(integer=True) + + p = fluxanalysis.FluxBalanceProblem(self._mm, solver) + + if self._args.tfba: + p.add_thermodynamic() + + threshold = self._args.threshold + if threshold.relative: + p.maximize(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.maximize(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..443eb564 --- /dev/null +++ b/psamm/commands/robustness.py @@ -0,0 +1,124 @@ +# 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 + +from six.moves import range + + +class RobustnessCommand(SolverCommandMixin, Command): + """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 + 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. + """ + + @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""" + + if self._args.no_tfba: + solver = self._get_solver() + else: + 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)) + + # 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) + flux_var = p.get_flux_var(varying_reaction) + c, = p.prob.add_linear_constraints(flux_var == fixed_flux) + + try: + 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, + p.get_flux(other_reaction))) + except fluxanalysis.FluxBalanceError: + pass + finally: + c.delete() diff --git a/psamm/commands/sbmlexport.py b/psamm/commands/sbmlexport.py new file mode 100644 index 00000000..f7f7f697 --- /dev/null +++ b/psamm/commands/sbmlexport.py @@ -0,0 +1,30 @@ +# 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.""" + + 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..7fdc06cb --- /dev/null +++ b/psamm/commands/search.py @@ -0,0 +1,144 @@ +# 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 __future__ import print_function + +import re + +from ..command import Command +from ..reaction import Compound + + +class SearchCommand(Command): + """Search for reactions and compounds in the model.""" + + @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/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 diff --git a/psamm/datasource/sbml.py b/psamm/datasource/sbml.py index 890e043c..30b25d9c 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) @@ -193,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 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/lpsolver/generic.py b/psamm/lpsolver/generic.py index 0e69c945..1cf4257d 100644 --- a/psamm/lpsolver/generic.py +++ b/psamm/lpsolver/generic.py @@ -17,8 +17,9 @@ """Generic interface to LP solver instantiation.""" -from __future__ import absolute_import +from __future__ import absolute_import, print_function +import argparse import os import logging @@ -28,6 +29,7 @@ logger = logging.getLogger(__name__) _solvers = [] +_solver_import_errors = {} # Try to load Cplex solver @@ -40,8 +42,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,14 +55,38 @@ 'rational': True, 'priority': 5 }) -except ImportError: - pass +except ImportError as e: + _solver_import_errors['qsoptex'] = str(e) + +# Try to load Gurobi solver +try: + from . import gurobi + _solvers.append({ + 'class': gurobi.Solver, + 'name': 'gurobi', + 'integer': True, + 'rational': False, + 'priority': 9 + }) +except ImportError as e: + _solver_import_errors['gurobi'] = str(e) 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 @@ -73,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 = list(filter_solvers(_solvers, self._requirements)) # Obtain solver priority from environment variable, if specified. priority = {} @@ -126,3 +149,73 @@ def parse_solver_setting(s): value = float(value) return key, value + + +def list_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( + solver['integer'])) + print('Rational solution: {}'.format(solver['rational'])) + print('Class: {}'.format(solver['class'])) + print() + else: + print('No solvers fullfil the requirements!') + print() + + 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)) 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/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 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 diff --git a/psamm/tests/test_fluxanalysis.py b/psamm/tests/test_fluxanalysis.py index c68358d0..ed864971 100755 --- a/psamm/tests/test_fluxanalysis.py +++ b/psamm/tests/test_fluxanalysis.py @@ -45,30 +45,63 @@ 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) + + 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): @@ -94,11 +127,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 +160,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 +178,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): 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 606768af..83ad67c3 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,91 @@ 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. + + 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, + stderr=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. diff --git a/setup.py b/setup.py index 8494ef75..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', @@ -37,13 +37,31 @@ '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' ], 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' + ], + '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' ] }, diff --git a/tox.ini b/tox.ini index b6faa4e7..18ddddc7 100644 --- a/tox.ini +++ b/tox.ini @@ -1,7 +1,8 @@ [tox] envlist = - py27-{nosolver,cplex,qsoptex}, - py34-{nosolver,qsoptex}, + py27-{nosolver,cplex,qsoptex,gurobi}, + py33-{nosolver,qsoptex}, + py34-{nosolver,cplex,qsoptex}, coverage, flake, docs @@ -15,12 +16,16 @@ setenv = nosolver: PSAMM_SOLVER=nosolver cplex: PSAMM_SOLVER=cplex qsoptex: PSAMM_SOLVER=qsoptex + 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 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]