Merge pull request #5 from ModelSEED/dev

ups

Author: Fxe

modelseedpy/core/fbahelper.py

@@ -3,34 +3,31 @@
 import logging
 from chemicals import periodic_table
 import re
+from cobra.core import Gene, Metabolite, Model, Reaction   # !!! Gene, Metabolite, and Model are never used
+from cobra.util import solver as sutil  # !!! sutil is never used
 import time
-from cobra.core import Gene, Metabolite, Model, Reaction
-from cobra.util import solver as sutil
 from modelseedpy.biochem import from_local
-from scipy.odr.odrpack import Output
+from scipy.odr.odrpack import Output  # !!! Output is never used
+from chemw import ChemMW
+from warnings import warn
 #from Carbon.Aliases import false
 
 logger = logging.getLogger(__name__)
 
-elementmass = {}
-for element in periodic_table:
-    elementmass[element.symbol] = element.MW
-
-
 class FBAHelper:
 
     @staticmethod
-    def add_autodrain_reactions_to_community_model(model,auto_sink = ["cpd02701", "cpd11416", "cpd15302"]):
+    def add_autodrain_reactions_to_community_model(model,auto_sink = ["cpd02701", "cpd15302"]):
         #Adding missing drains in the base model
         drain_reactions = []
         for metabolite in model.metabolites:
             msid = FBAHelper.modelseed_id_from_cobra_metabolite(metabolite)
             if msid in auto_sink:
-                if msid != "cpd11416" or metabolite.compartment == "c0":
+                if metabolite.compartment == "c0":
                     met_id = metabolite.id
-                    if all(rxn not in model.reactions for rxn in [f"EX_{met_id}", f"DM_{met_id}", f"SK_{met_id}"]):
+                    if all([rxn not in model.reactions for rxn in [f"EX_{met_id}", f"DM_{met_id}", f"SK_{met_id}"]]):
                         drain_reaction = FBAHelper.add_drain_from_metabolite_id(model,metabolite.id,0,100,"DM_")
-                        if drain_reaction != None:
+                        if not drain_reaction:
                             logger.info("Adding "+met_id+" DM")
                             drain_reactions.append(drain_reaction)
         model.add_reactions(drain_reactions)
@@ -48,25 +45,26 @@ def add_drain_from_metabolite_id(model, cpd_id, uptake, excretion, prefix='EX_',
         """
         if cpd_id in model.metabolites:
             cobra_metabolite = model.metabolites.get_by_id(cpd_id)
-            drain_reaction = Reaction(id=f'{prefix}{cpd_id}',
-                                      name=prefix_name + cobra_metabolite.name,
-                                      lower_bound=-1*uptake, 
-                                      upper_bound=excretion)
+            drain_reaction = Reaction(
+                id=f'{prefix}{cpd_id}',
+                name=prefix_name + cobra_metabolite.name,
+                lower_bound = -uptake,
+                upper_bound = excretion)
             drain_reaction.add_metabolites({cobra_metabolite : -1})
             drain_reaction.annotation["sbo"] = 'SBO:0000627'    
             #model.add_reactions([drain_reaction])
             return drain_reaction
         return None
-    
+
     @staticmethod
-    def set_reaction_bounds_from_direction(reaction, direction, add=0):
+    def set_reaction_bounds_from_direction(reaction, direction, add=False):
         if direction == "<":
             reaction.lower_bound = -100
-            if add == 0:
+            if not add:
                 reaction.upper_bound = 0
         if direction == ">":
             reaction.upper_bound = 100
-            if add == 0:
+            if not add:
                 reaction.lower_bound = 0
         reaction.update_variable_bounds()
 
@@ -76,9 +74,7 @@ def set_objective_from_target_reaction(model,target_reaction,minimize = False):
         sense = "max"
         if minimize:
             sense = "min"
-        model.objective = model.problem.Objective(
-            1 * target_reaction.flux_expression,
-            direction=sense)
+        model.objective = model.problem.Objective(target_reaction.flux_expression, direction=sense)
         return target_reaction
 
     @staticmethod
@@ -87,8 +83,7 @@ def modelseed_id_from_cobra_metabolite(metabolite):
             m = re.search('^(cpd\d+)', metabolite.id)
             return m[1]
         #TODO: should check to see if ModelSEED ID is in the annotations for the compound
-        else:
-            return None
+        return None
 
     @staticmethod
     def modelseed_id_from_cobra_reaction(reaction):
@@ -101,18 +96,16 @@ def modelseed_id_from_cobra_reaction(reaction):
 
     @staticmethod
     def metabolite_mw(metabolite):
-        mw = 0
-        elements = metabolite.elements
-        for element in elements:
-            if element not in elementmass:
-                print("Missing mass for element "+element+" in compound "+metabolite.id+". Element will be ignored when computing MW")
-            else:
-                mw += elements[element]*elementmass[element]
-        return mw
+        try:
+            chem_mw = ChemMW()
+            chem_mw.mass(metabolite.formula)
+            return chem_mw.raw_mw
+        except:
+            warn("The compound "+metabolite.id+" possesses an unconventional formula {metabolite.formula}; hence, the MW cannot be computed.")
 
     @staticmethod    
     def elemental_mass():
-        return elementmass
+        return {element.symbol:element.MW for element in periodic_table}
 
     @staticmethod
     def get_modelseed_db_api(modelseed_path):
@@ -131,8 +124,8 @@ def is_biomass(reaction):
         return reaction.id[0:3] == "bio"
 
     @staticmethod
-    def exchange_hash(model):
-        exchange_hash = {}
+    def exchange_hash(model):  #!!! This function is pointless?
+        exchange_hash = {}  # !!! this variable is never used
         for reaction in model.reactions:
             if len(reaction.metabolites) == 1:
                 for metabolite in reaction.metabolites:
@@ -141,102 +134,88 @@ def exchange_hash(model):
 
     @staticmethod
     def find_reaction(model,stoichiometry):
-        output = FBAHelper.stoichiometry_to_string(stoichiometry)
-        atpstring = output[0]
+        reaction_strings = FBAHelper.stoichiometry_to_string(stoichiometry)
+        atpstring = reaction_strings[0]
         rxn_hash = FBAHelper.rxn_hash(model)
         if atpstring in rxn_hash:
             return rxn_hash[atpstring]
         return None
 
     @staticmethod
-    def msid_hash(model): 
+    def msid_hash(model):
         output = {}
-        for cpd in model.metabolites:
-            msid = FBAHelper.modelseed_id_from_cobra_metabolite(cpd)
+        for met in model.metabolites:
+            msid = FBAHelper.modelseed_id_from_cobra_metabolite(met)
             if msid != None:
                 if msid not in output:
                     output[msid] = []
-                output[msid].append(cpd)
+                output[msid].append(met)
         return output
 
     @staticmethod
     def rxn_hash(model): 
         output = {}
         for rxn in model.reactions:
-            strings = FBAHelper.stoichiometry_to_string(rxn.metabolites)
-            output[strings[0]] = [rxn,1]
-            output[strings[1]] = [rxn,-1]
+            reaction_strings = FBAHelper.stoichiometry_to_string(rxn.metabolites)
+            output[reaction_strings[0]] = [rxn,1]
+            output[reaction_strings[1]] = [rxn,-1]
         return output
 
     @staticmethod
     def rxn_compartment(reaction): 
         compartments = list(reaction.compartments)
         if len(compartments) == 1:
             return compartments[0]
-        cytosol = None
-        othercomp = None
+        cytosol = othercomp = None
         for comp in compartments:
-            if comp[0:1] != "e":
-                if comp[0:1] == "c":
-                    cytosol = comp
-                else:
-                    othercomp = comp
-        if othercomp is not None:
-            return othercomp
-        return cytosol
+            if comp[0:1] == "c":
+                cytosol = comp
+            elif comp[0:1] != "e":
+                othercomp = comp
+        return othercomp or cytosol
 
     @staticmethod
     def stoichiometry_to_string(stoichiometry):
-        reactants = []
-        products = []
+        reactants, products = [], []
         for met in stoichiometry:
-            coef = stoichiometry[met]
+            stoich = stoichiometry[met]
             if not isinstance(met, str):
-                if FBAHelper.modelseed_id_from_cobra_metabolite(met) == "cpd00067":
-                    met = None
-                else:
-                    met = met.id
-            if met != None:
-                if coef < 0:
+                met = None if FBAHelper.modelseed_id_from_cobra_metabolite(met) == "cpd00067" else met.id
+            if met:
+                if stoich < 0:
                     reactants.append(met)
                 else:
                     products.append(met)
-        reactants.sort()
-        products.sort()
-        return ["+".join(reactants)+"="+"+".join(products),"+".join(products)+"="+"+".join(reactants)]
+        return ["+".join(sorted(reactants))+"="+"+".join(sorted(products)),"+".join(sorted(products))+"="+"+".join(sorted(reactants))]
 
     @staticmethod
     def add_atp_hydrolysis(model,compartment):
-        #Searching for ATP hydrolysis compounds
-        coefs = {"cpd00002":[-1,compartment],"cpd00001":[-1,compartment],"cpd00008":[1,compartment],"cpd00009":[1,compartment],"cpd00067":[1,compartment]}
-        msids = ["cpd00002","cpd00001","cpd00008","cpd00009","cpd00067"]
+        # Searching for ATP hydrolysis compounds
+        coefs = {"cpd00002": [-1,compartment], "cpd00001": [-1,compartment], "cpd00008": [1,compartment],
+            "cpd00009": [1,compartment], "cpd00067": [1,compartment]}
         stoichiometry = {}
         id_hash = FBAHelper.msid_hash(model)
-        for msid in msids:
+        for msid, content in coefs.items():
             if msid not in id_hash:
                 logger.warning("Compound "+msid+" not found in model!")
                 return None
             else:
                 for cpd in id_hash[msid]:
-                    if cpd.compartment == coefs[msid][1]:
-                        stoichiometry[cpd] = coefs[msid][0]
+                    if cpd.compartment == content[1]:
+                        stoichiometry[cpd] = content[0]
         output = FBAHelper.find_reaction(model,stoichiometry)
-        if output != None and output[1] == ">":
+        if output and output[1] == 1:  # !!! the second element of the output is 1/0 and not a direction string
             return {"reaction":output[0],"direction":">","new":False}
-        cobra_reaction = Reaction("rxn00062_"+compartment, 
-                                  name="ATP hydrolysis", 
-                                  lower_bound=0, 
-                                  upper_bound=1000)
-        cobra_reaction.annotation["sbo"] = "SBO:0000176" #biochemical reaction
-        cobra_reaction.annotation["seed.reaction"] = "rxn00062"
+        cobra_reaction = Reaction("rxn00062_"+compartment, name="ATP hydrolysis", lower_bound=0, upper_bound=1000)
+        cobra_reaction.annotation.update({"sbo":"SBO:0000176", "seed.reaction":"rxn00062"}) #biochemical reaction
         cobra_reaction.add_metabolites(stoichiometry)
         model.add_reactions([cobra_reaction])
         return {"reaction":cobra_reaction,"direction":">","new":True}
 
     @staticmethod
-    def parse_id(object):
-        if re.search('(.+)_([a-z])(\d+)$', object.id) != None:
-            m = re.search('(.+)_([a-z])(\d+)$', object.id)
+    def parse_id(cobra_obj):
+        if re.search('(.+)_([a-z])(\d+)$', cobra_obj.id):
+            m = re.search('(.+)_([a-z])(\d+)$', cobra_obj.id)
             return (m[1],m[2],int(m[3]))
         return None
 
@@ -247,11 +226,26 @@ def medianame(media):
         return media.id
 
     @staticmethod
-    def validate_dictionary(dictionary,required_keys,optional_keys):
+    def validate_dictionary(dictionary,required_keys, optional_keys={}):
         for item in required_keys:
             if item not in dictionary:
                 raise ValueError('Required key '+item+' is missing!')
         for key in optional_keys:
             if key not in dictionary:
-                dictionary[key] = defaults[key]
+                dictionary[key] = optional_keys[key]
         return dictionary
+    
+    @staticmethod
+    def get_reframed_model(kbase_model,):
+        from reframed import from_cobrapy
+        
+        reframed_model = from_cobrapy(kbase_model)
+        if hasattr(kbase_model, 'id'):
+            reframed_model.id = kbase_model.id
+        reframed_model.compartments.e0.external = True
+        return reframed_model
+    
+    @staticmethod
+    def parse_df(df):
+        from numpy import array
+        return array(dtype=object, object=[array(df.index), array(df.columns), df.to_numpy()])

modelseedpy/core/msatpcorrection.py

@@ -9,7 +9,7 @@
 from modelseedpy.core.msmodel import get_gpr_string, get_reaction_constraints_from_direction
 from cobra.core import Gene, Metabolite, Model, Reaction
 from modelseedpy.core.msmodelutl import MSModelUtil
-from modelseedpy.core import FBAHelper, MSGapfill
+from modelseedpy.core import FBAHelper, MSGapfill, MSMedia
 from modelseedpy.fbapkg.mspackagemanager import MSPackageManager
 
 logger = logging.getLogger(__name__)
@@ -43,7 +43,12 @@ def __init__(self, model, core_template, atp_medias, compartment="c0",
         else:
             output = self.modelutl.add_atp_hydrolysis(compartment)
             self.atp_hydrolysis = output["reaction"]
-        self.atp_medias = atp_medias
+        self.atp_medias = []
+        for media in atp_medias:
+            if isinstance(media,MSMedia):
+                self.atp_medias.append([media,0.01])
+            else:
+                self.atp_medias.append(media)        
         self.max_gapfilling = max_gapfilling
         self.gapfilling_delta = gapfilling_delta
         self.coretemplate = core_template
@@ -162,19 +167,21 @@ def evaluate_growth_media(self):
             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 in self.atp_medias:
+            for media_tuple in self.atp_medias:
+                media = media_tuple[0]
                 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 == 0 or solution.status != 'optimal':
-                    self.media_gapfill_stats[media] = self.msgapfill.run_gapfilling(media, self.atp_hydrolysis.id)
+                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 gapfilling fails - need to activate and penalize the noncore and try again
-                elif solution.objective_value > 0 or solution.status == 'optimal':
+                elif solution.objective_value >= media_tuple[1]:
                     self.media_gapfill_stats[media] = {'reversed': {}, 'new': {}}
+                logger.debug('gapfilling stats:',json.dumps(self.media_gapfill_stats[media],indent=2))
 
         if MSATPCorrection.DEBUG:
             with open('debug.json', 'w') as outfile: