diff --git a/.gitignore b/.gitignore index 97feb8f..2df6bcf 100644 --- a/.gitignore +++ b/.gitignore @@ -1,9 +1,13 @@ # Pycharm folders .idea/* +.idea/ cobra_docker/work/.idea/* *___jb_old___ *___jb_tmp___ +# Docker temp files +docker/temp/ + ########################################### ### Standard Python .gitignore ### ########################################### diff --git a/doc/Makefile b/doc/Makefile index 22a223c..89ad0f1 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -3,10 +3,11 @@ # You can set these variables from the command line. SPHINXOPTS = -SPHINXBUILD = python -msphinx +SPHINXBUILD = python3 -msphinx SPHINXPROJ = pytfa SOURCEDIR = . BUILDDIR = _build +AUTODIR = _autogen # Put it first so that "make" without argument is like "make help". help: @@ -17,4 +18,4 @@ help: # Catch-all target: route all unknown targets to Sphinx using the new # "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). %: Makefile - @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/doc/conf.py b/doc/conf.py index 3055f8d..dc8f954 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -21,7 +21,34 @@ import sys sys.path.insert(0, os.path.abspath('..')) +# In order to build documentation that requires libraries to import +class Mock(object): + def __init__(self, *args, **kwargs): + return + + def __call__(self, *args, **kwargs): + return Mock() + + @classmethod + def __getattr__(cls, name): + if name in ('__file__', '__path__'): + return '/dev/null' + else: + return Mock() + + +# These modules should correspond to the importable Python packages. +MOCK_MODULES = [ + 'cobra', + 'numpy', + 'scipy', + 'pandas', + 'optlang', 'optlang.interface', 'optlang.symbolics', + 'scipy', 'scipy.io', +] +for mod_name in MOCK_MODULES: + sys.modules[mod_name] = Mock() # -- General configuration ------------------------------------------------ # If your documentation needs a minimal Sphinx version, state it here. @@ -35,7 +62,13 @@ 'sphinx.ext.intersphinx', 'sphinx.ext.coverage', 'sphinx.ext.viewcode', - 'sphinx.ext.mathjax'] + 'sphinx.ext.mathjax', + 'autoapi.extension'] + + +# Document Python Code +autoapi_type = 'python' +autoapi_dirs = ['../pytfa'] # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -59,10 +92,10 @@ # built documents. # # The short X.Y version. -version = '0.8' +version = '0.9' # The full version, including alpha/beta/rc tags. -release = '0.8.1' +release = '0.9.0-b0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/pytfa.analysis.rst b/doc/pytfa.analysis.rst deleted file mode 100644 index 004f92f..0000000 --- a/doc/pytfa.analysis.rst +++ /dev/null @@ -1,38 +0,0 @@ -pytfa.analysis package -====================== - -Submodules ----------- - -pytfa.analysis.manipulation module ----------------------------------- - -.. automodule:: pytfa.analysis.manipulation - :members: - :undoc-members: - :show-inheritance: - -pytfa.analysis.sampling module ------------------------------- - -.. automodule:: pytfa.analysis.sampling - :members: - :undoc-members: - :show-inheritance: - -pytfa.analysis.variability module ---------------------------------- - -.. automodule:: pytfa.analysis.variability - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: pytfa.analysis - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/pytfa.core.rst b/doc/pytfa.core.rst deleted file mode 100644 index a6eba7b..0000000 --- a/doc/pytfa.core.rst +++ /dev/null @@ -1,22 +0,0 @@ -pytfa.core package -================== - -Submodules ----------- - -pytfa.core.model module ------------------------ - -.. automodule:: pytfa.core.model - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: pytfa.core - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/pytfa.io.rst b/doc/pytfa.io.rst deleted file mode 100644 index 0fcfbd9..0000000 --- a/doc/pytfa.io.rst +++ /dev/null @@ -1,62 +0,0 @@ -pytfa.io package -================ - -Submodules ----------- - -pytfa.io.base module --------------------- - -.. automodule:: pytfa.io.base - :members: - :undoc-members: - :show-inheritance: - -pytfa.io.dict module --------------------- - -.. automodule:: pytfa.io.dict - :members: - :undoc-members: - :show-inheritance: - -pytfa.io.enrichment module --------------------------- - -.. automodule:: pytfa.io.enrichment - :members: - :undoc-members: - :show-inheritance: - -pytfa.io.json module --------------------- - -.. automodule:: pytfa.io.json - :members: - :undoc-members: - :show-inheritance: - -pytfa.io.plotting module ------------------------- - -.. automodule:: pytfa.io.plotting - :members: - :undoc-members: - :show-inheritance: - -pytfa.io.viz module -------------------- - -.. automodule:: pytfa.io.viz - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: pytfa.io - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/pytfa.optim.rst b/doc/pytfa.optim.rst deleted file mode 100644 index 18fead7..0000000 --- a/doc/pytfa.optim.rst +++ /dev/null @@ -1,62 +0,0 @@ -pytfa.optim package -=================== - -Submodules ----------- - -pytfa.optim.constraints module ------------------------------- - -.. automodule:: pytfa.optim.constraints - :members: - :undoc-members: - :show-inheritance: - -pytfa.optim.debugging module ----------------------------- - -.. automodule:: pytfa.optim.debugging - :members: - :undoc-members: - :show-inheritance: - -pytfa.optim.reformulation module --------------------------------- - -.. automodule:: pytfa.optim.reformulation - :members: - :undoc-members: - :show-inheritance: - -pytfa.optim.relaxation module ------------------------------ - -.. automodule:: pytfa.optim.relaxation - :members: - :undoc-members: - :show-inheritance: - -pytfa.optim.utils module ------------------------- - -.. automodule:: pytfa.optim.utils - :members: - :undoc-members: - :show-inheritance: - -pytfa.optim.variables module ----------------------------- - -.. automodule:: pytfa.optim.variables - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: pytfa.optim - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/pytfa.rst b/doc/pytfa.rst deleted file mode 100644 index 0654a39..0000000 --- a/doc/pytfa.rst +++ /dev/null @@ -1,22 +0,0 @@ -pytfa package -============= - -Subpackages ------------ - -.. toctree:: - - pytfa.analysis - pytfa.core - pytfa.io - pytfa.optim - pytfa.thermo - pytfa.utils - -Module contents ---------------- - -.. automodule:: pytfa - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/pytfa.utils.rst b/doc/pytfa.utils.rst deleted file mode 100644 index 9f02d46..0000000 --- a/doc/pytfa.utils.rst +++ /dev/null @@ -1,38 +0,0 @@ -pytfa.utils package -=================== - -Submodules ----------- - -pytfa.utils.logger module -------------------------- - -.. automodule:: pytfa.utils.logger - :members: - :undoc-members: - :show-inheritance: - -pytfa.utils.numerics module ---------------------------- - -.. automodule:: pytfa.utils.numerics - :members: - :undoc-members: - :show-inheritance: - -pytfa.utils.str module ----------------------- - -.. automodule:: pytfa.utils.str - :members: - :undoc-members: - :show-inheritance: - - -Module contents ---------------- - -.. automodule:: pytfa.utils - :members: - :undoc-members: - :show-inheritance: diff --git a/doc/requirements.txt b/doc/requirements.txt new file mode 100644 index 0000000..70a2ca7 --- /dev/null +++ b/doc/requirements.txt @@ -0,0 +1,2 @@ +Sphinx +sphinx-autoapi \ No newline at end of file diff --git a/docker/run b/docker/run old mode 100644 new mode 100755 diff --git a/models/.gitattributes b/models/.gitattributes new file mode 100644 index 0000000..4d5d840 --- /dev/null +++ b/models/.gitattributes @@ -0,0 +1 @@ +*.mat filter=lfs diff=lfs merge=lfs -text diff --git a/models/glycolysis.mat b/models/glycolysis.mat index 1d28901..cbe3681 100644 Binary files a/models/glycolysis.mat and b/models/glycolysis.mat differ diff --git a/models/small_ecoli.mat b/models/small_ecoli.mat index 7c1b3e1..b9b59b9 100644 Binary files a/models/small_ecoli.mat and b/models/small_ecoli.mat differ diff --git a/pytfa/analysis/sampling.py b/pytfa/analysis/sampling.py index 1f942e5..b57e3e3 100644 --- a/pytfa/analysis/sampling.py +++ b/pytfa/analysis/sampling.py @@ -13,7 +13,7 @@ import numpy as np from sympy.core.singleton import S from time import time -from cobra.flux_analysis.sampling import OptGPSampler, ACHRSampler, HRSampler,\ +from cobra.sampling import OptGPSampler, ACHRSampler, HRSampler,\ shared_np_array from optlang.interface import OPTIMAL @@ -29,6 +29,8 @@ def __init__(self, model, thinning, nproj=None, seed=None): # This currently has to be done to reset the solver basis which is # required to get deterministic warmup point generation # (in turn required for a working `seed` arg) + HRSampler.__init__(self, model, thinning, seed=seed) + if model.solver.is_integer: raise TypeError("sampling does not work with integer problems :(") self.model = model.copy() @@ -200,4 +202,4 @@ def sample(model, n, method="optgp", thinning=100, processes=1, seed=None): else: raise ValueError("method must be 'optgp' or 'achr'!") - return sampler.sample(n, fluxes = False) \ No newline at end of file + return sampler.sample(n, fluxes = False) diff --git a/pytfa/analysis/variability.py b/pytfa/analysis/variability.py index de9730e..79d4c92 100644 --- a/pytfa/analysis/variability.py +++ b/pytfa/analysis/variability.py @@ -19,6 +19,7 @@ from cobra.core import Reaction from optlang.exceptions import SolverError from optlang.interface import INFEASIBLE +from tqdm import tqdm from ..optim import DeltaG from ..optim.constraints import ForbiddenProfile @@ -155,7 +156,7 @@ def variability_analysis(tmodel, kind='reactions', proc_num = BEST_THREAD_RATIO) results = {'min':{}, 'max':{}} for sense in ['min','max']: - for k,var in these_vars.items(): + for k,var in tqdm(these_vars.items(), desc=sense+'imizing'): tmodel.logger.debug(sense + '-' + k) results[sense][k] = _variability_analysis_element(tmodel,var,sense) diff --git a/pytfa/core/model.py b/pytfa/core/model.py index 7093bbf..ea156ac 100644 --- a/pytfa/core/model.py +++ b/pytfa/core/model.py @@ -19,8 +19,9 @@ from cobra.core.solution import Solution from ..utils.str import camel2underscores -from ..optim.variables import GenericVariable -from ..optim.utils import get_primal +from ..optim.variables import GenericVariable, ReactionVariable, MetaboliteVariable +from ..optim.constraints import ReactionConstraint, MetaboliteConstraint +from ..optim.utils import get_primal, get_all_subclasses import time @@ -42,7 +43,7 @@ def timed(self, *args, **kw): message = '%r (%r, %r) %2.2f sec' % (method.__name__, args, kw, te-ts) try: - self.logger.info(message) + self.logger.debug(message) except AttributeError: print(message) return result @@ -148,6 +149,49 @@ def add_constraint(self, kind, hook, expr, queue=False,**kwargs): return cons + def remove_reactions(self, reactions, remove_orphans=False): + # Remove the constraints and variables associated to these reactions + all_cons_subclasses = get_all_subclasses(ReactionConstraint) + all_var_subclasses = get_all_subclasses(ReactionVariable) + + self._remove_associated_consvar(all_cons_subclasses, all_var_subclasses, + reactions) + + Model.remove_reactions(self,reactions,remove_orphans) + + def remove_metabolites(self, metabolite_list, destructive=False): + # Remove the constraints and variables associated to these reactions + all_cons_subclasses = get_all_subclasses(MetaboliteConstraint) + all_var_subclasses = get_all_subclasses(MetaboliteVariable) + + self._remove_associated_consvar(all_cons_subclasses, all_var_subclasses, + metabolite_list) + + Model.remove_metabolites(self, metabolite_list, destructive) + + def _remove_associated_consvar(self, all_cons_subclasses, all_var_subclasses, + collection): + + if not hasattr(collection, '__iter__'): + collection = [collection] + + strfy = lambda x:x if isinstance(x, str) else x.id + + for cons_type in all_cons_subclasses: + for element in collection: + try: + cons = self._cons_kinds[cons_type.__name__].get_by_id(strfy(element)) + self.remove_constraint(cons) + except KeyError: + pass + for var_type in all_var_subclasses: + for element in collection: + try: + var = self._var_kinds[var_type.__name__].get_by_id(strfy(element)) + self.remove_variable(var) + except KeyError: + pass + def remove_variable(self, var): """ Removes a variable diff --git a/pytfa/optim/utils.py b/pytfa/optim/utils.py index 6a259fd..f69badb 100644 --- a/pytfa/optim/utils.py +++ b/pytfa/optim/utils.py @@ -25,6 +25,22 @@ SYMPY_ADD_CHUNKSIZE = 100 INTEGER_VARIABLE_TYPES = ('binary','integer') +def get_all_subclasses(cls): + """ + Given a variable or constraint class, get all the subclassses + that inherit from it + + :param cls: + :return: + """ + all_subclasses = [] + + for subclass in cls.__subclasses__(): + all_subclasses.append(subclass) + all_subclasses.extend(get_all_subclasses(subclass)) + + return all_subclasses + def chunk_sum(variables): """ This functions handles the sum of many sympy variables by chunks, which diff --git a/pytfa/redgem/__init__.py b/pytfa/redgem/__init__.py new file mode 100644 index 0000000..40a96af --- /dev/null +++ b/pytfa/redgem/__init__.py @@ -0,0 +1 @@ +# -*- coding: utf-8 -*- diff --git a/pytfa/redgem/debugging.py b/pytfa/redgem/debugging.py new file mode 100644 index 0000000..3c2a671 --- /dev/null +++ b/pytfa/redgem/debugging.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python + +# -*- coding: utf-8 -*- + +""" +.. module:: redgem + :platform: Unix, Windows + :synopsis: RedGEM Algorithm + +.. moduleauthor:: pyTFA team + +Debugging +""" + +from cobra import Reaction +from pandas import Series + +def make_sink(met, ub=100, lb=0): + rid = 'sink_' + met.id + try: + # if the sink already exists + # (happens when debugging several times a model) + new = met.model.reactions.get_by_id(rid) + except KeyError: + new = Reaction(id = rid) + new.add_metabolites({met:-1}) + + new.lower_bound = lb + new.upper_bound = ub + + return new + +def add_BBB_sinks(model,biomass_rxn_id, ub=100, lb=0): + + bio_rxn = model.reactions.get_by_id(biomass_rxn_id) + + all_BBBs = bio_rxn.reactants + all_sinks = list() + + for the_BBB in all_BBBs: + new = make_sink(the_BBB, ub=ub, lb=lb) + all_sinks.append(new) + + model.add_reactions([x for x in all_sinks if not x.id in model.reactions]) + return [model.reactions.get_by_id(x.id) for x in all_sinks] + +def check_BBB_production(model, biomass_rxn_id, verbose = False): + + all_sinks = add_BBB_sinks(model, biomass_rxn_id, lb = 0) + + prod = dict() + + for the_sink in all_sinks: + with model: + model.objective = the_sink.id + prod[the_sink.id] = model.slim_optimize() + + ret = Series(prod) + if verbose: + print(ret) + return ret + +def min_BBB_uptake(model,biomass_rxn_id, min_growth_value, verbose=False): + + + with model: + all_sinks = add_BBB_sinks(model, biomass_rxn_id, ub = 0, lb = -100) + # Uptake is negative + # Min absolute uptake = Max uptake + bio_rxn = model.reactions.get_by_id(biomass_rxn_id) + bio_rxn.lower_bound = min_growth_value + model.objective = sum( - 1* s.reverse_variable for s in all_sinks) + model.objective_direction = 'max' + model.optimize() + ret = Series({r:r.flux for r in all_sinks}) + + if verbose: + print(ret) + return ret \ No newline at end of file diff --git a/pytfa/redgem/lumpgem.py b/pytfa/redgem/lumpgem.py new file mode 100644 index 0000000..2a9fdb3 --- /dev/null +++ b/pytfa/redgem/lumpgem.py @@ -0,0 +1,524 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +from cobra import Reaction + +from ..optim.utils import symbol_sum +from ..thermo.utils import is_exchange, check_transport_reaction +from .utils import trim_epsilon_mets + +from ..optim.variables import ReactionVariable, BinaryVariable, get_binary_type +from ..optim.constraints import ReactionConstraint, ForbiddenProfile + +from numpy import sum, round + +from optlang.interface import INFEASIBLE, TIME_LIMIT, OPTIMAL + +from tqdm import tqdm + +from collections import defaultdict, namedtuple + +CPLEX = 'optlang-cplex' +GUROBI = 'optlang-gurobi' +GLPK = 'optlang-glpk' + + +# Transforms (OnePerBBB --> oneperbbb), (one_per_bbb --> oneperbbb), etc ... +disambiguate = lambda s:s.lower().replace('_','') + +Lump = namedtuple('Lump', ['id_', 'metabolites', 'subnetwork', 'gene_reaction_rule']) + +class InfeasibleExcept(Exception): + def __init__(self, status, feasibility): + self.status = status + self.feasibility = feasibility + + +class TimeoutExcept(Exception): + def __init__(self, time_limit): + self.time_limit = time_limit + + +class FluxKO(ReactionVariable, BinaryVariable): + prefix = 'KO_' + + def __init__(self, reaction, **kwargs): + ReactionVariable.__init__(self, reaction, + type=get_binary_type(), + **kwargs) + +# Define a new constraint type: +class UseOrKOInt(ReactionConstraint): + prefix = 'UKI_' +# Define a new constraint type: +class UseOrKOFlux(ReactionConstraint): + prefix = 'UKF_' + + +class LumpGEM: + """ + A class encapsulating the LumpGEM algorithm + """ + def __init__(self, tfa_model, additional_core_reactions, params): + """ + :param tfa_model: The GEM (associated with the thermodynamics constraints) that lumpGEM must work on + :type tfa_model: pytfa model + + :param biomass_rxns: list of biomass reactions + :type biomass_rxns: [GEM.biomass_rxn.id] + + :param core_subsystems: list of Core subsystems names + :type core_subsystems: [string] + + :param growth_rate: theoretical maximum specific growth rate in 1/hr units + :type growth_rate: float + + :param timeout_limit: the maximum amount of time allowed to compute each optimization. Default is 3600s (1 hour) + :type timeout_limit: float (seconds) + """ + + self._tfa_model = tfa_model + + self._param_dict = params + self.init_params() + + # Set containing every BBB reaction + self._rBBB = list() + # Set containing every exchange reaction + self._exchanges = list() + # Set containing every transport reaction + self._transports = list() + # Set containing every core reaction + self._rcore = list() + # Set containing every non-core reaction + self._rncore = list() + + # For each reaction + for rxn in self._tfa_model.reactions: + # If it's a BBB reaction + if rxn.id in self.biomass_rxns: + self._rBBB.append(rxn) + # If it is an exchange reaction + elif is_exchange(rxn): + self._exchanges.append(rxn) + # If it is a transport reaction + elif check_transport_reaction(rxn): + self._transports.append(rxn) + # If it's a core reaction + elif rxn.subsystem in self.core_subsystems: + self._rcore.append(rxn) + # If it is part of the intrasubsystem expansion + elif rxn.id in additional_core_reactions: + self._rcore.append(rxn) + # If it's neither BBB nor core, then it's non-core + else: + self._rncore.append(rxn) + + # Growth rate + self._growth_rate = self.growth_rate + + # TODO : solver choice + # TODO default : solver du modele + #self._solver = 'optlang-cplex' + + self._tfa_model.solver.configuration.timeout = self.timeout_limit + print("Timeout limit is {}s".format(self.timeout_limit)) + + # lumpgem binary variables to deactivate non-core reactions. The reaction is deactivated when the value of + # the variable is 1 + self._activation_vars = {rxn: self._tfa_model.add_variable(kind=FluxKO, + hook=rxn, + lb=0, + ub=1, + queue=False) + for rxn in self._rncore} + + self._generate_usage_constraints() + self._generate_objective() + self._sinks = self._prepare_sinks() + + def init_params(self): + self.core_subsystems = self._param_dict["core_subsystems"] + self.extracellular_system = self._param_dict["extracellular_system"] + self.biomass_rxns = self._param_dict["biomass_rxns"] + + self.growth_rate = self._param_dict["growth_rate"] + + self.small_metabolites = self._param_dict["small_metabolites"] + self.cofactor_pairs = self._param_dict["cofactor_pairs"] + # Flatten cofactor_pairs list + self.cofactors = [cofactor for pair in self.cofactor_pairs for cofactor in pair] + self.inorganics = self._param_dict["inorganics"] + + self.timeout_limit = self._param_dict["timeout"] + + self.constraint_method = self._param_dict["constraint_method"] + + def _generate_usage_constraints(self): + """ + Generate carbon intake related constraints for each non-core reaction + For each reaction rxn : rxn.forward_variable + rxn.reverse_variable + activation_var * C_uptake < C_uptake + """ + flux_methods = ['flux', 'fluxes', 'both'] + int_methods = ['int', 'integer', 'both'] + + if self.constraint_method.lower() not in flux_methods + int_methods: + raise ArgumentError('{} is not a correct constraint method. ' + 'Choose among [Flux, Integer, Both]. ' + 'If you do not know what to choose, go for Flux.' + 'If it is too slow, go for integer.' + 'If you get strange lumps, go for both' + .format(self.constraint_method)) + + for rxn in self._rncore: + activation_var = self._activation_vars[rxn] + if self.constraint_method.lower() in flux_methods: + bigM = 100 + reac_var = rxn.forward_variable + rxn.reverse_variable + activation_var * bigM + # adding the constraint to the model + self._tfa_model.add_constraint(kind=UseOrKOFlux, + hook=rxn, + expr=reac_var, + ub=bigM, + lb=0, + queue=True) + if self.constraint_method.lower() in int_methods: + fu = self._tfa_model.forward_use_variable .get_by_id(rxn.id) + bu = self._tfa_model.backward_use_variable.get_by_id(rxn.id) + reac_var = fu + bu + activation_var + # adding the constraint to the model + self._tfa_model.add_constraint(kind=UseOrKOInt, + hook=rxn, + expr=reac_var, + ub=1, + lb=0, + queue=True) + + + # push constraints in one bulk (faster) + self._tfa_model._push_queue() + # refresh constraint fields + self._tfa_model.repair() + + def get_cofactor_adjusted_stoich(self,rxn): + stoich_dict = {x.id:v for x,v in rxn.metabolites.items()} + + for a,b in self.cofactor_pairs: + try: + na = stoich_dict[a] # looks like -54 atp_c + nb = stoich_dict[b] # looks like +53 adp_c + + n = na+nb # looks like -1 + + if n == 0: + self._tfa_model.logger.warn( + 'Cofactor pair {}/{} is equimolar in reaction {}' + .format(a,b,rxn.id)) + elif n > 0: + n = -n + self._tfa_model.logger.warn( + 'Cofactor pair {}/{} looks inverted in reaction {}' + .format(a,b,rxn.id)) + + stoich_dict[a] = n # looks like 1 + stoich_dict[b] = -n # looks like -1 + except KeyError: + pass + return stoich_dict + + + def _prepare_sinks(self): + """ + For each BBB (reactant of the biomass reactions), generate a sink, i.e an unbalanced reaction BBB -> + of which purpose is to enable the BBB to be output of the GEM + :return: the dict {BBB: sink} containing every BBB (keys) and their associated sinks + """ + all_sinks = {} + print("Preparing sinks...") + + for bio_rxn in self._rBBB: + stoich_dict = self.get_cofactor_adjusted_stoich(bio_rxn) + for met in bio_rxn.reactants: + stoech_coeff = stoich_dict[met.id] + # stoech_coeff < 0 indicates that the metabolite is a reactant + if (stoech_coeff < 0) and (met not in all_sinks.keys()): + sink = Reaction("Sink_" + bio_rxn.id + "_" + met.id) + sink.name = "Sink_" + bio_rxn.name + "_" + met.name + # Subsystem specific to BBB sinks + sink.subsystem = "Demand" + + # A sink is simply a reaction which consumes the BBB + sink.add_metabolites({met: -1}) + + # The sinks will be activated later (cf compute_lumps), one at a time + # sink.knock_out() + + # The stoechiometric coefficients will be used to define the lower bound of the sink, + # thus it must be stored + all_sinks[met] = (sink.id, -stoech_coeff) + self._tfa_model.add_reactions([sink]) + + # reactant already seen + elif stoech_coeff < 0: + # The BBB has already been associated to a sink, so we simply increase the bound of the sink + all_sinks[met][1] -= stoech_coeff + + # Must be called before changing the reaction.thermo['computed'] values + self._tfa_model.prepare() + for ncrxn in self._rncore: + ncrxn.thermo['computed'] = False + + return all_sinks + + def _generate_objective(self): + """ + Generate and add the maximization objective : set as many activation variables as possible to 1 + When an activation variable is set to 1, the corresponding non-core reaction is deactivated + """ + # Sum of binary variables to be maximized + objective_sum = symbol_sum(list(self._activation_vars.values())) + # Set the sum as the objective function + self._tfa_model.objective = self._tfa_model.problem.Objective(objective_sum, direction='max') + + def compute_lumps(self, force_solve=False, method='OnePerBBB'): + """ + For each BBB (reactant of the biomass reaction), add the corresponding sink to the model, then optimize and + lump the result into one lumped reaction + :param force_solve: Indicates whether the computations must continue when one lumping yields a status "infeasible" + :return: The dict {BBB: lump} containing every lumped reactions, associated to their BBBs + """ + + # Must be called before optimization + self._tfa_model.convert() + # self._tfa_model.objective_direction = 'min' + + epsilon = self._tfa_model.solver.configuration.tolerances.feasibility + + the_method = disambiguate(method) + print('Lumping method detected: {}'.format(the_method)) + + # dict: {metabolite: lumped_reaction} + lumps = {} + + self._tfa_model.objective_direction = 'max' + + sink_iter = tqdm(self._sinks.items(), desc = 'met') + + for met_BBB, (sink_id, stoech_coeff) in sink_iter: + + # Cute stuff + sink_iter.set_description('met={}'.format(met_BBB.id[:10])) + sink_iter.refresh() + + sink = self._tfa_model.reactions.get_by_id(sink_id) + # Activate reaction by setting its lower bound + prev_lb = sink.lower_bound + min_prod = self._growth_rate * stoech_coeff + sink.lower_bound = min_prod - epsilon + + if the_method == 'oneperbbb': + this_lump = self._lump_one_per_bbb(met_BBB, sink, force_solve) + lumped_reactions = [this_lump] if this_lump is not None else list() + elif the_method.startswith('min+'): + try: + p = int(the_method.replace('min+','')) + except ValueError: + raise ValueError('Min+p method must have p as an integer') + lumped_reactions = self._lump_min_plus_p(met_BBB, sink, p, force_solve) + elif the_method.startswith('min'): + lumped_reactions = self._lump_min_plus_p(met_BBB, sink, 0, force_solve) + else: + raise ValueError('Lumping method not recognized: {}. ' + 'Valid methods are ' + 'OnePerBBB, Min, Min+p, p natural integer' + .format(the_method)) + + + if not lumped_reactions: + continue + + lumps[met_BBB] = lumped_reactions + + # Deactivating reaction by setting both bounds to 0 + sink.lower_bound = prev_lb + # sink.knock_out() + + self.lumps = lumps + return lumps + + def _lump_one_per_bbb(self, met_BBB, sink, force_solve): + """ + + :param met_BBB: + :param sink: + :param force_solve: + :return: + """ + + n_da = self._tfa_model.slim_optimize() + + try: + # Timeout reached + if self._tfa_model.solver.status == TIME_LIMIT: + raise TimeoutExcept(self._tfa_model.solver.configuration.timeout) + # Not optimal status -> infeasible + elif self._tfa_model.solver.status != OPTIMAL: + raise InfeasibleExcept(self._tfa_model.solver.status, + self._tfa_model.solver.configuration.tolerances.feasibility) + + except (TimeoutExcept, InfeasibleExcept) as err: + # If the user want to continue anyway, suits him + if force_solve: + # Raise a warning + return None + else: + raise err + + # print('Produced {}'.format(sink.flux), + # 'with {0:.0f} reactions deactivated'.format(n_da)) + + lumped_reaction = self._build_lump(met_BBB, sink) + + return lumped_reaction + + + def _lump_min_plus_p(self, met_BBB, sink, p, force_solve): + """ + + :param met_BBB: + :param sink: + :param force_solve: + :return: + """ + + epsilon = self._tfa_model.solver.configuration.tolerances.integrality + + try: + max_lumps =self._param_dict['max_lumps_per_BBB'] + except KeyError: + # TODO: Put a warning + max_lumps=10 + + lumps = list() + + with self._tfa_model as model: + activation_vars = model.get_variables_of_type(FluxKO) + + # Solve a first time, obtain minimal subnet + model.slim_optimize() + max_deactivated_rxns = model.objective.value + + # Add constraint forbidding subnets bigger than p + expr = symbol_sum(activation_vars) + + # The lower bound is the max number of deactivated, minus p + # Which allows activating the minimal number of reactions, plus p + lb = max_deactivated_rxns - p + model.add_constraint(kind=ForbiddenProfile, + hook = model, + id_ = 'MAX_DEACT_{}'.format(met_BBB.id), + expr = expr, + lb = lb, + ub = max_deactivated_rxns, + ) + + n_deactivated_reactions = max_deactivated_rxns + + # While loop, break on infeasibility + while len(lumps) 0: + self._intermediate_metabolites_id[subsystem_i][subsystem_j][d].add(path[i]) + + def find_min_distance_between_subsystems(self): + """ + Find minimal distance between each subsystems in both directions + + :return: Dict with distances + """ + for i in self._core_subsystems: + for j in self._core_subsystems: + for k in range(self._d+1): + # If there path of length d + if self._intermediate_paths[i][j][k]: + self._min_distance_sub_to_sub[i][j] = k + break + # If min distance os not found, then + if self._min_distance_sub_to_sub[i][j] == -1: + pass + return self._min_distance_sub_to_sub + + def breadth_search_extracellular_system_paths(self, subsystem, n): + """ + Breadth first search from each metabolite in the extracellular system with special stop + conditions during exploration for paths of length n. + + This function explores the graph through allowed paths only : this path can't go through + the extracellular system or the subsystem but must start in the extracellular system and + end in the subsystem. The length of each path found is n. + + :param subsystem: Destination subsystem + :param n: Path length desired + :return: None + """ + for metabolite_id in self._extracellular_system: + # Find metabolites at a distance n from metabolite_id + if metabolite_id not in self._graph: + continue + ancestors = {} + frontier = {metabolite_id} + explored = {metabolite_id} + for i in range(n): + new_nodes = set() + for current_node in frontier: + for new_node in set(self._graph.adj[current_node]): + if self.is_node_allowed_extracellular(new_node, i, explored, subsystem, n): + new_nodes.add(new_node) + # new_node can already be in ancestors if there are 2 paths of same + # length to it + if new_node in ancestors: + ancestors[new_node].append(current_node) + else: + ancestors[new_node] = [current_node] + explored = explored.union(new_nodes) + frontier = new_nodes + # Handle n = 0 case, since it didn't go through the loop + if n == 0 and metabolite_id not in self._subsystem_metabolites_id[subsystem]: + frontier = {} + # Retrieve and save metabolites, reactions and paths + for node in frontier: + paths = self.retrieve_all_paths(node, metabolite_id, ancestors) + self._intermediate_extracellular_paths[subsystem][n] = \ + self._intermediate_extracellular_paths[subsystem][n].union(set(paths)) + self.retrieve_intermediate_extracellular_metabolites_and_reactions(paths, subsystem, + n) + + def is_node_allowed_extracellular(self, node, i, explored, subsystem, n): + """ + Checks whether or not a metabolite is allowed for the current path. + + The new node is added if it is not already explored, if it is not in the extracellular + system, and if it is not in the destination subsystem except if it is the last round + of exploration + + :param node: Metabolite id + :param i: Current step + :param explored: Explored node for this path + :param subsystem: Destination subsystem + :param n: Path length desired + :return: Boolean answering the question + """ + if node in explored: + return False + if node in self._extracellular_system: + return False + if i < n-1 and node in self._subsystem_metabolites_id[subsystem]: + return False + if i == n-1 and node not in self._subsystem_metabolites_id[subsystem]: + return False + return True + + def retrieve_intermediate_extracellular_metabolites_and_reactions(self, paths, subsystem, n): + """ + Retrieves and stores intermediate metabolites and reactions for the extracellular system + + This function adds all reactions contained in these paths, and all metabolites between + + :param paths: List of paths + :param subsystem: Destination subsystem + :param n: Path length + :return: None + """ + for path in paths: + for i in range(len(path) - 1): + reaction = self._graph[path[i]][path[i + 1]]['rxn_id'] + self._intermediate_extracellular_reactions_id[subsystem][n].add(reaction) + if i > 0: + self._intermediate_extracellular_metabolites_id[subsystem][n].add(path[i]) + + def run_between_all_subsystems(self): + """ + Retrieve subsystem and intermediate reactions and metabolites. + + :return: None + """ + for subsystem in self._core_subsystems: + self.extract_subsystem_reactions(subsystem) + self.extract_subsystem_metabolites(subsystem) + + for subsystem_i in self._core_subsystems: + for subsystem_j in self._core_subsystems: + for k in range(self._d+1): + self.breadth_search_subsystems_paths_length_d(subsystem_i, subsystem_j, k) + + def run_extracellular_system(self): + """ + Retrieve intermediate reactions and metabolites for the extracellular system + + :return: None + """ + for subsystem in self._core_subsystems: + for k in range(self._n + 1): + self.breadth_search_extracellular_system_paths(subsystem, k) + + def extract_sub_network(self): + """ + Extracts the reduced gem. + + :return: None + """ + def extract_id(x): + return x.id + to_remove_metabolites = set(map(extract_id, self._redgem.metabolites)) + to_remove_reactions = set(map(extract_id, self._redgem.reactions)) + + # Keep subsystems reactions and metabolites + for name in self._core_subsystems: + to_remove_reactions = to_remove_reactions - self._subsystem_reactions_id[name] + to_remove_metabolites = to_remove_metabolites - self._subsystem_metabolites_id[name] + + # Keep intermediate reactions and metabolites + for i in self._core_subsystems: + for j in self._core_subsystems: + for k in range(self._d+1): + to_remove_reactions = to_remove_reactions \ + - self._intermediate_reactions_id[i][j][k] + to_remove_metabolites = to_remove_metabolites \ + - self._intermediate_metabolites_id[i][j][k] + + # Keep extracellular metabolites + to_remove_metabolites = to_remove_metabolites - set(self._extracellular_system) + + # Keep intermediate extracellular reactions and metabolites + for i in self._core_subsystems: + for k in range(self._d+1): + to_remove_reactions = to_remove_reactions \ + - self._intermediate_extracellular_reactions_id[i][k] + to_remove_metabolites = to_remove_metabolites \ + - self._intermediate_extracellular_metabolites_id[i][k] + + self._redgem.remove_reactions(to_remove_reactions, True) + + def run(self): + """ + Runs RedGEM. + + :return: None + """ + self.create_new_stoichiometric_matrix() + self.run_between_all_subsystems() + self.run_extracellular_system() + self.extract_sub_network() + + return self._redgem diff --git a/pytfa/redgem/redgem.py b/pytfa/redgem/redgem.py new file mode 100644 index 0000000..3a8b9d4 --- /dev/null +++ b/pytfa/redgem/redgem.py @@ -0,0 +1,191 @@ +#!/usr/bin/env python + +# -*- coding: utf-8 -*- + +""" +.. module:: redgem + :platform: Unix, Windows + :synopsis: RedGEM Algorithm + +.. moduleauthor:: pyTFA team + +Model class +""" + +from pytfa.redgem.network_expansion import NetworkExpansion +from pytfa.redgem.lumpgem import LumpGEM +from cobra import Reaction +from .utils import remove_blocked_reactions, set_medium +import yaml + +class RedGEM(): + def __init__(self, gem, parameters_path, inplace=False): + + self.read_parameters(parameters_path) + + # If inplace is True, no deepcopy is performed : the modifications are applied directly onto the gem + prepared_gem = set_medium(gem, self.params['medium'], inplace) + + self._gem = prepared_gem + # This one is used to perform the lumping + self._source_gem = prepared_gem.copy() + + self.logger = self._gem.logger + + self.fill_default_params() + + self.set_solver() + + def read_parameters(self, parameters_path): + with open(parameters_path, 'r') as stream: + try: + self.params = yaml.safe_load(stream) + print("Opened parameters file") + except yaml.YAMLError as exc: + print(exc) + + def fill_default_params(self): + # If auto is activated, automatically extracts inorganics from the gem + if "inorganic" not in self.params or self.params["inorganics"] == "auto": + self.logger.info("Automatically computing inorganics to use") + self.params["inorganics"] = self._extract_inorganics() + if "growth_rate" not in self.params or self.params["growth_rate"] == "auto": + self.logger.info("Setting minimal growth rate to 95% of the TFA solution") + obj_val = self._source_gem.slim_optimize() + self.logger.info("Setting minimal growth rate to {}".format(obj_val)) + self.params["growth_rate"] = 0.95*obj_val + if "force_solve" not in self.params: + self.params["force_solve"] = False + if "timeout" not in self.params: + self.logger.info("Using default timeout : 3600s") + self.params["timeout"] = 3600 + if "feasibility" not in self.params: + self.logger.info("Using default solver feasibility : 1e-9") + self.params["feasibility"] = 1e-9 + else: + # numbers like 1e-9 are detected as strings by yaml module + # to enable their use, we cast them into floats + try: + self.params["feasibility"] = float(self.params["feasibility"]) + except ValueError as v: + self.logger.error(v) + + def set_solver(self): + if "solver" not in self.params or self.params["solver"].lower() == "auto": + return None + elif 'gurobi' in self.params["solver"].lower(): + solver = 'gurobi' + elif 'cplex' in self.params["solver"].lower(): + solver = 'cplex' + elif 'glpk' in self.params["solver"].lower(): + solver = 'glpk' + else: + solver = self.params["solver"] + + self._gem.solver = solver + self._source_gem.solver = solver + + + def run(self): + # Extracting parameters + core_subsystems = self.params["core_subsystems"] + extracellular_system = self.params["extracellular_system"] + biomass_rxn_ids = self.params["biomass_rxns"] + + biomass_rxns = [self._gem.reactions.get_by_id(x) for x in biomass_rxn_ids] + main_bio_rxn = biomass_rxns[0] + + growth_rate = self.params["growth_rate"] + + small_metabolites = self.params["small_metabolites"] + cofactor_pairs = self.params["cofactor_pairs"] + # Flatten cofactor_pairs list + cofactors = [cofactor for pair in cofactor_pairs for cofactor in pair] + inorganics = self.params["inorganics"] + + d = self.params["d"] + n = self.params["n"] + lump_method = self.params["lump_method"] + + force_solve = self.params["force_solve"] + timeout = self.params["timeout"] + self._gem.solver.configuration.tolerances.feasibility = self.params["feasibility"] + self._gem.solver.configuration.tolerances.integrality = self.params["feasibility"] + self._source_gem.solver.configuration.tolerances.feasibility = self.params["feasibility"] + self._source_gem.solver.configuration.tolerances.integrality = self.params["feasibility"] + + self.logger.info("Computing network expansion...") + expander = NetworkExpansion(self._gem, core_subsystems, extracellular_system, + cofactors, small_metabolites, inorganics, + d, n) + reduced_gem = expander.run() + self.logger.info("Done.") + + # Add the expansion to core reactions + core_reactions = reduced_gem.reactions + + self.logger.info("Computing lumps...") + lumper = LumpGEM(self._source_gem, core_reactions, self.params) + lumps = lumper.compute_lumps(force_solve, method = lump_method) + self.logger.info("Done.") + + self.logger.info("Create final network...") + to_add = [x for x in biomass_rxns + +lumper._exchanges + +lumper._transports + +lumper._rcore + if not x.id in reduced_gem.reactions] + reduced_gem.add_reactions(to_add) + + for rxns in lumps.values(): + the_lumps = [add_lump(reduced_gem,rxn,id_suffix='_{}'.format(e)) + for e,rxn in enumerate(rxns)] + # reduced_gem.add_reactions(rxns) + self.logger.info("Done.") + + reduced_gem.objective = main_bio_rxn + reduced_gem.reactions.get_by_id(main_bio_rxn.id).lower_bound = growth_rate + + if self.params['remove_blocked_reactions']: + self.logger.info('Detecting blocked reactions') + # Remove blocked reactions + nrxn_1 = len(reduced_gem.reactions) + self.removed_reactions = remove_blocked_reactions(reduced_gem) + nrxn_2 = len(reduced_gem.reactions) + self.logger.info('Removed {} blocked reaction with ' + 'FVA post-processing'.format(nrxn_1-nrxn_2)) + + if main_bio_rxn.id not in reduced_gem.reactions: + raise RuntimeError('Main Biomass reaction appears blocked') + + + # For debugging purposes + self.lumper = lumper + main_bio_rxn.lower_bound = 0 + + return reduced_gem + + def _extract_inorganics(self): + """ + Extract inorganics from self._gem based on their formula + + :return: list of inorganics metabolites + """ + + inorganics = [] + for met in self._gem.metabolites: + if not met.elements == {}: # Edge case + # met is inorganic if it has 0 carbon in its formula + if (not 'C' in met.elements) or met.elements['C'] <= 0: + inorganics.append(met.id) + + return inorganics + +def add_lump(model, lump_object, id_suffix=''): + new = Reaction(id = lump_object.id_+id_suffix) + model.add_reaction(new) + new.add_metabolites(lump_object.metabolites) + new.gene_reaction_rule = lump_object.gene_reaction_rule + new.subnetwork = lump_object.subnetwork + + return new diff --git a/pytfa/redgem/utils.py b/pytfa/redgem/utils.py new file mode 100644 index 0000000..df4f7b1 --- /dev/null +++ b/pytfa/redgem/utils.py @@ -0,0 +1,42 @@ +import numpy as np +# from cobra.flux_analysis import flux_variability_analysis +from pytfa.analysis.variability import variability_analysis + +def remove_blocked_reactions(model): + epsilon = model.solver.configuration.tolerances.feasibility + # fva = flux_variability_analysis(model) + fva = variability_analysis(model, kind='reactions') + # Blocked reactions have max and min at 0 + df = fva[ (fva.max(axis=1).abs()<1*epsilon) + & (fva.min(axis=1).abs()<1*epsilon)] + rid_to_rm = df.index + + model.remove_reactions(rid_to_rm) + + return df + + +def trim_epsilon_mets(met_dict, epsilon): + + n = int(-1*np.log10(epsilon)) + round_dict = {x:np.round(v,n) for x,v in met_dict.items()} + met_dict.update(round_dict) + + rm_list = [x for x,v in met_dict.items() if abs(v) <= epsilon] + [met_dict.pop(x) for x in rm_list] + + return met_dict + +def set_medium(model, medium_dict, inplace): + if inplace: + new = model + else: + new = model.copy() + + if medium_dict is None or not medium_dict: + return new + + for rxn_id, lb in medium_dict.items(): + new.reactions.get_by_id(rxn_id).lower_bound = lb + + return new \ No newline at end of file diff --git a/pytfa/thermo/tmodel.py b/pytfa/thermo/tmodel.py index a568bc2..d80c96f 100644 --- a/pytfa/thermo/tmodel.py +++ b/pytfa/thermo/tmodel.py @@ -136,11 +136,11 @@ def _prepare_metabolite(self,met): # Which index of the reaction DB do you correspond to ? if not 'seed_id' in met.annotation: # raise Exception("seed_id missing for " + met.name) - self.logger.warning("Metabolite {} ({}) has no seed_id".\ + self.logger.debug("Metabolite {} ({}) has no seed_id".\ format(met.id, met.name)) metData = None elif not met.annotation['seed_id'] in self.compounds_data: - self.logger.warning("Metabolite {} ({}) not present in thermoDB" + self.logger.debug("Metabolite {} ({}) not present in thermoDB" .format(met.annotation['seed_id'], met.name)) metData = None else: @@ -201,7 +201,7 @@ def _prepare_reaction(self,reaction): or len(reaction.metabolites) >= 100 or balanceResult in ['missing atoms', 'drain flux']): - self.logger.info('{} : thermo constraint NOT created'.format(reaction.id)) + self.logger.debug('{} : thermo constraint NOT created'.format(reaction.id)) reaction.thermo['computed'] = False reaction.thermo['deltaGR'] = BIGM_DG reaction.thermo['deltaGRerr'] = BIGM_DG @@ -346,7 +346,7 @@ def _convert_metabolite(self, met, add_potentials, verbose): metDeltaGF) else: - self.logger.info('NOT generating thermo variables for {}'.format(met.id)) + self.logger.debug('NOT generating thermo variables for {}'.format(met.id)) if LC != None: # Register the variable to find it more easily diff --git a/pytfa/thermo/utils.py b/pytfa/thermo/utils.py index 6441cbf..3b93430 100644 --- a/pytfa/thermo/utils.py +++ b/pytfa/thermo/utils.py @@ -185,4 +185,8 @@ def check_transport_reaction(reaction): def is_same_stoichiometry(this_reaction, that_reaction): this_met_dict = {k.id:v for k,v in this_reaction.metabolites.items()} that_met_dict = {k.id:v for k,v in that_reaction.metabolites.items()} - return this_met_dict == that_met_dict \ No newline at end of file + return this_met_dict == that_met_dict + + +def is_exchange(rxn): + return len(rxn.metabolites) == 1 \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 172a54a..ad2fdc4 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,8 @@ cobra>=0.11.3 bokeh>=0.12.1 +networkx optlang sympy pytest scipy -tqdm +tqdm \ No newline at end of file diff --git a/setup.py b/setup.py index 58d57f0..e73435c 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,7 @@ # return reqs -version_tag = '0.8.1' +version_tag = '0.9.0-b0' setup(name='pytfa', version=version_tag, @@ -30,6 +30,7 @@ download_url='https://github.com/EPFL-LCSB/pytfa/archive/'+version_tag+'.tar.gz', install_requires=['cobra>0.11', 'bokeh>=0.12.1', + 'networkx', 'optlang', 'pytest', 'scipy', diff --git a/tests/redgem_params.yml b/tests/redgem_params.yml new file mode 100644 index 0000000..3f1dcb3 --- /dev/null +++ b/tests/redgem_params.yml @@ -0,0 +1,64 @@ +--- +core_subsystems: + # - "Citric Acid Cycle" + # - "Pentose Phoshate Pathway" + - "Glycolysis/Gluconeogenesis" + # - "Pyruvate Metabolism" + # - "Glyoxylate Metabolism" + - 'Oxidative Phosphorylation' + +extracellular_system: + - DM_ac_e + - DM_etoh_e + - DM_glyc_e + - DM_lac-D_e + - DM_akg_e + - DM_for_e + - DM_pyr_e + - DM_fum_e + - DM_co2_e + - DM_mal-L_e + +biomass_rxns: + - Ec_biomass_iJO1366_WT_53p95M + +# Define minimal growth to reach. If auto, sets to 95% of max TFA growth +growth_rate: auto + +# Define the medium. Any other reaction will be kept as provided by the model +medium: + DM_cbl1_e: -0.01 + DM_glc_e: -8.16 + DM_glycogenn1_c: -10.0 + +small_metabolites: + - h_c + - h2o_c + - co2_c + - o2_c + - pi_c + - ppi_c + - nh4_c + - h2o2_c + +inorganics: auto + +cofactor_pairs: + - ['atp_c', 'adp_c'] + - ['nad_c', 'nadh_c'] + - ['nadp_c', 'nadph_c'] + +# Reduction parameters +lump_method: Min+1 # OnePerBBB, Min, Min+p +max_lumps_per_BBB: 2 # Only used if method = Min or Min+p +remove_blocked_reactions: no # Remove reactions that cannot carry flux + # at the end of the reduction + +d: 2 +n: 2 + +solver: auto +force_solve: False +timeout: 300 +feasibility: 1e-9 +constraint_method: integer # flux, integer, or both diff --git a/tests/redgem_params_textbook.yaml b/tests/redgem_params_textbook.yaml new file mode 100644 index 0000000..1211d86 --- /dev/null +++ b/tests/redgem_params_textbook.yaml @@ -0,0 +1,56 @@ +--- +core_subsystems: + - "Glycolysis" + - "Cofactor Regeneration" + +extracellular_system: + - EX_co2_e + - EX_glc__D_e + - EX_h_e + - EX_h2o_e + - EX_nh4_e + - EX_o2_e + - EX_pi_e + +biomass_rxns: + - Biomass_Ecoli_core + +# Define minimal growth to reach. If auto, sets to 95% of max TFA growth +growth_rate: auto + +# Define the medium. Any other reaction will be kept as provided by the model +medium: + EX_glc__D_e: -10 + +small_metabolites: + - h_c + - h2o_c + - co2_c + - o2_c + - pi_c + - ppi_c + - nh4_c + - h2o2_c + +inorganics: auto + +cofactor_pairs: + - ['atp_c', 'adp_c'] + - ['nad_c', 'nadh_c'] + - ['nadp_c', 'nadph_c'] + - ['accoa_c', 'coa_c'] + +# Reduction parameters +lump_method: Min+1 # OnePerBBB, Min, Min+p +max_lumps_per_BBB: 2 # Only used if method = Min or Min+p +remove_blocked_reactions: no # Remove reactions that cannot carry flux + # at the end of the reduction + +d: 2 +n: 2 + +solver: auto +force_solve: False +timeout: 300 +feasibility: 1e-9 +constraint_method: flux # flux, integer, or both diff --git a/tests/test_redgem.py b/tests/test_redgem.py new file mode 100644 index 0000000..00285fb --- /dev/null +++ b/tests/test_redgem.py @@ -0,0 +1,79 @@ +from pytfa.redgem.redgem import RedGEM + +from pytfa.io import import_matlab_model +from pytfa.io.base import load_thermoDB +from pytfa.thermo.tmodel import ThermoModel +from pytfa.io import read_compartment_data, apply_compartment_data, read_lexicon, annotate_from_lexicon +from settings import this_directory +from os.path import join + +import pytest + +# Check if we are running on Travis CI, to make the run lighter +import os +is_travis = 'TRAVIS' in os.environ \ + or 'CI' in os.environ \ + or os.environ.get('USER')=='travis' +DEBUG = False + +if is_travis or DEBUG: + # # Remove 80% of the mets of the biomass reaction so that less lumps need to be computed: + # print('Travis env detected. Trimming the biomass reaction') + # bio_rxn = model.reactions.get_by_id('Ec_biomass_iJO1366_WT_53p95M') + # bio_rxn.add_metabolites({k:-v for e,(k,v) in + # enumerate(bio_rxn.metabolites.items()) + # if e != 1}) + # # if e%10 != 0}) + from cobra.test import create_test_model + model = create_test_model('textbook') + gly_rxns = ['ENO','FBA','FBP','GAPD','PDH','PFK','PGI','PGK','PGM','PPS', + 'PYK','TPI'] + cofactor_regen_rxns = ['NADTRHD','NADH16','NADTRHD','ATPM','ATPS4r'] + for x in gly_rxns: + model.reactions.get_by_id(x).subsystem = 'Glycolysis' + for x in cofactor_regen_rxns: + model.reactions.get_by_id(x).subsystem = 'Cofactor Regeneration' + + path_to_params = join(this_directory, '..', 'tests/redgem_params_textbook.yaml') + + +else: + path_to_model = join(this_directory, '..', 'models/small_ecoli.mat') + model = import_matlab_model(path_to_model) + path_to_params = join(this_directory, '..', 'tests/redgem_params.yml') + +thermoDB = join(this_directory, '..', 'data/thermo_data.thermodb') +path_to_lexicon = join(this_directory, '..', 'models/iJO1366/lexicon.csv') +path_to_compartment_data = join(this_directory, '..', 'models/iJO1366/compartment_data.json') + + +# Scaling to avoid numerical errors with bad lumps +for rxn in model.reactions: + if rxn.id.startswith('LMPD_'): + rxn.add_metabolites({x:v*(0.1 - 1) for x,v in rxn.metabolites.items()}) + +thermo_data = load_thermoDB(thermoDB) +lexicon = read_lexicon(path_to_lexicon) +compartment_data = read_compartment_data(path_to_compartment_data) + +tfa_model = ThermoModel(thermo_data, model) +annotate_from_lexicon(tfa_model, lexicon) +apply_compartment_data(tfa_model, compartment_data) + +tfa_model.name = 'Lumped Model' +tfa_model.prepare() +tfa_model.convert() + +# tfa_model.solver.configuration.verbosity = True +tfa_model.logger.setLevel = 30 + +def test_redgem(): + redgem = RedGEM(tfa_model, path_to_params, False) + rgem = redgem.run() + obj_val = rgem.slim_optimize() + # assert(obj_val > 0) + return rgem, redgem + + +if __name__ == '__main__': + rgem, redgem = test_redgem() \ No newline at end of file diff --git a/tutorials/tuto_redgem_params.yaml b/tutorials/tuto_redgem_params.yaml new file mode 100644 index 0000000..ba6b67d --- /dev/null +++ b/tutorials/tuto_redgem_params.yaml @@ -0,0 +1,71 @@ +--- +core_subsystems: + - "Citric Acid Cycle" + - "Pentose Phoshate Pathway" + - "Glycolysis/Gluconeogenesis" + - "Pyruvate Metabolism" + - "Glyoxylate Metabolism" + - 'Oxidative Phosphorylation' + +extracellular_system: + - DM_ac_e + - DM_etoh_e + - DM_glyc_e + - DM_lac-D_e + - DM_akg_e + - DM_for_e + - DM_pyr_e + - DM_fum_e + - DM_co2_e + - DM_mal-L_e + +biomass_rxns: + - Ec_biomass_iJO1366_WT_53p95M + +# Define minimal growth to reach. If auto, sets to 95% of max TFA growth +growth_rate: auto + +# Define the medium. Any other reaction will be kept as provided by the model +medium: +# DM_cbl1_e: -0.01 +# DM_glc_e: -8.16 +# DM_glycogenn1_c: -10.0 +# DM_bglycogenn1_c: -10.0 + +small_metabolites: + - h_c + - h2o_c + - co2_c + - o2_c + - pi_c + - ppi_c + - nh4_c + - h2o2_c + +inorganics: auto + +cofactor_pairs: + - ['atp_c', 'adp_c'] + - ['nad_c', 'nadh_c'] + - ['nadp_c', 'nadph_c'] + +# Reduction parameters +lump_method: Min+1 # OnePerBBB, Min, Min+p +max_lumps_per_BBB: 2 # Only used if method = Min or Min+p +remove_blocked_reactions: yes # Remove reactions that cannot carry flux + # at the end of the reduction + +d: 2 +n: 2 + +solver: auto +force_solve: False +timeout: 300 +feasibility: 1e-9 + +# Optimisation shenanigans +# Integer is faster with a good parallel solver, but might yield some numerical errors +# Flux is slower but more reliable +# Both is good if you can afford it +constraint_method: both # flux, integer, or both. + diff --git a/tutorials/tutorial_redgem.py b/tutorials/tutorial_redgem.py new file mode 100644 index 0000000..67b33e4 --- /dev/null +++ b/tutorials/tutorial_redgem.py @@ -0,0 +1,76 @@ +from pytfa.redgem.redgem import RedGEM + +from pytfa.io import import_matlab_model +from pytfa.io.base import load_thermoDB +from pytfa.thermo.tmodel import ThermoModel +from pytfa.io import read_compartment_data, apply_compartment_data, \ + read_lexicon, annotate_from_lexicon +from cobra.io import load_json_model + +from os.path import join + +# Paths +path_to_model = join('..','models','iJO1366.json') +thermoDB = join('..','data','thermo_data.thermodb') +path_to_lexicon = join('..','models','iJO1366','lexicon.csv') +path_to_compartment_data = join('..','models','iJO1366','compartment_data.json') + + +# FBA +model = load_json_model(path_to_model) + +fba_solution = model.optimize() +fba_value = fba_solution.objective_value + +# Thermo prep + +thermo_data = load_thermoDB(thermoDB) +lexicon = read_lexicon(path_to_lexicon) +compartment_data = read_compartment_data(path_to_compartment_data) + +# Initialize the cobra_model +tfa_model = ThermoModel(thermo_data, model) +tfa_model.name = 'Tutorial' + +# Annotate the cobra_model +annotate_from_lexicon(tfa_model, lexicon) +apply_compartment_data(tfa_model, compartment_data) + +tfa_model.prepare() +tfa_model.convert() + +# tfa_model.solver.configuration.verbosity = True +tfa_model.logger.setLevel = 30 + +tfa_solution = tfa_model.optimize() +tfa_value = tfa_solution.objective_value + +# It might happen that the model is infeasible. In this case, we can relax +# thermodynamics constraints: + +if tfa_value < 0.1: + from pytfa.optim.relaxation import relax_dgo + + biomass_rxn = 'Ec_biomass_iJO1366_WT_53p95M' + tfa_model.reactions.get_by_id(biomass_rxn).lower_bound = 0.9 * fba_value + relaxed_model, slack_model, relax_table = relax_dgo(tfa_model, in_place=True) + + original_model, tfa_model = tfa_model, relaxed_model + + print('Relaxation: ') + print(relax_table) + + tfa_solution = tfa_model.optimize() + tfa_value = tfa_solution.objective_value + +path_to_params = 'tuto_redgem_params.yaml' +redgem = RedGEM(tfa_model, path_to_params, False) +rgem = redgem.run() + +redgem_solution = rgem.optimize() +redgem_value = redgem_solution.objective_value + +# Report +print('FBA Solution found : {0:.5g}'.format(fba_value)) +print('TFA Solution found : {0:.5g}'.format(tfa_value)) +print('redGEM Solution found : {0:.5g}'.format(redgem_value))