From e7a96d9cdfa0aece53ec8befca064b6b3c0ab0eb Mon Sep 17 00:00:00 2001 From: Haris Zafeiropoulos Date: Thu, 25 Jan 2024 17:34:48 +0100 Subject: [PATCH] support of medium modification and cobra model loading (#76) * remove commented line * address issues #74 and #75 * import Dict type * test for boundary change and decreasing of fba when glc goes down * edit test for media indices * edit test for modify_medium() --- dingo/MetabolicNetwork.py | 97 +++++++++++++++++++++++++++++++++-- dingo/loading_models.py | 104 +++++++++++--------------------------- pyproject.toml | 1 - tests/fba.py | 43 +++++++++++++++- 4 files changed, 164 insertions(+), 81 deletions(-) diff --git a/dingo/MetabolicNetwork.py b/dingo/MetabolicNetwork.py index f7d1da69..eec1a8b5 100644 --- a/dingo/MetabolicNetwork.py +++ b/dingo/MetabolicNetwork.py @@ -8,7 +8,9 @@ import numpy as np import sys -from dingo.loading_models import read_json_file, read_mat_file, read_sbml_file +from typing import Dict +import cobra +from dingo.loading_models import read_json_file, read_mat_file, read_sbml_file, parse_cobra_model from dingo.fva import slow_fva from dingo.fba import slow_fba @@ -34,7 +36,7 @@ def __init__(self, tuple_args): except ImportError as e: self._parameters["fast_computations"] = False - if len(tuple_args) != 7: + if len(tuple_args) != 10: raise Exception( "An unknown input format given to initialize a metabolic network object." ) @@ -46,6 +48,9 @@ def __init__(self, tuple_args): self._reactions = tuple_args[4] self._biomass_index = tuple_args[5] self._biomass_function = tuple_args[6] + self._medium = tuple_args[7] + self._medium_indices = tuple_args[8] + self._exchanges = tuple_args[9] try: if self._biomass_index is not None and ( @@ -83,13 +88,22 @@ def from_mat(cls, arg): @classmethod def from_sbml(cls, arg): - if (not isinstance(arg, str)) or (arg[-3:] != "xml"): + if (not isinstance(arg, str)) and ((arg[-3:] == "xml") or (arg[-4] == "sbml")): raise Exception( "An unknown input format given to initialize a metabolic network object." ) return cls(read_sbml_file(arg)) + @classmethod + def from_cobra_model(cls, arg): + if (not isinstance(arg, cobra.core.model.Model)): + raise Exception( + "An unknown input format given to initialize a metabolic network object." + ) + + return cls(parse_cobra_model(arg)) + def fva(self): """A member function to apply the FVA method on the metabolic network.""" @@ -146,6 +160,14 @@ def biomass_index(self): def biomass_function(self): return self._biomass_function + @property + def medium(self): + return self._medium + + @property + def exchanges(self): + return self._exchanges + @property def parameters(self): return self._parameters @@ -160,6 +182,9 @@ def get_as_tuple(self): self._reactions, self._biomass_index, self._biomass_function, + self._medium, + self._inter_medium, + self._exchanges ) def num_of_reactions(self): @@ -196,6 +221,72 @@ def biomass_index(self, value): def biomass_function(self, value): self._biomass_function = value + + @medium.setter + def medium(self, medium: Dict[str, float]) -> None: + """Set the constraints on the model exchanges. + + `model.medium` returns a dictionary of the bounds for each of the + boundary reactions, in the form of `{rxn_id: rxn_bound}`, where `rxn_bound` + specifies the absolute value of the bound in direction of metabolite + creation (i.e., lower_bound for `met <--`, upper_bound for `met -->`) + + Parameters + ---------- + medium: dict + The medium to initialize. medium should be a dictionary defining + `{rxn_id: bound}` pairs. + """ + + def set_active_bound(reaction: str, reac_index: int, bound: float) -> None: + """Set active bound. + + Parameters + ---------- + reaction: cobra.Reaction + Reaction to set + bound: float + Value to set bound to. The bound is reversed and set as lower bound + if reaction has reactants (metabolites that are consumed). If reaction + has reactants, it seems the upper bound won't be set. + """ + if any(x < 0 for x in list(self._S[:, reac_index])): + self._lb[reac_index] = -bound + elif any(x > 0 for x in list(self._S[:, reac_index])): + self._ub[reac_index] = bound + + # Set the given media bounds + media_rxns = [] + exchange_rxns = frozenset(self.exchanges) + for rxn_id, rxn_bound in medium.items(): + if rxn_id not in exchange_rxns: + logger.warning( + f"{rxn_id} does not seem to be an an exchange reaction. " + f"Applying bounds anyway." + ) + media_rxns.append(rxn_id) + + reac_index = self._reactions.index(rxn_id) + + set_active_bound(rxn_id, reac_index, rxn_bound) + + frozen_media_rxns = frozenset(media_rxns) + + # Turn off reactions not present in media + for rxn_id in exchange_rxns - frozen_media_rxns: + """ + is_export for us, needs to check on the S + order reactions to their lb and ub + """ + # is_export = rxn.reactants and not rxn.products + reac_index = self._reactions.index(rxn_id) + products = np.any(self._S[:,reac_index] > 0) + reactants_exist = np.any(self._S[:,reac_index] < 0) + is_export = True if not products and reactants_exist else False + set_active_bound( + rxn_id, reac_index, min(0.0, -self._lb[reac_index] if is_export else self._ub[reac_index]) + ) + def set_fast_mode(self): try: diff --git a/dingo/loading_models.py b/dingo/loading_models.py index 4acf9ca1..84efceb5 100644 --- a/dingo/loading_models.py +++ b/dingo/loading_models.py @@ -23,80 +23,15 @@ def read_json_file(input_file): input_file -- a json file that contains the information about a mettabolic network, for example see http://bigg.ucsd.edu/models """ - with open(input_file, "r") as f: - - distros_dict = json.load(f) - - reactions_list = distros_dict["reactions"] - - metabolites = [] - reactions = [] - - for reaction in reactions_list: - - metabolites_dic = reaction["metabolites"] - reaction_name = reaction["id"] - reactions.append(reaction_name) - - for metabolite in metabolites_dic.keys(): - if metabolite not in metabolites: - metabolites.append(metabolite) - - list_of_reaction_lists = [] - vector_of_ubs = [] - vector_of_lbs = [] - - for reaction in reactions_list: - - ub = float(reaction["upper_bound"]) - vector_of_ubs.append(ub) - lb = float(reaction["lower_bound"]) - vector_of_lbs.append(lb) - - metabolites_dic = reaction["metabolites"] - reaction_vector = [] - - for term in range(len(metabolites)): - - metabolite = metabolites[term] - - if metabolite in metabolites_dic.keys(): - - reaction_vector.append(metabolites_dic[metabolite]) - else: - reaction_vector.append(0) - - list_of_reaction_lists.append(reaction_vector) - - # Build function's output; - - # lower and upper flux bounds - lb = np.asarray(vector_of_lbs) - ub = np.asarray(vector_of_ubs) - - lb = np.asarray(lb, dtype="float") - lb = np.ascontiguousarray(lb, dtype="float") - - ub = np.asarray(ub, dtype="float") - ub = np.ascontiguousarray(ub, dtype="float") - - # The stoichiometric martrix S - S = np.asarray(list_of_reaction_lists) - S = np.transpose(S) - - S = np.asarray(S, dtype="float") - S = np.ascontiguousarray(S, dtype="float") + try: + cobra.io.load_matlab_model( input_file ) + except: + cobra_config = cobra.Configuration() + cobra_config.solver = 'glpk' - # Get biomass function if there - biomass_function = np.zeros(S.shape[1]) - biomass_index = None - for i in reactions: - j = i.casefold() - if "biom" in j: - biomass_index = reactions.index(i) - biomass_function[biomass_index] = 1 + model = cobra.io.load_json_model( input_file ) - return lb, ub, S, metabolites, reactions, biomass_index, biomass_function + return (parse_cobra_model( model )) def read_mat_file(input_file): """A Python function based on the to read a .mat file and returns, @@ -111,7 +46,7 @@ def read_mat_file(input_file): input_file -- a mat file that contains a MATLAB structure with the information about a mettabolic network, for example see http://bigg.ucsd.edu/models """ try: - cobra.io.load_matlab_model( input_file ) + cobra.io.load_matlab_model( input_file ) except: cobra_config = cobra.Configuration() cobra_config.solver = 'glpk' @@ -137,7 +72,7 @@ def read_sbml_file(input_file): https://github.com/VirtualMetabolicHuman/AGORA/blob/master/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/Abiotrophia_defectiva_ATCC_49176.xml """ try: - cobra.io.read_sbml_model( input_file ) + cobra.io.read_sbml_model( input_file ) except: cobra_config = cobra.Configuration() cobra_config.solver = 'glpk' @@ -188,4 +123,23 @@ def parse_cobra_model(cobra_model): ub = np.asarray(ub, dtype="float") ub = np.ascontiguousarray(ub, dtype="float") - return lb, ub, S, metabolites, reactions, biomass_index, biomass_function + medium = cobra_model.medium + inter_medium = {} + + for index, reaction in enumerate(cobra_model.reactions): + for ex_reaction in medium.keys(): + if ex_reaction == reaction.id: + inter_medium[ex_reaction] = index + + exchanges_cobra_reactions = cobra_model.exchanges + exchanges = [] + for reac in exchanges_cobra_reactions: + exchanges.append(reac.id) + + + return lb, ub, S, metabolites, reactions, biomass_index, biomass_function, medium, inter_medium, exchanges + + + + + diff --git a/pyproject.toml b/pyproject.toml index 35b38b59..2e83c04e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,6 @@ argparse = "^1.4.0" matplotlib = "^3.4.1" cobra = "^0.26.0" plotly = "^5.11.0" -#kaleido = "0.2.1" kaleido = "0.2.1" [tool.poetry.dev-dependencies] diff --git a/tests/fba.py b/tests/fba.py index 9f993d3c..94dc5b6f 100644 --- a/tests/fba.py +++ b/tests/fba.py @@ -10,7 +10,7 @@ from dingo import MetabolicNetwork class TestFba(unittest.TestCase): - + def test_fba_json(self): input_file_json = os.getcwd() + "/ext_data/e_coli_core.json" @@ -21,7 +21,7 @@ def test_fba_json(self): self.assertTrue(abs(res[1] - 0.8739215067486387) < 1e-03) def test_fba_mat(self): - + input_file_mat = os.getcwd() + "/ext_data/e_coli_core.mat" model = MetabolicNetwork.from_mat(input_file_mat) model.set_slow_mode() @@ -40,6 +40,45 @@ def test_fba_sbml(self): self.assertTrue(abs(res[1] - 0.8739215067486387) < 1e-03) + def test_modify_medium(self): + + input_file_sbml = os.getcwd() + "/ext_data/e_coli_core.xml" + model = MetabolicNetwork.from_sbml(input_file_sbml) + model.set_slow_mode() + + initial_medium = model.medium + initial_fba = model.fba()[-1] + + e_coli_core_medium_compound_indices = { + "EX_co2_e" : 46, + "EX_glc__D_e" : 51, + "EX_h_e" : 54, + "EX_h2o_e" : 55, + "EX_nh4_e" : 58, + "EX_o2_e" : 59, + "EX_pi_e" : 60 + } + + glc_index = model.reactions.index("EX_glc__D_e") + o2_index = model.reactions.index("EX_o2_e") + + new_media = initial_medium.copy() + new_media["EX_glc__D_e"] = 1.5 + new_media["EX_o2_e"] = -0.5 + + model.medium = new_media + + updated_media = model.medium + updated_medium_indices = {} + for reac in updated_media: + updated_medium_indices[reac] = model.reactions.index(reac) + + self.assertTrue(updated_medium_indices == e_coli_core_medium_compound_indices) + + self.assertTrue(model.lb[glc_index] == -1.5 and model.lb[o2_index] == 0.5) + + self.assertTrue(initial_fba - model.fba()[-1] > 0) + if __name__ == "__main__": unittest.main()