Merge pull request #30 from KristianJensen/devel

Better FSEOF output and more tests

cameo/core/reaction.py

@@ -185,7 +185,7 @@ def lower_bound(self, value):
 
             elif self._upper_bound <= 0: # reverse irreversible
                 if value > 0:
-                    reverse_variable.ub = 0
+                    reverse_variable.lb = 0
                     reverse_variable.ub = 0
                     forward_variable.ub = value
                     self._upper_bound = value

cameo/network_analysis/util.py

@@ -22,9 +22,9 @@ def distance_based_on_molecular_formula(metabolite1, metabolite2, normalize=True
     for element in elements:
          distance += abs(metabolite1.formula.elements.get(element, 0) - metabolite2.formula.elements.get(element, 0))
     if normalize:
-        try:
-            return distance / sum(list(metabolite1.formula.elements.values()) + list(metabolite2.formula.elements.values()))
-        except:
-            print(metabolite1, metabolite2, metabolite1.formula.elements, metabolite2.formula.elements)
+        #try:
+        return distance / sum(list(metabolite1.formula.elements.values()) + list(metabolite2.formula.elements.values()))
+        #except:
+        #    print(metabolite1, metabolite2, metabolite1.formula.elements, metabolite2.formula.elements)
     else:
         return distance

cameo/strain_design/deterministic/flux_variability_based.py

@@ -24,14 +24,15 @@
 from pandas import DataFrame, pandas
 from progressbar import ProgressBar
 from progressbar.widgets import ETA, Bar
-from copy import copy
 
+from cameo.flux_analysis import fba
 from cameo import config, flux_variability_analysis, Reaction
 from cameo.parallel import SequentialView, MultiprocessingView
 from cameo.core.solver_based_model import Reaction
 from cameo.strain_design import StrainDesignMethod
 from cameo.flux_analysis.analysis import phenotypic_phase_plane
 from cameo.util import TimeMachine
+import cameo
 
 import logging
 import six
@@ -193,8 +194,8 @@ def run(self, surface_only=False, improvements_only=True, view=None):
                 tm(do=int, undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
                 tm(do=int, undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
             target_reaction = self.design_space_model.reactions.get_by_id(self.objective)
-            tm(do=int, undo=partial(setattr, target_reaction, 'lower_bound', reaction.lower_bound))
-            tm(do=int, undo=partial(setattr, target_reaction, 'upper_bound', reaction.upper_bound))
+            tm(do=int, undo=partial(setattr, target_reaction, 'lower_bound', target_reaction.lower_bound))
+            tm(do=int, undo=partial(setattr, target_reaction, 'upper_bound', target_reaction.upper_bound))
 
             if view is None:
                 view = config.default_view
@@ -271,7 +272,7 @@ def _set_bounds(self, point):
         target_reaction.lower_bound, target_reaction.upper_bound = target_bound, target_bound
 
 
-def fseof(model, enforced_reaction, max_enforced_flux=0.9, granularity=10, primary_objective=None, exclude=[]):
+def fseof(model, enforced_reaction, max_enforced_flux=0.9, granularity=10, primary_objective=None, solution_method=fba, exclude=[]):
     """
     Performs a Flux Scanning based on Enforced Objective Flux (FSEOF) analysis.
     :param model: SolverBasedModel
@@ -282,37 +283,36 @@ def fseof(model, enforced_reaction, max_enforced_flux=0.9, granularity=10, prima
     :param exclude: Iterable of reactions or reaction ids that will not be included in the output.
     :return: List of reactions that correlate with enforced flux.
     """
-
-    model = model.copy()
     ndecimals = config.ndecimals
+    with TimeMachine() as tm:
+
+        # Convert enforced reaction to Reaction object
+        if not isinstance(enforced_reaction, Reaction):
+            enforced_reaction = model.reactions.get_by_id(enforced_reaction)
+        primary_objective = primary_objective or model.objective
+
+        # Exclude list
+        exclude += model.exchanges
+        exclude_ids = [enforced_reaction.id]
+        for reaction in exclude:
+            if isinstance(reaction, Reaction):
+                exclude_ids.append(reaction.id)
+            else:
+                exclude_ids.append(reaction)
 
-    # Convert enforced reaction to Reaction object
-    if not isinstance(enforced_reaction, Reaction):
-        enforced_reaction = model.reactions.get_by_id(enforced_reaction)
-    primary_objective = primary_objective or copy(model.objective)
-
-    # Exclude list
-    exclude_ids = []
-    for reaction in exclude:
-        if isinstance(reaction, Reaction):
-            exclude_ids.append(reaction.id)
-        else:
-            exclude_ids.append(reaction)
-
-    original_objective = copy(model.objective)
-    original_lb = enforced_reaction.lower_bound
-    original_ub = enforced_reaction.upper_bound
+        tm(do=int, undo=partial(setattr, model, "objective", model.objective))
+        tm(do=int, undo=partial(setattr, enforced_reaction, "lower_bound", enforced_reaction.lower_bound))
+        tm(do=int, undo=partial(setattr, enforced_reaction, "upper_bound", enforced_reaction.upper_bound))
 
-    try:
         # Find initial flux of enforced reaction
         model.objective = primary_objective
-        initial_solution = model.solve()
+        initial_solution = solution_method(model)
         initial_fluxes = initial_solution.fluxes
         initial_flux = round(initial_fluxes[enforced_reaction.id], ndecimals)
 
         # Find theoretical maximum of enforced reaction
         model.objective = enforced_reaction
-        max_theoretical_flux = round(model.solve().fluxes[enforced_reaction.id], ndecimals)
+        max_theoretical_flux = round(solution_method(model).fluxes[enforced_reaction.id], ndecimals)
 
         max_flux = max_theoretical_flux * max_enforced_flux
 
@@ -327,23 +327,64 @@ def fseof(model, enforced_reaction, max_enforced_flux=0.9, granularity=10, prima
         for enforcement in enforcements:
             enforced_reaction.lower_bound = enforcement
             enforced_reaction.upper_bound = enforcement
-            solution = model.solve()
+            solution = solution_method(model)
             for reaction_id, flux in solution.fluxes.items():
                 results[reaction_id].append(round(flux, config.ndecimals))
 
-    finally:
-        # Reset objective and bounds
-        enforced_reaction.lower_bound = original_lb
-        enforced_reaction.upper_bound = original_ub
-        model.objective = original_objective
-
     # Test each reaction
     fseof_reactions = []
     for reaction_id, fluxes in results.items():
-        if abs(fluxes[-1]) > abs(fluxes[0]) and min(fluxes)*max(fluxes) >= 0:
-            fseof_reactions.append(reaction_id)
+        if reaction_id not in exclude_ids and abs(fluxes[-1]) > abs(fluxes[0]) and min(fluxes)*max(fluxes) >= 0:
+            fseof_reactions.append(model.reactions.get_by_id(reaction_id))
+
+    return FseofResult(fseof_reactions, enforced_reaction, model)
+
 
-    return fseof_reactions
+class FseofResult(cameo.core.result.Result):
+    """
+    Object for holding FSEOF results.
+    """
+    def __init__(self, reactions, objective, model, *args, **kwargs):
+        super(FseofResult, self).__init__(*args, **kwargs)
+        self._reactions = reactions
+        self._objective = objective
+        self._model = model
+
+    def __iter__(self):
+        return iter(self.reactions)
+
+    def __eq__(self, other):
+        return isinstance(other, self.__class__) and self.objective == other.objective and self.reactions == other.reactions
+
+    @property
+    def reactions(self):
+        return self._reactions
+
+    @property
+    def model(self):
+        return self._model
+
+    @property
+    def objective(self):
+        return self._objective
+
+    def _repr_html_(self):
+        template = """
+<table>
+     <tr>
+        <td><b>Enforced objective</b></td>
+        <td>%(objective)s</td>
+    </tr>
+    <tr>
+        <td><b>Reactions</b></td>
+        <td>%(reactions)s</td>
+    <tr>
+</table>"""
+        return template % {'objective': self.objective.nice_id, 'reactions': "<br>".join(reaction.id for reaction in self.reactions)}
+
+    @property
+    def data_frame(self):
+        return pandas.DataFrame((r.id for r in self.reactions), columns=["Reaction id"])
 
 
 if __name__ == '__main__':

tests/test_parallel.py

@@ -20,6 +20,7 @@
 from cameo.parallel import SequentialView
 import subprocess
 from time import sleep
+from multiprocessing import cpu_count
 
 from IPython.parallel import Client, interactive
 from six.moves import range
@@ -62,6 +63,15 @@ def test_apply(self):
             for i in range(100):
                 self.assertEqual(self.view.apply(to_the_power_of_2, i), SOLUTION[i])
 
+        def test_length(self):
+            self.assertEqual(len(self.view), cpu_count())
+
+            view = MultiprocessingView(4)
+            self.assertEqual(len(view), 4)
+
+            view = MultiprocessingView(processes=3)
+            self.assertEqual(len(view), 3)
+
 except ImportError:
     print("Skipping MultiprocessingView tests ...")
 

tests/test_strain_design_deterministic.py

@@ -18,10 +18,9 @@
 import unittest
 
 from cameo import load_model
-from cameo.util import RandomGenerator as Random
-from cameo.strain_design.deterministic import fseof
-from cameo.parallel import SequentialView, MultiprocessingView
+from cameo.strain_design.deterministic.flux_variability_based import fseof, FseofResult
 from six.moves import range
+from pandas import DataFrame
 
 TESTDIR = os.path.dirname(__file__)
 iJO_MODEL = load_model(os.path.join(TESTDIR, 'data/iJO1366.xml'), sanitize=False)
@@ -33,8 +32,17 @@ def setUp(self):
         self.model.solver = 'glpk'
 
     def test_fseof(self):
+        objective = self.model.objective
         fseof_result = fseof(self.model, enforced_reaction="EX_succ_LPAREN_e_RPAREN_")
-        self.assertIsInstance(fseof_result, list)
+        self.assertIsInstance(fseof_result, FseofResult)
+        self.assertIs(objective, self.model.objective)
+
+    def test_fseof_result(self):
+        fseof_result = fseof(self.model, self.model.reactions.EX_ac_LPAREN_e_RPAREN_, 0.8, exclude=["PGI"])
+        self.assertIsInstance(fseof_result.data_frame, DataFrame)
+        self.assertIs(fseof_result.objective, self.model.reactions.EX_ac_LPAREN_e_RPAREN_)
+        self.assertIs(fseof_result.model, self.model)
+        self.assertEqual(list(fseof_result), list(fseof_result.reactions))
 
 
 if __name__ == "__main__":

tests/test_util.py

@@ -16,10 +16,13 @@
 
 import unittest
 from functools import partial
+from itertools import chain
 
 from cameo.util import TimeMachine, generate_colors, Singleton, partition
+from cameo.network_analysis.util import distance_based_on_molecular_formula
 import six
 from six.moves import range
+from cobra import Metabolite
 
 
 class TimeMachineTestCase(unittest.TestCase):
@@ -71,18 +74,34 @@ def test_partition(self):
             test_output = partition(fixture, chunks)
             self.assertEqual(len(fixture), sum(map(len, test_output)))
             self.assertEqual(len(test_output), chunks)
+            self.assertEqual(list(fixture), list(chain(*test_output)))
             for out_chunk in test_output:
                 self.assertTrue(set(out_chunk).issubset(set(fixture)))
 
         bad_input = 5
         self.assertRaises(TypeError, partition, bad_input, chunks)
 
+    @unittest.skip  # Development API changes to cobra.core.Metabolite
+    def test_distance_based_on_molecular_formula(self):  # from network_analysis.util
+        met1 = Metabolite("H2O", formula="H2O")
+        met2 = Metabolite("H2O2", formula="H2O2")
+        met3 = Metabolite("C6H12O6", formula="C6H12O6")
+
+        self.assertEqual(distance_based_on_molecular_formula(met1, met2, normalize=False), 1)
+        self.assertEqual(distance_based_on_molecular_formula(met1, met2, normalize=True), 1./7)
+
+        self.assertEqual(distance_based_on_molecular_formula(met2, met3, normalize=False), 20)
+        self.assertEqual(distance_based_on_molecular_formula(met2, met3, normalize=True), 20./28)
+
+        self.assertEqual(distance_based_on_molecular_formula(met1, met3, normalize=False), 21)
+        self.assertEqual(distance_based_on_molecular_formula(met1, met3, normalize=True), 21./27)
+
 
 class TestSingleton(unittest.TestCase):
     def test_singleton(self):
         s1 = Singleton()
         s2 = Singleton()
-        self.assertEqual(s1, s2)
+        self.assertIs(s1, s2)
 
 
 if __name__ == "__main__":