Work on bilevel optimization and model utilities module

Author: cshenry

modelseedpy/core/msmodelutl.py

@@ -1,6 +1,7 @@
 import logging
 import re
 from cobra import Model, Reaction, Metabolite
+from modelseedpy.fbapkg.mspackagemanager import MSPackageManager
 
 logger = logging.getLogger(__name__)
 
@@ -14,6 +15,7 @@ class MSModelUtil:
 
     def __init__(self,model):
         self.model = model
+        self.pkgmgr = MSPackageManager.get_pkg_mgr(model)
         self.metabolite_hash = None
         self.search_metabolite_hash = None
 
@@ -40,6 +42,8 @@ def add_name_to_metabolite_hash(self,name,met):
         self.search_metabolite_hash[sname].append(met)
 
     def find_met(self,name):
+        if self.metabolite_hash == None:
+            self.build_metabolite_hash()
         if name in self.metabolite_hash:
             return self.metabolite_hash[name]
         sname = search_name(name)
@@ -107,4 +111,125 @@ def add_exchanges_for_metabolites(self,cpds,uptake=0,excretion=0,prefix='EX_', p
         return drains
 
     def reaction_scores(self):
-        return {}
+        return {}
+    
+    #Required this function to add gapfilled compounds to a KBase model for saving gapfilled model    
+    def convert_cobra_compound_to_kbcompound(self,cpd,kbmodel,add_to_model = 1):
+        refid = "cpd00000"
+        if re.search('cpd\d+_[a-z]+',cpd.id):
+            refid = cpd.id
+            refid = re.sub("_[a-z]\d+$","",refid)
+        cpd_data = {
+            "aliases": [],
+            "charge": cpd.charge,
+            "compound_ref": "~/template/compounds/id/"+refid,
+            "dblinks": {},
+            "formula": cpd.formula,
+            "id": cpd.id,
+            "modelcompartment_ref": "~/modelcompartments/id/"+cpd.id.split("_").pop(),
+            "name": cpd.name,
+            "numerical_attributes": {},
+            "string_attributes": {}
+        }
+        if add_to_model == 1:
+            kbmodel["modelcompounds"].append(cpd_data)
+        return cpd_data
+
+    #Required this function to add gapfilled reactions to a KBase model for saving gapfilled model    
+    def convert_cobra_reaction_to_kbreaction(self,rxn,kbmodel,cpd_hash,direction = "=",add_to_model = 1,reaction_genes = None):
+        rxnref = "~/template/reactions/id/rxn00000_c"
+        if re.search('rxn\d+_[a-z]+',rxn.id):
+            rxnref = "~/template/reactions/id/"+rxn.id
+            rxnref = re.sub("\d+$","",rxnref)
+        rxn_data = {
+            "id": rxn.id,
+            "aliases": [],
+            "dblinks": {},
+            "direction": direction,
+            "edits": {},
+            "gapfill_data": {},
+            "maxforflux": 1000000,
+            "maxrevflux": 1000000,
+            "modelReactionProteins": [],
+            "modelReactionReagents": [],
+            "modelcompartment_ref": "~/modelcompartments/id/"+rxn.id.split("_").pop(),
+            "name": rxn.name,
+            "numerical_attributes": {},
+            "probability": 0,
+            "protons": 0,
+            "reaction_ref": rxnref,
+            "string_attributes": {}
+        }
+        for cpd in rxn.metabolites:
+            if cpd.id not in kbmodel["modelcompounds"]:
+                cpd_hash[cpd.id] = self.convert_cobra_compound_to_kbcompound(cpd,kbmodel,1)
+            rxn_data["modelReactionReagents"].append({
+                "coefficient" : rxn.metabolites[cpd],
+                "modelcompound_ref" : "~/modelcompounds/id/"+cpd.id
+            })
+        if reaction_genes != None and rxn.id in reaction_genes:
+            best_gene = None
+            for gene in reaction_genes[rxn.id]:
+                if best_gene == None or reaction_genes[rxn.id][gene] > reaction_genes[rxn.id][best_gene]:
+                    best_gene = gene
+            rxn_data["modelReactionProteins"] = [{"note":"Added from gapfilling","modelReactionProteinSubunits":[],"source":"Unknown"}]
+            rxn_data["modelReactionProteins"][0]["modelReactionProteinSubunits"] = [{"note":"Added from gapfilling","optionalSubunit":0,"triggering":1,"feature_refs":["~/genome/features/id/"+best_gene],"role":"Unknown"}]
+        if add_to_model == 1:
+            kbmodel["modelreactions"].append(rxn_data)
+        return rxn_data
+    
+    def add_gapfilling_solution_to_kbase_model(self,newmodel,gapfilled_reactions,gfid=None,media_ref = None,reaction_genes = None):
+        rxn_table = []
+        gapfilling_obj = None
+        if gfid == None:
+            largest_index = 0
+            for gapfilling in newmodel["gapfillings"]:
+                current_index = int(gapfilling["id"].split(".").pop())
+                if largest_index == 0 or largest_index < current_index:
+                    largest_index = current_index
+            largest_index += 1
+            gfid = "gf."+str(largest_index)
+        else:
+            for gapfilling in newmodel["gapfillings"]:
+                if gapfilling["id"] == gfid:
+                    gapfilling_obj = gapfilling
+        if gapfilling_obj == None:    
+            gapfilling_obj = {
+                "gapfill_id": newmodel["id"]+"."+gfid,
+                "id": gfid,
+                "integrated": 1,
+                "integrated_solution": "0",
+                "media_ref": media_ref
+            }
+            newmodel["gapfillings"].append(gapfilling_obj)
+        cpd_hash = {}
+        for cpd in newmodel["modelcompounds"]:
+            cpd_hash[cpd["id"]] = cpd
+        for rxn in gapfilled_reactions["new"]:
+            reaction = self.model.reactions.get_by_id(rxn)
+            kbrxn = self.convert_cobra_reaction_to_kbreaction(reaction,newmodel,cpd_hash,gapfilled_reactions["new"][rxn],1,reaction_genes)
+            kbrxn["gapfill_data"][gfid] = dict()
+            kbrxn["gapfill_data"][gfid]["0"] = [gapfilled_reactions["new"][rxn],1,[]]
+            rxn_table.append({
+                'id':kbrxn["id"],
+                'name':kbrxn["name"],
+                'direction':format_direction(kbrxn["direction"]),
+                'gene':format_gpr(kbrxn),
+                'equation':format_equation(kbrxn,cpd_hash),
+                'newrxn':1
+            })
+        for rxn in gapfilled_reactions["reversed"]:
+            for kbrxn in newmodel["modelreactions"]:
+                if kbrxn["id"] == rxn:
+                    kbrxn["direction"] = "="
+                    rxn_table.append({
+                        'id':kbrxn["id"],
+                        'name':kbrxn["name"],
+                        'direction':format_direction(kbrxn["direction"]),
+                        'gene':format_gpr(kbrxn),
+                        'equation':format_equation(kbrxn,cpd_hash),
+                        'newrxn':0
+                    })
+                    kbrxn["gapfill_data"][gfid] = dict()
+                    kbrxn["gapfill_data"][gfid]["0"] = [gapfilled_reactions["reversed"][rxn],1,[]]
+        return rxn_table

modelseedpy/fbapkg/bilevelpkg.py

@@ -2,14 +2,15 @@
 
 from __future__ import absolute_import
 
+import re
 from optlang.symbolics import Zero, add
 from cobra.core import Gene, Metabolite, Model, Reaction
 from modelseedpy.fbapkg.basefbapkg import BaseFBAPkg
 
 #Base class for FBA packages
 class BilevelPkg(BaseFBAPkg):
     def __init__(self,model):
-        BaseFBAPkg.__init__(self,model,"reaction use",{"dualconst":"string","dualub":"string","duallb":"string"},{"dualvar":"string","objective":"string"})
+        BaseFBAPkg.__init__(self,model,"reaction use",{"dualconst":"string","dualub":"string","duallb":"string"},{"dualvar":"string","objective":"string","dualbin":"string"})
 
     def build_package(self,filter = None,binary_variable_count = 0):
         self.validate_parameters({}, [], {
@@ -38,7 +39,6 @@ def build_package(self,filter = None,binary_variable_count = 0):
         if self.parameters["binary_variable_count"] > 0:
             for reaction in self.model.reactions: 
                 var = self.build_variable("bflxcmp",reaction,None)
-                const = self.build_constraint("bflxcmp",reaction,None,None,None)                
         #Now implementing dual variables and constraints
         varhash = {}
         for var in variables:
@@ -80,12 +80,20 @@ def build_variable(self,type,object,obj_coef):
             var = BaseFBAPkg.build_variable(self,type,lb,ub,"continuous",object.name)
             obj_coef[var] = coef
             return var
-        if type == "dualub":
+        if type == "dualub":#Add a constraint that makes this variable zero when binary variable is zero
             var = BaseFBAPkg.build_variable(self,type,0,1000000,"continuous",object.name)
+            if re.search('(.+)_(fflxcmp\d+)$', object.name) is not None:
+                m = re.search('(.+)_(fflxcmp\d+)$', object.name)
+                bvar = self.variables[m[2]][m[1]]
+                BaseFBAPkg.build_constraint(self,"dualbin",None,0,{var:1,bvar:-1000000},object.name)
             obj_coef[var] = object.ub
             return var
         if type == "duallb":
             var = BaseFBAPkg.build_variable(self,type,-1000000,0,"continuous",object.name)
+            #if re.search('(.+)_(fflxcmp\d+)$', object.name) is not None:
+                #m = re.search('(.+)_(fflxcmp\d+)$', metabolite.id)
+                #bvar = self.variables[m[2]][m[1]]
+                #BaseFBAPkg.build_constraint(self,object.name+"_lbdualbin",None,0,{var:-1,bvar:-1000000},object)
             obj_coef[var] = object.lb
             return var
         if type == "flxcmp" and self.parameters["binary_variable_count"] > 0:

modelseedpy/fbapkg/gapfillingpkg.py

@@ -502,6 +502,29 @@ def run_test_conditions(self,condition_list,solution = None,max_iterations = 10)
             return None
         return solution
 
+    def filter_database_based_on_tests(self,test_conditions):
+        filetered_list = []
+        with self.model:
+            rxnlist = []
+            for reaction in self.model.reactions:
+                if reaction.id in self.gapfilling_penalties:
+                    if "reverse" in self.gapfilling_penalties[reaction.id]:
+                        reaction.lower_bound = 0
+                        rxnlist.append([reaction,"<"])
+                    if "forward" in self.gapfilling_penalties[reaction.id]:
+                        reaction.upper_bound = 0
+                        rxnlist.append([reaction,">"])
+                    reaction.update_variable_bounds()
+            self.pkgmgr.getpkg("ObjConstPkg").constraints["objc"]["1"].lb = 0
+            filtered_list = FBAHelper.reaction_expansion_test(self.model,rxnlist,test_conditions,self.pkgmgr)
+        #Now constraining filtered reactions to zero
+        for item in filtered_list:
+            print("Filtering:",item[0].id,item[1])
+            if item[1] == ">":
+                self.model.reactions.get_by_id(item[0].id).upper_bound = 0
+            else:
+                self.model.reactions.get_by_id(item[0].id).lower_bound = 0
+    
     def compute_gapfilled_solution(self, flux_values=None):
         if flux_values is None:
             flux_values = FBAHelper.compute_flux_values_from_variables(self.model)

modelseedpy/multiomics/msexpression.py

@@ -55,10 +55,11 @@ def __init__(self, id):
 
 class MSExpressionFeature:
 
-    def __init__(self, feature):
+    def __init__(self,feature,parent):
         self.id = feature.id
         self.feature = feature
         self.values = {}
+        self.parent = parent
 
     def add_value(self,condition,value):
         if condition in self.values:
@@ -72,6 +73,11 @@ def add_value(self,condition,value):
         self.values[condition] = value
 
     def get_value(self,condition,normalization = None):
+        if isinstance(condition,str):
+            if condition not in self.parent.conditions:
+                logger.warning("Condition "+condition+" not found in expression object!")
+                return None
+            condition = self.parent.conditions.get_by_id(condition)
         if condition not in self.values:
             logger.info("Condition "+condition.id+" has no value in "+self.feature.id)
             return None
@@ -133,11 +139,11 @@ def add_feature(self,id,create_gene_if_missing = False):
             if id in self.object.reactions:
                 feature = self.object.reactions.get_by_id(id)
         if feature == None:
-            logger.warning("Gene referred by expression "+id+" not found in genome object!")
+            logger.warning("Feature referred by expression "+id+" not found in genome object!")
             return None
         if feature.id in self.features:
             return self.features.get_by_id(feature.id)
-        protfeature = MSExpressionFeature(feature)
+        protfeature = MSExpressionFeature(feature,self)
         self.features.append(protfeature)
         return protfeature
 
@@ -147,11 +153,6 @@ def get_value(self,feature,condition,normalization = None):
                 logger.warning("Feature "+feature+" not found in expression object!")
                 return None
             feature = self.features.get_by_id(feature)
-        if isinstance(condition,str):
-            if condition not in self.conditions:
-                logger.warning("Condition "+condition+" not found in expression object!")
-                return None
-            condition = self.conditions.get_by_id(condition)
         return feature.get_value(condition,normalization)
 
     def build_reaction_expression(self,model,default):