Merge pull request #78 from Fxe/dev

Fxe · web-flow · commit 8a83e2ab6bb9 · 2022-07-20T06:21:40.000-05:00
Unit tests ATP method and Gap fill
diff --git a/.gitignore b/.gitignore
@@ -1,4 +1,5 @@
 .DS_Store
+@eaDir
 
 # Byte-compiled / optimized / DLL files
 __pycache__/
diff --git a/modelseedpy/core/msatpcorrection.py b/modelseedpy/core/msatpcorrection.py
@@ -19,7 +19,7 @@ class MSATPCorrection:
 
     DEBUG = False
 
-    def __init__(self, model, core_template, atp_medias, compartment="c0",
+    def __init__(self, model, core_template, atp_medias: list, compartment="c0",
                  max_gapfilling=None, gapfilling_delta=0, atp_hydrolysis_id=None):
         """
 
@@ -45,8 +45,8 @@ def __init__(self, model, core_template, atp_medias, compartment="c0",
             self.atp_hydrolysis = output["reaction"]
         self.atp_medias = []
         for media in atp_medias:
-            if isinstance(media,MSMedia):
-                self.atp_medias.append([media,0.01])
+            if isinstance(media, MSMedia):
+                self.atp_medias.append([media, 0.01])
             else:
                 self.atp_medias.append(media)        
         self.max_gapfilling = max_gapfilling
@@ -163,25 +163,26 @@ def evaluate_growth_media(self):
         with self.model:
             self.model.objective = self.atp_hydrolysis.id
             #self.model.objective = self.model.problem.Objective(Zero,direction="max")
-            #self.atp_hydrolysis.update_variable_bounds()
+
             logger.debug(f'ATP bounds: ({self.atp_hydrolysis.lower_bound}, {self.atp_hydrolysis.upper_bound})')
             #self.model.objective.set_linear_coefficients({self.atp_hydrolysis.forward_variable:1})
             pkgmgr = MSPackageManager.get_pkg_mgr(self.model)
-            for media_tuple in self.atp_medias:
-                media = media_tuple[0]
+            for media, minimum_obj in self.atp_medias:
                 logger.debug('evaluate media %s', media)
                 pkgmgr.getpkg("KBaseMediaPkg").build_package(media)
                 logger.debug('model.medium %s', self.model.medium)
                 solution = self.model.optimize()
                 logger.debug('evaluate media %s - %f (%s)', media.id, solution.objective_value, solution.status)
                 self.media_gapfill_stats[media] = None
                 output[media.id] = solution.objective_value
-                if solution.objective_value < media_tuple[1] or solution.status != 'optimal':
-                    self.media_gapfill_stats[media] = self.msgapfill.run_gapfilling(media, self.atp_hydrolysis.id,media_tuple[1])
+                if solution.objective_value < minimum_obj or solution.status != 'optimal':
+                    self.media_gapfill_stats[media] = self.msgapfill.run_gapfilling(media,
+                                                                                    self.atp_hydrolysis.id,
+                                                                                    minimum_obj)
                     #IF gapfilling fails - need to activate and penalize the noncore and try again
-                elif solution.objective_value >= media_tuple[1]:
+                elif solution.objective_value >= minimum_obj:
                     self.media_gapfill_stats[media] = {'reversed': {}, 'new': {}}
-                logger.debug('gapfilling stats:',json.dumps(self.media_gapfill_stats[media],indent=2))
+                logger.debug('gapfilling stats: %s', json.dumps(self.media_gapfill_stats[media], indent=2))
 
         if MSATPCorrection.DEBUG:
             with open('debug.json', 'w') as outfile:
@@ -335,10 +336,13 @@ def run_atp_correction(self):
         self.evaluate_growth_media()
         self.determine_growth_media()
         self.apply_growth_media_gapfilling()
+        self.evaluate_growth_media()
         self.expand_model_to_genome_scale()
+        return self.build_tests()
     
     @staticmethod
-    def atp_correction(model,coretemplate,atp_medias = None,atp_objective = "bio2",max_gapfilling = None,gapfilling_delta = 0):
-        msatpobj = MSATPCorrection(model,coretemplate,atp_medias,atp_objective,max_gapfilling,gapfilling_delta)
-        msatpobj.run_atp_correction()
-        return msatpobj
+    def atp_correction(model, core_template, atp_medias=None, atp_objective="bio2",
+                       max_gapfilling=None, gapfilling_delta=0):
+        atp_correction = MSATPCorrection(model, core_template, atp_medias, atp_hydrolysis_id=atp_objective,
+                                         max_gapfilling=max_gapfilling, gapfilling_delta=gapfilling_delta)
+        return atp_correction.run_atp_correction()
diff --git a/tests/core/test_msatpcorreption.py b/tests/core/test_msatpcorreption.py
@@ -11,10 +11,16 @@ def template():
         return MSTemplateBuilder.from_dict(json.load(fh)).build()
 
 
+@pytest.fixture
+def template_genome_scale():
+    with open('./tests/test_data/template_genome_scale_bigg.json', 'r') as fh:
+        return MSTemplateBuilder.from_dict(json.load(fh)).build()
+
+
 @pytest.fixture
 def get_model():
 
-    def _method(ko=None):
+    def _method(ko=None, added_compounds=None, added_reactions=None):
         if ko is None:
             ko = []
         with open('./tests/test_data/e_coli_core.json', 'r') as fh:
@@ -36,6 +42,13 @@ def _method(ko=None):
                     r['id'] += '_' + 'c0'  # hack cause there is only combo between e0 and c0
 
             model_json['reactions'] = [x for x in model_json['reactions'] if x['id'] not in ko]
+
+            if added_compounds:
+                for o in added_compounds:
+                    model_json['metabolites'].append(o)
+            if added_reactions:
+                for o in added_reactions:
+                    model_json['reactions'].append(o)
             model = cobra.io.from_json(json.dumps(model_json))
             model.reactions.ATPM_c0.lower_bound = 0
             model.reactions.ATPM_c0.upper_bound = 1000
@@ -68,25 +81,143 @@ def media_acetate_aerobic():
     return media
 
 
+@pytest.fixture
+def media_genome_scale_glucose_aerobic():
+    media = MSMedia.from_dict({
+        'glc__D': (-10, 1000),
+        'o2': (-1000, 1000),
+        'h': (-1000, 1000),
+        'h2o': (-1000, 1000),
+        'pi': (-1000, 1000),
+        'co2': (-1000, 1000),
+        'nh4': (-1000, 1000),
+        'k': (-1000, 1000),
+    })
+    return media
+
+
 @pytest.fixture
 def media_all_aerobic(media_glucose_aerobic, media_acetate_aerobic):
     return [media_glucose_aerobic, media_acetate_aerobic]
 
 
+@pytest.fixture
+def get_model_with_infinite_atp_loop(get_model, template_genome_scale):
+
+    def _method(ko=None):
+        added_compounds = [{
+            'id': 'k_c0',
+            'name': 'K [c]',
+            'compartment': 'c0',
+            'charge': 1,
+            'formula': 'K',
+            'notes': {},
+            'annotation': {'sbo': 'SBO:0000247'}
+        }, {
+            'id': 'k_e0',
+            'name': 'K [e]',
+            'compartment': 'e0',
+            'charge': 1,
+            'formula': 'K',
+            'notes': {},
+            'annotation': {'sbo': 'SBO:0000247'}
+        }]
+        added_reactions = [{
+            'id': 'EX_k_e0',
+            'name': 'K exchange',
+            'metabolites': {'k_e0': -1.0},
+            'lower_bound': -1000,
+            'upper_bound': 1000.0,
+            'gene_reaction_rule': '',
+            'subsystem': 'Extracellular exchange',
+            'notes': {},
+            'annotation': {}
+        }]
+        model = get_model(ko, added_compounds, added_reactions)
+        model.reactions.get_by_id('BIOMASS_Ecoli_core_w_GAM_c0').add_metabolites({
+            model.metabolites.get_by_id('k_c0'): 0.01
+        })
+        model.add_reactions([template_genome_scale.reactions.Kt2r_c.to_reaction(model)])
+        model.add_reactions([template_genome_scale.reactions.Ktex_c.to_reaction(model)])
+        model.medium = {
+            'EX_co2_e0': 1000.0,
+            'EX_glc__D_e0': 10.0,
+            'EX_h_e0': 1000.0,
+            'EX_h2o_e0': 1000.0,
+            'EX_nh4_e0': 1000.0,
+            'EX_o2_e0': 1000.0,
+            'EX_pi_e0': 1000.0,
+            'EX_k_e0': 1000.0
+        }
+
+        return model
+
+    return _method
+
+
+def test_infinite_atp_model_growth_boost(get_model_with_infinite_atp_loop, template, template_genome_scale, media_glucose_aerobic):
+    model = get_model_with_infinite_atp_loop()
+    solution = model.optimize()
+    assert solution.objective_value > 1.2  # should be around 1.316
+    assert solution.status == 'optimal'
+
+
 def test_ms_atp_correction1(get_model, template, media_all_aerobic):
-    model = get_model(["GLCpts_c0", "NADH16_c0", "GLCpts_c0", "O2t_c0"])
+    model = get_model(["GLCpts_c0", "NADH16_c0", "CYTBD_c0", "O2t_c0"])
     atp_correction = MSATPCorrection(model, template, media_all_aerobic, atp_hydrolysis_id='ATPM_c0')
     atp_correction.evaluate_growth_media()
-    print('non_core_reactions', len(atp_correction.noncore_reactions),
-          'other_compartments', len(atp_correction.other_compartments),
-          'original_bounds', len(atp_correction.original_bounds))
-    for media in atp_correction.media_gapfill_stats:
-        print(media, atp_correction.media_gapfill_stats[media])
+    assert len(atp_correction.noncore_reactions) == 1  # the biomass
+    assert len(atp_correction.other_compartments) == 0  # none
+    assert len(atp_correction.original_bounds) == 70  # 70 reactions
+
+    """
+    glc/o2 {'reversed': {}, 'new': {'GLCpts_c0': '>'}}
+    ac/o2 {'reversed': {}, 'new': {'CYTBD_c0': '>', 'NADH16_c0': '>', 'O2t_c0': '>'}}
+    """
     atp_correction.determine_growth_media()
-    print(atp_correction.selected_media)
+
+    assert len(atp_correction.selected_media) == 1  # selects glucose
+
+    atp_correction.apply_growth_media_gapfilling()
+
     media_eval = atp_correction.evaluate_growth_media()
-    print(media_eval)
+
     atp_correction.expand_model_to_genome_scale()
     tests = atp_correction.build_tests()
-    print(tests)
-    assert True
+
+    assert tests
+    assert len(tests) == 1
+    assert tests[0]['threshold'] > 0
+    assert tests[0]['objective'] == 'ATPM_c0'
+
+
+def test_ms_atp_correction_and_gap_fill1(get_model_with_infinite_atp_loop, template, template_genome_scale,
+                                         media_glucose_aerobic, media_genome_scale_glucose_aerobic):
+    from modelseedpy import MSGapfill
+
+    model = get_model_with_infinite_atp_loop(['GLCpts_c0', 'GLUSy_c0', 'GLUDy_c0'])
+    model.reactions.ATPM_c0.lower_bound = 0
+    model.reactions.ATPM_c0.upper_bound = 1000
+
+    atp_correction = MSATPCorrection(model, template, [media_glucose_aerobic], atp_hydrolysis_id='ATPM_c0')
+    tests = atp_correction.run_atp_correction()
+
+    # expected tests = [{'media': MSMedia object, 'is_max_threshold': True, 'threshold': 21.0, 'objective': 'ATPM_c0'}]
+
+    assert tests
+    assert len(tests) == 1
+    assert tests[0]['threshold'] > 0
+    assert tests[0]['objective'] == 'ATPM_c0'
+
+    gap_fill = MSGapfill(model, [template_genome_scale], [], tests, {}, [])
+    result = gap_fill.run_gapfilling(media_genome_scale_glucose_aerobic, 'BIOMASS_Ecoli_core_w_GAM_c0', minimum_obj=0.1)
+
+    # either GLUSy_c0 or GLUDy_c0 should be gap filled for glutamate
+
+    assert result
+    assert len(result['new']) == 1
+    assert 'GLUSy_c0' in result['new'] or 'GLUDy_c0' in result['new']
+
+    model = gap_fill.integrate_gapfill_solution(result)
+
+    assert model
diff --git a/tests/core/test_msgapfill.py b/tests/core/test_msgapfill.py
@@ -70,7 +70,7 @@ def test_ms_gap_fill1(template, get_model, media_glucose_aerobic):
     """
     Test gap filling with glucose aerobic requesting minimal growth value
     """
-    model = get_model(["GLCpts_c0", "NADH16_c0", "GLCpts_c0", "O2t_c0"])
+    model = get_model(["CYTBD_c0", "NADH16_c0", "GLCpts_c0", "O2t_c0"])
     gap_fill = MSGapfill(model, [template])
     result = gap_fill.run_gapfilling(media_glucose_aerobic, 'BIOMASS_Ecoli_core_w_GAM_c0', minimum_obj=0.1)
     assert result
diff --git a/tests/test_data/template_genome_scale_bigg.json b/tests/test_data/template_genome_scale_bigg.json

Original file line number	Diff line number	Diff line change
`@@ -1,4 +1,5 @@`
`1`	`1`	`.DS_Store`
	`2`	`+@eaDir`
`2`	`3`
`3`	`4`	`# Byte-compiled / optimized / DLL files`
`4`	`5`	`__pycache__/`