Merge branch 'flux-coupling' into develop

docs/api.rst

@@ -16,6 +16,7 @@ PSAMM API
    expression_boolean
    fastcore
    fluxanalysis
+   fluxcoupling
    formula
    gapfill
    lpsolver_lp

docs/commands.rst

@@ -169,6 +169,21 @@ 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.
 
+Flux coupling analysis (``fluxcoupling``)
+-----------------------------------------
+
+The flux coupling analysis 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).
+
+.. code-block:: shell
+
+    $ psamm-model fluxcoupling
+
 Stoichiometric consistency check (``masscheck``)
 ------------------------------------------------
 

docs/fluxcoupling.rst

@@ -0,0 +1,6 @@
+
+``psamm.fluxcoupling`` -- Flux coupling analysis
+================================================
+
+.. automodule:: psamm.fluxcoupling
+   :members:

docs/references.rst

@@ -5,6 +5,9 @@ References
 .. [Steffensen15] Steffensen JL, Dufault-Thompson K, Zhang Y. PSAMM: A Portable
     System for the Analysis of Metabolic Models. Submitted.
 
+.. [Burgard04] Burgard AP, Nikolaev E V, Schilling CH, Maranas CD. Flux
+    coupling analysis of genome-scale metabolic network reconstructions.
+    Genome Res. 2004;14: 301–312. :doi:`10.1101/gr.1926504`.
 .. [Gevorgyan08] Gevorgyan A, Poolman MG, Fell DA. Detection of stoichiometric
     inconsistencies in biomolecular models. Bioinformatics. 2008;24: 2245–2251.
     :doi:`10.1093/bioinformatics/btn425`.

psamm/command.py

@@ -45,7 +45,7 @@
 from .reaction import Compound
 from .datasource.native import NativeModel
 from .datasource import sbml
-from . import fluxanalysis, massconsistency, fastcore
+from . import fluxanalysis, fluxcoupling, massconsistency, fastcore
 from .lpsolver import generic
 
 from six import add_metaclass, iteritems
@@ -553,6 +553,121 @@ def run(self):
                         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."""
 

psamm/fluxcoupling.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 <http://www.gnu.org/licenses/>.
+#
+# Copyright 2015  Jon Lund Steffensen <jon_steffensen@uri.edu>
+
+"""Flux coupling analysis
+
+Described in [Burgard04]_.
+"""
+
+import logging
+
+from six import iteritems
+
+from psamm.lpsolver import lp
+
+
+logger = logging.getLogger(__name__)
+
+
+class CouplingClass(object):
+    """Enumeration of coupling types."""
+
+    Inconsistent = object()
+    """Reaction is flux inconsistent (v1 is always zero)."""
+
+    Uncoupled = object()
+    """Uncoupled reactions."""
+
+    DirectionalForward = object()  # v1 -> v2
+    """Directionally coupled from reaction 1 to 2."""
+
+    DirectionalReverse = object()  # v2 -> v1
+    """Directionally coupled from reaction 2 to 1."""
+
+    Partial = object()  # v1 <-> v2
+    """Partially coupled reactions."""
+
+    Full = object()  # v1 <=> v2
+    """Fully coupled reactions."""
+
+
+class FluxCouplingProblem(object):
+    """A specific flux coupling analysis applied to a metabolic model.
+
+    Args:
+        model: MetabolicModel to apply the flux coupling problem to
+        bounded: Dictionary of reactions with minimum flux values
+        solver: LP solver instance to use.
+    """
+
+    def __init__(self, model, bounded, solver):
+        self._prob = solver.create_problem()
+
+        # Define t variable
+        self._prob.define('t')
+        t = self._prob.var('t')
+
+        # Define flux variables
+        for reaction_id in model.reactions:
+            self._prob.define(('vbow', reaction_id))
+
+            lower, upper = model.limits[reaction_id]
+            if reaction_id in bounded:
+                lower = bounded[reaction_id]
+            flux_bow = self._prob.var(('vbow', reaction_id))
+            self._prob.add_linear_constraints(
+                flux_bow >= t * lower, flux_bow <= t * upper)
+
+        # Define mass balance constraints
+        massbalance_lhs = {compound: 0 for compound in model.compounds}
+        for (compound, reaction_id), value in iteritems(model.matrix):
+            flux_bow = self._prob.var(('vbow', reaction_id))
+            massbalance_lhs[compound] += flux_bow * value
+        for compound, lhs in iteritems(massbalance_lhs):
+            self._prob.add_linear_constraints(lhs == 0)
+
+        self._reaction_constr = None
+
+    def solve(self, reaction_1, reaction_2):
+        """Return the flux coupling between two reactions
+
+        The flux coupling is returned as a tuple indicating the minimum and
+        maximum value of the v1/v2 reaction flux ratio. A value of None as
+        either the minimum or maximum indicates that the interval is unbounded
+        in that direction.
+        """
+        # Update objective for reaction_1
+        self._prob.set_linear_objective(self._prob.var(('vbow', reaction_1)))
+
+        # Update constraint for reaction_2
+        if self._reaction_constr is not None:
+            self._reaction_constr.delete()
+
+        reaction_2_vbow = self._prob.var(('vbow', reaction_2))
+        self._reaction_constr, = self._prob.add_linear_constraints(
+            reaction_2_vbow == 1)
+
+        results = []
+        for sense in (lp.ObjectiveSense.Minimize, lp.ObjectiveSense.Maximize):
+            result = self._prob.solve(sense)
+            if not result:
+                results.append(None)
+            else:
+                results.append(result.get_value(('vbow', reaction_1)))
+
+        return tuple(results)
+
+
+def classify_coupling(coupling):
+    """Return a constant indicating the type of coupling.
+
+    Depending on the type of coupling, one of the constants from
+    :class:`.CouplingClass` is returned.
+
+    Args:
+        coupling: Tuple of minimum and maximum flux ratio
+    """
+    lower, upper = coupling
+
+    if lower is None and upper is None:
+        return CouplingClass.Uncoupled
+    elif lower is None or upper is None:
+        return CouplingClass.DirectionalReverse
+    elif lower == 0.0 and upper == 0.0:
+        return CouplingClass.Inconsistent
+    elif lower <= 0.0 and upper >= 0.0:
+        return CouplingClass.DirectionalForward
+    elif abs(lower - upper) < 1e-6:
+        return CouplingClass.Full
+    else:
+        return CouplingClass.Partial