corrections to the PreProcess PR

dingo/preprocess.py

@@ -7,10 +7,27 @@
 
 class PreProcess:
 
-    def __init__(self, model, open_exchanges=False):
+    def __init__(self, model, tol=1e-6, open_exchanges=False):
 
+        """
+        model parameter gets a cobra model as input
+        
+        tol parameter gets a cutoff value used to classify 
+        zero-flux and mle reactions and compare FBA solutions
+        before and after reactions removal
+        
+        open_exchanges parameter is used in the function that identifies blocked reactions
+        It controls whether or not to open all exchange reactions to very high flux ranges
+        """
+
         self._model = model
+        self._tol = tol
+
+        if self._tol > 1e-6:
+            print("Tolerance value set to",self._tol,"while default value is 1e-6. A looser check will be performed")
+
         self._open_exchanges = open_exchanges
+
         self._objective = self._objective_function()
         self._initial_reactions = self._initial()
         self._reaction_bounds_dict = self._reaction_bounds_dictionary()
@@ -82,7 +99,7 @@ def _zero_flux(self):
         when running a FVA analysis with the fraction of optimum set to 90%
         """
 
-        tol = 1e-6
+        tol = self._tol
 
         fva = cobra.flux_analysis.flux_variability_analysis(self._model, fraction_of_optimum=0.9)
         zero_flux = fva.loc[ (abs(fva['minimum']) < tol ) & (abs(fva['maximum']) < tol)]
@@ -112,7 +129,7 @@ def _metabolically_less_efficient(self):
         set to 0.95, the reactions that have no flux are the metabolically less efficient.
         """
 
-        tol = 1e-6
+        tol = self._tol
 
         fba_solution = self._model.optimize()
 
@@ -177,7 +194,7 @@ def reduce(self, extend=False):
         remained_reactions = list((Counter(self._initial_reactions)-Counter(self._removed_reactions)).elements())
         remained_reactions = list((Counter(remained_reactions)-Counter(self._essential_reactions)).elements())
 
-        tol = 1e-6
+        tol = self._tol
 
         if extend != False and extend != True:
             raise Exception("Wrong Input to extend parameter")
@@ -222,22 +239,15 @@ def reduce(self, extend=False):
             fba_solution_final = self._model.optimize().objective_value
 
 
-            # if FBA solution after removal is infesible or altered
+            # if FBA solution after removal is infeasible or altered
             # restore the initial reactions bounds
-            if (fba_solution_final == None):
+            if (fba_solution_final == None) | (abs(fba_solution_final - fba_solution_initial) > tol):
                 for reaction in additional_removed_reactions_list:
                     self._model.reactions.get_by_id(reaction).bounds = self._reaction_bounds_dict[reaction]
                     self._removed_reactions.remove(reaction)
                 print(len(self._removed_reactions), "of the", len(self._initial_reactions), \
                 "reactions were removed from the model with extend set to", extend)
 
-            elif(abs(fba_solution_final - fba_solution_initial) > tol):
-                for reaction in additional_removed_reactions_list:
-                    self._model.reactions.get_by_id(reaction).bounds = self._reaction_bounds_dict[reaction]
-                    self._removed_reactions.remove(reaction)
-                print(len(self._removed_reactions), "of the", len(self._initial_reactions), \
-                "reactions were removed from the model with extend set to", extend) 
-
             else:
                 print(len(self._removed_reactions), "of the", len(self._initial_reactions), \
                 "reactions were removed from the model with extend set to", extend)

tests/preprocess.py

@@ -1,64 +1,66 @@
 
 from cobra.io import load_json_model
+from dingo import MetabolicNetwork
 from dingo.preprocess import PreProcess
 import unittest
 import numpy as np
-import cobra
 
 class TestPreprocess(unittest.TestCase):
 
     def test_preprocess(self):
 
         # load cobra model
         cobra_model = load_json_model("ext_data/e_coli_core.json")
-
-        # find reaction ids from the loaded model
-        initial_reactions_ids = []    
-        for reaction in cobra_model.reactions:
-            reaction_id = reaction.id
-            initial_reactions_ids.append(reaction_id)
+        
+        # convert cobra to dingo model
+        initial_dingo_model = MetabolicNetwork.from_cobra_model(cobra_model)
+
+        # perform an FBA to find the initial FBA solution
+        initial_fba_solution = initial_dingo_model.fba()[1]
 
 
         # call the reduce function from the PreProcess class 
         # with extend=False to remove reactions from the model        
-        obj = PreProcess(cobra_model, open_exchanges=False)  
-        removed_reactions, dingo_model = obj.reduce(extend=False)      
+        obj = PreProcess(cobra_model, tol=1e-5, open_exchanges=False)  
+        removed_reactions, final_dingo_model = obj.reduce(extend=False)      
 
         # calculate the count of removed reactions with extend set to False        
         removed_reactions_count = len(removed_reactions)
         self.assertTrue( 46 - removed_reactions_count == 0 )
 
         # calculate the count of reactions with bounds equal to 0 
         # with extend set to False from the dingo model
-        dingo_removed_reactions = np.sum((dingo_model.lb == 0) & (dingo_model.ub == 0))
+        dingo_removed_reactions = np.sum((final_dingo_model.lb == 0) & (final_dingo_model.ub == 0))
         self.assertTrue( 46 - dingo_removed_reactions == 0 )
 
-        # perform an FBA to check the result after reactions removal
-        res = dingo_model.fba()
-        self.assertTrue(abs(res[1] - 0.8739215067486387) < 1e-03)       
-
+        # perform an FBA to check the solution after reactions removal
+        final_fba_solution = final_dingo_model.fba()[1]
+        self.assertTrue(abs(final_fba_solution - initial_fba_solution) < 1e-03)       
 
 
         # load models in cobra and dingo format again to restore bounds
-        cobra_model = load_json_model("ext_data/e_coli_core.json")        
-
+        cobra_model = load_json_model("ext_data/e_coli_core.json")
+
+        # convert cobra to dingo model
+        initial_dingo_model = MetabolicNetwork.from_cobra_model(cobra_model)
+
         # call the reduce function from the PreProcess class 
         # with extend=True to remove additional reactions from the model        
-        obj = PreProcess(cobra_model, open_exchanges=False)        
-        removed_reactions, dingo_model = obj.reduce(extend=True)        
+        obj = PreProcess(cobra_model, tol=1e-6, open_exchanges=False)        
+        removed_reactions, final_dingo_model = obj.reduce(extend=True)        
 
         # calculate the count of removed reactions with extend set to True        
         removed_reactions_count = len(removed_reactions)
         self.assertTrue( 47 - removed_reactions_count == 0 )
 
         # calculate the count of reactions with bounds equal to 0 
         # with extend set to True from the dingo model
-        dingo_removed_reactions = np.sum((dingo_model.lb == 0) & (dingo_model.ub == 0))
+        dingo_removed_reactions = np.sum((final_dingo_model.lb == 0) & (final_dingo_model.ub == 0))
         self.assertTrue( 47 - dingo_removed_reactions == 0 )
 
         # perform an FBA to check the result after reactions removal
-        res = dingo_model.fba()
-        self.assertTrue(abs(res[1] - 0.8739215067486387) < 1e-03)   
+        final_fba_solution = final_dingo_model.fba()[1]
+        self.assertTrue(abs(final_fba_solution - initial_fba_solution) < 1e-03)   
 
 
 if __name__ == "__main__":