Merge pull request #44 from biosustain/opt_gene

OptGene as a method

cameo/api/designer.py

@@ -14,8 +14,6 @@
 
 from __future__ import absolute_import, print_function
 
-__all__ = ['design']
-
 import re
 
 import numpy as np
@@ -29,8 +27,8 @@
 from cameo.api.hosts import hosts, Host
 from cameo.api.products import products
 from cameo.exceptions import SolveError
-from cameo.strain_design.heuristic import GeneKnockoutOptimization
-from cameo.strain_design.heuristic.objective_functions import biomass_product_coupled_yield
+from cameo.strain_design.deterministic import OptKnock
+from cameo.strain_design.heuristic import OptGene
 from cameo.ui import notice, searching, stop_loader
 from cameo.strain_design import pathway_prediction
 from cameo.util import TimeMachine
@@ -41,6 +39,10 @@
 
 import logging
 
+
+__all__ = ['design']
+
+
 logger = logging.getLogger(__name__)
 
 # TODO: implement cplex preference (if available)
@@ -51,11 +53,16 @@ def __call__(self, strategy, *args, **kwargs):
         (host, model, pathway) = (strategy[0], strategy[1], strategy[2])
         with TimeMachine() as tm:
             pathway.plug_model(model, tm)
-            objective = biomass_product_coupled_yield(model.biomass,
-                                                      pathway.product,
-                                                      model.carbon_source)
-            opt = GeneKnockoutOptimization(model=model, objective_function=objective, progress=True, plot=False)
-            return opt.run(product=pathway.product.id, max_evaluations=10000)
+            opt_gene = OptGene(model=model, plot=False)
+            opt_gene_designs = opt_gene.run(target=pathway.product.id, biomass=model.biomass,
+                                            substrate=model.carbon_source, max_evaluations=10000)
+
+            opt_knock = OptKnock(model=model)
+            opt_knock_designs = opt_knock.run(5, target=pathway.product.id, max_results=5)
+
+            designs = opt_gene_designs + opt_knock_designs
+
+            return designs
 
 
 class DesignerResult(Result):

cameo/core/pathway.py

@@ -120,7 +120,7 @@ def plug_model(self, model, tm=None):
             Optionally, a TimeMachine object can be added to the operation
 
         """
-        metabolites = reduce(lambda x, y: x+y, [r.metabolites.keys() for r in self.reactions], [])
+        metabolites = reduce(lambda x, y: x+y, [list(r.metabolites.keys()) for r in self.reactions], [])
         exchanges = [model.add_demand(m, prefix="EX_", time_machine=tm) for m in metabolites
                      if m not in model.metabolites]
         for exchange in exchanges:

cameo/core/result.py

@@ -61,5 +61,5 @@ def meta_information(self):
     def data_frame(self):
         raise NotImplementedError
 
-    def plot(self, grid=None, width=None, height=None, title=None):
+    def plot(self, grid=None, width=None, height=None, title=None, *args, **kwargs):
         raise NotImplementedError

cameo/core/solver_based_model.py

@@ -19,7 +19,8 @@
 from __future__ import absolute_import, print_function
 from functools import partial
 
-__all__ = ['to_solver_based_model', 'SolverBasedModel']
+from cobra.core import Metabolite
+
 
 import six
 
@@ -47,6 +48,10 @@
 
 import logging
 
+
+__all__ = ['to_solver_based_model', 'SolverBasedModel']
+
+
 logger = logging.getLogger(__name__)
 
 add = Add._from_args
@@ -149,7 +154,7 @@ def objective(self):
 
     @objective.setter
     def objective(self, value):
-        if isinstance(value, str):
+        if isinstance(value, six.string_types):
             value = self.reactions.get_by_id(value)
         if isinstance(value, Reaction):
             self.solver.objective = self.solver.interface.Objective(value.flux_expression, sloppy=True)
@@ -183,7 +188,7 @@ def solver(self, value):
         not_valid_interface = ValueError(
             '%s is not a valid solver interface. Pick from %s, or specify an optlang interface (e.g. optlang.glpk_interface).' % (
                 value, list(config.solvers.keys())))
-        if isinstance(value, str):
+        if isinstance(value, six.string_types):
             try:
                 interface = config.solvers[value]
             except KeyError:
@@ -309,12 +314,12 @@ def add_demand(self, metabolite, prefix="DM_", time_machine=None):
         demand_reaction.upper_bound = 1000
         if time_machine is not None:
             time_machine(do=partial(self.add_reactions, [demand_reaction]),
-                         undo=partial(self.remove_reactions, [demand_reaction]))
+                         undo=partial(self.remove_reactions, [demand_reaction], delete=False))
         else:
             self.add_reactions([demand_reaction])
         return demand_reaction
 
-    def fix_objective_as_constraint(self, time_machine=None):
+    def fix_objective_as_constraint(self, time_machine=None, fraction=1):
         """Fix current objective as an additional constraint (e.g., ..math`c^T v >= max c^T v`).
 
         Parameters
@@ -326,8 +331,8 @@ def fix_objective_as_constraint(self, time_machine=None):
         -------
         None
         """
-        objective_value = self.solve().objective_value
-        constraint = self.solver.interface.Constraint(self.objective.expression, lb=objective_value,
+        objective_value = self.solve().objective_value * fraction
+        constraint = self.solver.interface.Constraint(self.objective.expression,
                                                       name='Fixed_objective_{}'.format(self.objective.name))
         if self.objective.direction == 'max':
             constraint.lb = objective_value
@@ -543,7 +548,7 @@ def load_medium(self, medium, copy=False):
             model._load_medium_from_dict(model, medium)
         elif isinstance(medium, pandas.DataFrame):
             model._load_medium_from_dataframe(model, medium)
-        elif isinstance(medium, str):
+        elif isinstance(medium, six.string_types):
             model._load_medium_from_file(model, medium)
         else:
             raise AssertionError("input type (%s) is not valid" % type(medium))
@@ -554,14 +559,38 @@ def _ids_to_reactions(self, reactions):
         """Translate reaction IDs into reactions (skips reactions)."""
         clean_reactions = list()
         for reaction in reactions:
-            if isinstance(reaction, str):
+            if isinstance(reaction, six.string_types):
                 clean_reactions.append(self.reactions.get_by_id(reaction))
             elif isinstance(reaction, Reaction):
                 clean_reactions.append(reaction)
             else:
                 raise Exception('%s is not a reaction or reaction ID.' % reaction)
         return clean_reactions
 
+    def reaction_for(self, value, time_machine=None):
+        if isinstance(value, Reaction):
+            value = self.reactions.get_by_id(value.id)
+
+        if isinstance(value, six.string_types):
+            try:
+                value = self.reactions.get_by_id(value)
+            except KeyError:
+                try:
+                    value = self.metabolites.get_by_id(value)
+                except KeyError:
+                    raise KeyError("Invalid target %s." % value)
+
+        if isinstance(value, cobra.core.Metabolite):
+            try:
+                value = self.reactions.get_by_id("EX_%s" % value.id)
+            except KeyError:
+                try:
+                    value = self.reactions.get_by_id("DM_%s" % value.id)
+                except KeyError:
+                    value = self.add_demand(value, time_machine=time_machine)
+
+        return value
+
     @staticmethod
     def _load_medium_from_dict(model, medium):
         for rid, values in six.iteritems(medium):

cameo/flux_analysis/analysis.py

@@ -148,7 +148,7 @@ def phenotypic_phase_plane(model, variables=[], objective=None, points=20, view=
     References
     ----------
     [1] Edwards, J. S., Ramakrishna, R. and Palsson, B. O. (2002). Characterizing the metabolic phenotype: a phenotype
-    phase plane analysis. Biotechnology and Bioengineering, 77(1), 27–36. doi:10.1002/bit.10047
+        phase plane analysis. Biotechnology and Bioengineering, 77(1), 27–36. doi:10.1002/bit.10047
     """
     if isinstance(variables, str):
         variables = [variables]
@@ -159,11 +159,11 @@ def phenotypic_phase_plane(model, variables=[], objective=None, points=20, view=
         view = config.default_view
     with TimeMachine() as tm:
         if objective is not None:
-            if isinstance(objective, Metabolite):
-                try:
-                    objective = model.reactions.get_by_id("DM_%s" % objective.id)
-                except KeyError:
-                    objective = model.add_demand(objective, time_machine=tm)
+            try:
+                objective = model.reaction_for(objective, time_machine=tm)
+            except KeyError:
+                pass
+
             tm(do=partial(setattr, model, 'objective', objective),
                undo=partial(setattr, model, 'objective', model.objective))
 
@@ -416,48 +416,53 @@ def flux_balance_impact_degree(model, knockouts, view=config.default_view, metho
 
     Returns
     -------
-    int: perturbation
-        The number of changes in reachable reactions (reactions that can carry flux)
+    FluxBalanceImpactDegreeResult: perturbation
+        The changes in reachable reactions (reactions that can carry flux)
     """
 
     if method == "fva":
-        _fbid_fva(model, knockouts, view)
+        reachable_reactions, perturbed_reactions =_fbid_fva(model, knockouts, view)
     elif method == "em":
         raise NotImplementedError("Elementary modes approach is not implemented")
     else:
         raise ValueError("%s method is not valid to compute Flux Balance Impact Degree" % method)
 
+    return FluxBalanceImpactDegreeResult(reachable_reactions, perturbed_reactions, method)
+
 
 def _fbid_fva(model, knockouts, view):
-    tm = TimeMachine()
-    for reaction in model.reactions:
-        if reaction.reversibility:
-            tm(do=partial(setattr, reaction, 'lower_bound', -1),
-               undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
-            tm(do=partial(setattr, reaction, 'upper_bound', 1),
-               undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
-        else:
-            tm(do=partial(setattr, reaction, 'lower_bound', 0),
-               undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
-            tm(do=partial(setattr, reaction, 'upper_bound', 1),
-               undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
+    with TimeMachine() as tm:
+
+        for reaction in model.reactions:
+            if reaction.reversibility:
+                tm(do=partial(setattr, reaction, 'lower_bound', -1),
+                   undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
+                tm(do=partial(setattr, reaction, 'upper_bound', 1),
+                   undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
+            else:
+                tm(do=partial(setattr, reaction, 'lower_bound', 0),
+                   undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
+                tm(do=partial(setattr, reaction, 'upper_bound', 1),
+                   undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
+
+        wt_fva = flux_variability_analysis(model, view=view, remove_cycles=False)
+        wt_fva._data_frame = wt_fva._data_frame.apply(numpy.round)
+
+        reachable_reactions = wt_fva.data_frame.query("lower_bound != 0 | upper_bound != 0")
 
-    wt_fva = flux_variability_analysis(model, view)
-    for reaction in knockouts:
-        tm(do=partial(setattr, reaction, 'upper_bound', 0),
-           undo=partial(setattr, reaction, 'upper_bound', reaction.upper_bound))
-        tm(do=partial(setattr, reaction, 'lower_bound', 0),
-           undo=partial(setattr, reaction, 'lower_bound', reaction.lower_bound))
+        for reaction in model._ids_to_reactions(knockouts):
+            reaction.knock_out(tm)
 
-    mt_fva = flux_variability_analysis(model, view)
+        mt_fva = flux_variability_analysis(model, reactions=reachable_reactions.index, view=view, remove_cycles=False)
+        mt_fva._data_frame = mt_fva._data_frame.apply(numpy.round)
 
-    perturbation = 0
-    for reaction in model.reactions:
-        if wt_fva[reaction.id] != 0 and mt_fva[reaction.id] == 0:
-            perturbation += 1
+        perturbed_reactions = []
+        for reaction in reachable_reactions.index:
+            if wt_fva.upper_bound(reaction) != mt_fva.upper_bound(reaction) or \
+              wt_fva.lower_bound(reaction) != wt_fva.lower_bound(reaction):
+                perturbed_reactions.append(reaction)
 
-    tm.reset()
-    return perturbation
+        return list(reachable_reactions.index), perturbed_reactions
 
 
 class PhenotypicPhasePlaneResult(Result):
@@ -471,12 +476,12 @@ def __init__(self, phase_plane, variable_ids, objective, *args, **kwargs):
     def data_frame(self):
         return pandas.DataFrame(self._phase_plane)
 
-    def plot(self, grid=None, width=None, height=None, title=None, axis_font_size=None, **kwargs):
+    def plot(self, grid=None, width=None, height=None, title=None, axis_font_size=None, color="lightblue", **kwargs):
         if len(self.variable_ids) > 1:
             notice("Multi-dimensional plotting is not supported")
             return
         plotting.plot_production_envelope(self._phase_plane, objective=self.objective, key=self.variable_ids[0],
-                                          grid=grid, width=width, height=height, title=title,
+                                          grid=grid, width=width, height=height, title=title, color=color,
                                           axis_font_size=axis_font_size, **kwargs)
 
     def __getitem__(self, item):
@@ -507,15 +512,64 @@ def __init__(self, data_frame, *args, **kwargs):
     def data_frame(self):
         return self._data_frame
 
-    def plot(self, index=None, grid=None, width=None, height=None, title=None, axis_font_size=None):
+    def plot(self, index=None, grid=None, width=None, height=None, title=None, axis_font_size=None, color="lightblue",
+             **kwargs):
         if index is None:
             index = self.data_frame.index[0:10]
         fva_result = self.data_frame.loc[index]
         plotting.plot_flux_variability_analysis(fva_result, grid=grid, width=width, height=height, title=title,
-                                                axis_font_size=axis_font_size)
+                                                axis_font_size=axis_font_size, color=color)
 
     def __getitem__(self, item):
         return self._data_frame[item]
 
+    def upper_bound(self, item):
+        if isinstance(item, Reaction):
+            item = item.id
+        return self['upper_bound'][item]
+
+    def lower_bound(self, item):
+        if isinstance(item, Reaction):
+            item = item.id
+        return self['lower_bound'][item]
+
     def iterrows(self):
         return self._data_frame.iterrows()
+
+
+class FluxBalanceImpactDegreeResult(Result):
+
+    def __init__(self, reachable_reactions, perturbed_reactions, method, *args, **kwargs):
+        super(FluxBalanceImpactDegreeResult, self).__init__(*args, **kwargs)
+        self._method = method
+        self._reachable_reactions = reachable_reactions
+        self._perturbed_reactions = perturbed_reactions
+
+    def __contains__(self, item):
+        if isinstance(item, Reaction):
+            item = item.id
+        return item in self._reachable_reactions
+
+    def _repr_html_(self):
+        return """
+        <h3>Flux Balance Impact Degree</h3>
+        <ul>
+            <li> Degree: %i</li>
+            <li> Reactions: %i</li>
+        </ul>
+        %s
+        """ % (self.degree, len(self._reachable_reactions), self.data_frame._repr_html_())
+
+    @property
+    def degree(self):
+        return len(self._perturbed_reactions)
+
+    @property
+    def data_frame(self):
+        data_frame = pandas.DataFrame(columns=["perturbed"])
+        for reaction in self._reachable_reactions:
+            data_frame.loc[reaction] = [reaction in self._perturbed_reactions]
+        return data_frame
+
+    def plot(self, grid=None, width=None, height=None, title=None):
+        pass