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

@@ -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:

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
+
+
+
+
+

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]

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()