Merge branch 'release'

NEWS.md

@@ -1,4 +1,22 @@
 
+v0.30 (2017-06-23)
+------------------
+
+- Adds the new command `primarypairs` for predicting reactant/product element
+  transfers using the new FindPrimaryPairs method as well as the MapMaker
+  method.
+- A new option has been added to the `genedelete` command which allows use of
+  _minimization of metabolic adjustments_ to simulate biomass production for
+  gene knockouts.
+- Fixes a bug where the epsilon parameter was accidentally ignored by
+  `fastgapfill`.
+- Fixes a bug where the `psamm-sbml-model` command did not ignore boundary
+  species. With this change, the boundary species are also ignored by default
+  when using the API to read SBML models.
+- Fixes a performance issues with gap-filling that made the `gapfill` and
+  `fastgapfill` commands take much longer time to run on large models than
+  necessary.
+
 v0.29 (2017-04-18)
 ------------------
 

docs/api/moma.rst

@@ -0,0 +1,6 @@
+
+``psamm.moma`` -- Minimization of metabolic adjustments
+=======================================================
+
+.. automodule:: psamm.moma
+   :members:

docs/commands.rst

@@ -209,6 +209,48 @@ The file gene_file.txt would contain the following lines::
     ExGene1
     ExGene2
 
+To delete genes using different algorithms use ``--method`` to specify
+which algorithm for the solver to use. The default method for the command is
+FBA. To delete genes using the Minimization of Metabolic Adjustment (MOMA)
+algorithm use the command line argument ``--method moma``. MOMA is based on
+the assumption that the knockout organism has not had time to adjust its gene
+regulation to maximize biomass production so fluxes will be close to wildtype
+fluxes.
+
+.. code-block:: shell
+
+    $ psamm-model genedelete --gene ExGene1 --method moma
+
+There are four variations of MOMA available in PSAMM, defined in the following
+way (where :math:`\bar{v}` is the wild type fluxes and :math:`\bar{u}` is the
+knockout fluxes):
+
+MOMA (``--method moma``)
+    Finds the reaction fluxes in the knockout, such that the difference in flux
+    from the wildtype is minimized. Minimization is performed with the
+    Euclidean distance: :math:`\sum_j (v_j - u_j)^2`. The wildtype fluxes are
+    obtained from the wildtype model (i.e. before genes are deleted) by FBA
+    with L1 minimization. L1 minimization is performed on the FBA result to
+    remove loops and make the result disregard internal loop fluxes. [Segre02]_
+
+Linear MOMA (``--method lin_moma``)
+    Finds the reaction fluxes in the knockout, such that the difference in flux
+    from the wildtype is minimized. Minimization is performed with the
+    Manhattan distance: :math:`\sum_j \|v_j - u_j\|`. The wildtype fluxes are
+    obtained from the wildtype model (i.e. before genes are deleted) by FBA
+    with L1 minimization. L1 minimization is performed on the FBA result to
+    remove loops and make the result disregard internal loop fluxes. [Mo09]_
+
+Combined-model MOMA (``--method moma2``) (Experimental)
+    Similar to ``moma``, but this implementation solves for the wild type
+    fluxes at the same time as the knockout fluxes to ensure not to rely on the
+    arbitrary flux vector found with FBA.
+
+Combined-model linear MOMA (``--method lin_moma2``) (Experimental)
+    Similar to ``lin_moma``, but this implementation solves for the wild type
+    fluxes at the same time as the knockout fluxes to ensure not to rely on the
+    arbitrary flux vector found with FBA.
+
 Flux coupling analysis (``fluxcoupling``)
 -----------------------------------------
 
@@ -457,6 +499,25 @@ minimal solution.
 
     $ psamm-model fastgapfill --penalty penalty.tsv
 
+Predict primary pairs (``primarypairs``)
+----------------------------------------------------------------------
+
+This command is used to predict element-transferring reactant/product pairs
+in the reactions of the model. This can be used to determine the flow of
+elements through reactions. Two methods for predicting the pairs are available:
+FindPrimaryPairs (``fpp``) [Steffensen17]_ and
+MapMaker (``mapmaker``) [Tervo16]_. The ``--method`` option can used to select
+which prediction method to use:
+
+.. code-block:: shell
+
+    $ psamm-model primarypairs --method=fpp
+
+The result is reported as a table of four columns. The first column is the
+reactions ID, the second and third columns contain the compound ID of the
+reactant and product. The fourth column contains the predicted transfer of
+elements.
+
 SBML Export (``sbmlexport``)
 ----------------------------
 

docs/conf.py

@@ -67,8 +67,8 @@
 
 # General information about the project.
 project = u'PSAMM'
-copyright = u'2015, Jon Lund Steffensen'
-author = u'Jon Lund Steffensen'
+copyright = u'2015-2017, PSAMM developers'
+author = u'PSAMM developers'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -249,7 +249,7 @@
 #  author, documentclass [howto, manual, or own class]).
 latex_documents = [
   (master_doc, 'psamm.tex', u'PSAMM Documentation',
-   u'Jon Lund Steffensen', 'manual'),
+   u'PSAMM developers', 'manual'),
 ]
 
 # The name of an image file (relative to this directory) to place at the top of

docs/references.rst

@@ -20,6 +20,9 @@ References
 .. [Mahadevan03] Mahadevan R, Schilling CH. The effects of alternate optimal
     solutions in constraint-based genome-scale metabolic models. Metab Eng.
     2003;5: 264–276. :doi:`10.1016/j.ymben.2003.09.002`.
+.. [Mo09] Mo ML, Palsson BØ, Herrgård MJ. Connecting extracellular metabolomic
+    measurements to intracellular flux states in yeast. BMC Systems Biology.
+    2009;3(1):3-37. :doi:`10.1186/1752-0509-3-37`.
 .. [Muller13] Müller AC, Bockmayr A. Fast thermodynamically constrained flux
     variability analysis. Bioinformatics. 2013;29: 903–909.
     :doi:`10.1093/bioinformatics/btt059`.
@@ -30,6 +33,16 @@ References
     definition of metabolic pathways and their use in interpreting metabolic
     function from a pathway-oriented perspective. J Theor Biol. 2000;203:
     229–48. :doi:`10.1006/jtbi.2000.1073`.
+.. [Segre02] Segrè D, Vitkup D, Church GM. Analysis of optimality in natural
+    and perturbed metabolic networks. Proceedings of the National Academy of
+    Sciences of the United States of America. 2002;99(23):15112-15117.
+    :doi:`10.1073/pnas.232349399`.
+.. [Steffensen17] Steffensen JL, Zhang Y. FindPrimaryPairs: An efficient
+    algorithm for predicting element-transferring reactant/product pairs in
+    metabolic networks. *Submitted*
+.. [Tervo16] Tervo CJ, Reed JL. MapMaker and PathTracer for tracking carbon in
+    genome-scale metabolic models. Biotechnol J. 2016; 1–23.
+    :doi:`10.1002/biot.201400305`.
 .. [Thiele14] Thiele I, Vlassis N, Fleming RMT. fastGapFill: efficient gap
     filling in metabolic networks. Bioinformatics. 2014;30: 2529–31.
     :doi:`10.1093/bioinformatics/btu321`.
@@ -43,4 +56,4 @@ References
 .. [Orth11] Orth JD, Conrad TM, Na J, Lerman JA, Nam H, Feist AM, et al. A
     comprehensive genome-scale reconstruction of Escherichia coli
     metabolism--2011. Mol Syst Biol. EMBO Press; 2011;7: 535.
-    :doi:`10.1038/msb.2011.65`.
+    :doi:`10.1038/msb.2011.65`.

psamm/balancecheck.py

@@ -23,7 +23,7 @@
 
 from six.moves import reduce
 
-from .formula import Formula
+from .formula import Formula, ParseError
 
 logger = logging.getLogger(__name__)
 
@@ -106,10 +106,12 @@ def formula_balance(model):
             try:
                 f = Formula.parse(compound.formula).flattened()
                 compound_formula[compound.id] = f
-            except ValueError:
-                logger.warning(
-                    'Error parsing formula for compound {}: {}'.format(
-                        compound.id, compound.formula), exc_info=True)
+            except ParseError as e:
+                msg = 'Error parsing formula for compound {}:\n{}\n{}'.format(
+                    compound.id, e, compound.formula)
+                if e.indicator is not None:
+                    msg += '\n{}'.format(e.indicator)
+                logger.warning(msg)
 
     for reaction in model.reactions:
         yield reaction, reaction_formula(reaction.equation, compound_formula)

psamm/command.py

@@ -560,6 +560,7 @@ def main(command_class=None, args=None):
 
 
 def main_sbml(command_class=None, args=None):
+    """Run the SBML command line interface."""
     # Set up logging for the command line interface
     if 'PSAMM_DEBUG' in os.environ:
         level = getattr(logging, os.environ['PSAMM_DEBUG'].upper(), None)
@@ -575,16 +576,15 @@ def main_sbml(command_class=None, args=None):
             base_logger.addHandler(handler)
             base_logger.propagate = False
 
-    logger.warning(
-        'This command is experimental. It currently only fully parses level 3'
-        ' SBML files!')
-
     title = 'Metabolic modeling tools (SBML)'
     if command_class is not None:
         title, _, _ = command_class.__doc__.partition('\n\n')
 
     parser = argparse.ArgumentParser(description=title)
     parser.add_argument('model', metavar='file', help='SBML file')
+    parser.add_argument('--merge-compounds', action='store_true',
+                        help=('Merge identical compounds occuring in various'
+                              ' compartments.'))
     parser.add_argument(
         '-V', '--version', action='version',
         version='%(prog)s ' + package_version)
@@ -622,6 +622,9 @@ def main_sbml(command_class=None, args=None):
     context = FilePathContext(parsed_args.model)
     with context.open('r') as f:
         model = sbml.SBMLReader(f, context=context).create_model()
+        sbml.convert_sbml_model(model)
+        if parsed_args.merge_compounds:
+            sbml.merge_equivalent_compounds(model)
 
     # Instantiate command with model and run
     command = parsed_args.command(model, parsed_args)

psamm/commands/gapcheck.py

@@ -142,10 +142,11 @@ def run_sink_check(self, model, solver, threshold, implicit_sinks=True):
                 mass_balance_constrs[compound].delete()
 
             prob.set_objective(lhs)
-            result = prob.solve(lp.ObjectiveSense.Maximize)
-            if not result:
-                logger.warning('Failed to solve for compound: {}'.format(
-                    compound))
+            try:
+                result = prob.solve(lp.ObjectiveSense.Maximize)
+            except lp.SolverError as e:
+                logger.warning('Failed to solve for compound: {} ({})'.format(
+                    compound, e))
 
             if result.get_value(lhs) < threshold:
                 yield compound
@@ -190,11 +191,12 @@ def run_reaction_production_check(self, model, solver, threshold,
             prob.set_objective(v(reaction))
             for sense in (lp.ObjectiveSense.Maximize,
                           lp.ObjectiveSense.Minimize):
-                result = prob.solve(sense)
-                if not result:
+                try:
+                    result = prob.solve(sense)
+                except lp.SolverError as e:
                     self.fail(
                         'Failed to solve for compound, reaction: {}, {}:'
-                        ' {}'.format(compound, reaction, result.status))
+                        ' {}'.format(compound, reaction, e))
 
                 flux = result.get_value(v(reaction))
                 for compound, value in model.get_reaction_values(reaction):

psamm/commands/genedelete.py

@@ -16,17 +16,18 @@
 # Copyright 2014-2017  Jon Lund Steffensen <[email protected]>
 # Copyright 2015  Keith Dufault-Thompson <[email protected]>
 
-from __future__ import unicode_literals
+from __future__ import unicode_literals, division
 
 import time
 import logging
 
 from ..command import (Command, MetabolicMixin, ObjectiveMixin,
                        SolverCommandMixin, FilePrefixAppendAction)
 from .. import fluxanalysis
+from .. import moma
 from ..expression import boolean
 
-from six import string_types
+from six import string_types, text_type
 
 logger = logging.getLogger(__name__)
 
@@ -42,9 +43,16 @@ def init_parser(cls, parser):
         parser.add_argument(
             '--gene', metavar='genes', action=FilePrefixAppendAction,
             type=str, default=[], help='Delete multiple genes from model')
+        parser.add_argument(
+            '--method', metavar='method',
+            choices=['fba', 'lin_moma', 'lin_moma2', 'moma', 'moma2'],
+            type=text_type, default='fba',
+            help='Select which method to use. (fba, lin_moma, lin_moma2, '
+                 'moma, moma2)')
         super(GeneDeletionCommand, cls).init_parser(parser)
 
     def run(self):
+        """Delete the specified gene and solve using the desired method."""
         obj_reaction = self._get_objective()
 
         genes = set()
@@ -81,27 +89,60 @@ def run(self):
                     logger.info('Deleting reaction {}...'.format(reaction))
                     deleted_reactions.add(reaction)
 
-        solver = self._get_solver()
-        prob = fluxanalysis.FluxBalanceProblem(self._mm, solver)
+        if self._args.method in ['moma', 'moma2']:
+            solver = self._get_solver(quadratic=True)
+        else:
+            solver = self._get_solver()
 
-        try:
-            prob.maximize(obj_reaction)
-        except fluxanalysis.FluxBalanceError as e:
-            self.report_flux_balance_error(e)
-        wild = prob.get_flux(obj_reaction)
+        if self._args.method == 'fba':
+            logger.info('Solving using FBA...')
+            prob = fluxanalysis.FluxBalanceProblem(self._mm, solver)
+
+            try:
+                prob.maximize(obj_reaction)
+            except fluxanalysis.FluxBalanceError as e:
+                self.report_flux_balance_error(e)
 
-        for reaction in deleted_reactions:
-            flux_var = prob.get_flux_var(reaction)
-            prob.prob.add_linear_constraints(flux_var == 0)
+            wild = prob.get_flux(obj_reaction)
 
-        prob.maximize(obj_reaction)
-        deleteflux = prob.get_flux(obj_reaction)
+            for reaction in deleted_reactions:
+                flux_var = prob.get_flux_var(reaction)
+                prob.prob.add_linear_constraints(flux_var == 0)
+
+            prob.maximize(obj_reaction)
+            deleteflux = prob.get_flux(obj_reaction)
+        elif self._args.method in ['lin_moma', 'lin_moma2', 'moma', 'moma2']:
+            prob = moma.MOMAProblem(self._mm, solver)
+            wt_fluxes = prob.get_minimal_fba_flux(obj_reaction)
+            wild = wt_fluxes[obj_reaction]
+
+            for reaction in deleted_reactions:
+                flux_var = prob.get_flux_var(reaction)
+                prob.prob.add_linear_constraints(flux_var == 0)
+
+            try:
+                if self._args.method == 'moma':
+                    logger.info('Solving using MOMA...')
+                    prob.moma(wt_fluxes)
+                elif self._args.method == 'lin_moma':
+                    logger.info('Solving using linear MOMA...')
+                    prob.lin_moma(wt_fluxes)
+                elif self._args.method == 'moma2':
+                    logger.info('Solving using combined-model MOMA...')
+                    prob.moma2(obj_reaction, wild)
+                elif self._args.method == 'lin_moma2':
+                    logger.info('Solving using combined-model linear MOMA...')
+                    prob.lin_moma2(obj_reaction, wild)
+            except moma.MOMAError:
+                self.fail('Error computing the MOMA result.')
+
+            deleteflux = prob.get_flux(obj_reaction)
 
         logger.info(
             'Solving took {:.2f} seconds'.format(time.time() - start_time))
         logger.info(
             'Objective reaction after gene deletion has flux {}'.format(
-                abs(deleteflux) if deleteflux == 0 else deleteflux))
+                deleteflux + 0))
         if wild != 0:
             logger.info(
                 'Objective reaction has {:.2%} flux of wild type flux'.format(