Updating and improving code for growth phenotype analysis and omics integration

Author: cshenry

modelseedpy/__init__.py

@@ -48,7 +48,7 @@
 from modelseedpy.fbapkg import (
     BaseFBAPkg,RevBinPkg,ReactionUsePkg,SimpleThermoPkg,TotalFluxPkg,ElementUptakePkg,BilevelPkg,
     CommKineticPkg,KBaseMediaPkg,FluxFittingPkg,ProteomeFittingPkg,GapfillingPkg,MetaboFBAPkg,FlexibleBiomassPkg,
-    ProblemReplicationPkg,FullThermoPkg,MSPackageManager,ObjConstPkg
+    ProblemReplicationPkg,FullThermoPkg,MSPackageManager,ObjConstPkg,ChangeOptPkg
 )
 
 from modelseedpy.multiomics import (

modelseedpy/core/msgenome.py

@@ -38,7 +38,7 @@ def parse_fasta_str(faa_str, split='|', h_func=None):
             elif split:
                 header_data = line[1:].split(split, 1)
                 seq_id = header_data[0]
-                desc = header_data[1]
+                #desc = header_data[1]
 
             seq = MSFeature(seq_id, "", desc)
         else:

modelseedpy/core/msgrowthphenotypes.py

@@ -1,5 +1,6 @@
 import pandas as pd
 import logging
+import cobra
 from cobra.core.dictlist import DictList
 from modelseedpy.core.msmedia import MSMedia
 from modelseedpy.fbapkg.mspackagemanager import MSPackageManager
@@ -33,7 +34,9 @@ def build_media(self):
             full_media.merge(parent.base_media,overwrite_overlap = False)
         return full_media
 
-    def simulate(self,modelutl,growth_threshold=0.001,add_missing_exchanges=False,save_fluxes=False):        
+    def simulate(self,modelutl,growth_threshold=0.001,add_missing_exchanges=False,save_fluxes=False,pfba=False):        
+        if not isinstance(modelutl,MSModelUtil):
+            modelutl = MSModelUtil(modelutl)
         media = self.build_media()
         output = {"growth":None,"class":None,"missing_transports":[]}
         if add_missing_exchanges:
@@ -44,11 +47,13 @@ def simulate(self,modelutl,growth_threshold=0.001,add_missing_exchanges=False,sa
             if gene in modelutl.model.genes:
                 geneobj = modelutl.model.genes.get_by_id(gene)
                 geneobj.knock_out()
-        solution = modelutl.model.optimize()        
+        solution = modelutl.model.optimize()
+        output["growth"] = solution.objective_value
+        if solution.objective_value > 0 and pfba:
+            solution = cobra.flux_analysis.pfba(modelutl.model)
         if save_fluxes:
             output["fluxes"] = solution.fluxes
-        output["growth"] = solution.objective_value
-        if solution.objective_value >= growth_threshold:
+        if output["growth"] >= growth_threshold:
             if self.growth > 0:
                 output["class"] = "CP"
             else:
@@ -61,6 +66,8 @@ def simulate(self,modelutl,growth_threshold=0.001,add_missing_exchanges=False,sa
         return output        
 
     def gapfill_model_for_phenotype(self,modelutl,default_gapfill_templates,test_conditions,default_gapfill_models=[],blacklist=[],growth_threshold=0.001,add_missing_exchanges=False):
+        if not isinstance(modelutl,MSModelUtil):
+            modelutl = MSModelUtil(modelutl)
         self.gapfilling = MSGapfill(modelutl.model, default_gapfill_templates, default_gapfill_models, test_conditions, modelutl.reaction_scores(), blacklist) 
         media = self.build_media()
         if add_missing_exchanges:
@@ -173,9 +180,7 @@ def simulate_phenotypes(self,model,biomass,add_missing_exchanges=False,correct_f
         model.objective = biomass
         modelutl = MSModelUtil(model)
         summary = {"Label":["Accuracy","CP","CN","FP","FN"],"Count":[0,0,0,0,0]}
-        data = {"Phenotype":[],"Observed growth":[],"Simulated growth":[],"Class":[],"Transports missing":[]}
-        if correct_false_negatives:
-            data["Gapfilled reactions"] = []
+        data = {"Phenotype":[],"Observed growth":[],"Simulated growth":[],"Class":[],"Transports missing":[],"Gapfilled reactions":[]}
         for pheno in self.phenotypes:
             with model:
                 result = pheno.simulate(modelutl,growth_threshold,add_missing_exchanges,save_fluxes)#Result should have "growth" and "class"

modelseedpy/core/msmedia.py

@@ -32,7 +32,7 @@ def __init__(self, media_id):
     def from_dict(d):
         """
         Either dict with exchange bounds (example: {'cpd00027': (-10, 1000)}) or
-        just absolute value of update (example: {''cpd00027': 10})
+        just absolute value of uptake (example: {''cpd00027': 10})
         :param d:
         :return:
         """

modelseedpy/core/msmodelutl.py

@@ -6,7 +6,7 @@
 
 def search_name(name):
     name = name.lower()
-    name = re.sub(r'_[a-z]\d+$', '', name)
+    name = re.sub(r'_[a-z]\d*$', '', name)
     name = re.sub(r'\W+', '', name)
     return name
 
@@ -24,7 +24,11 @@ def build_metabolite_hash(self):
             self.add_name_to_metabolite_hash(met.id,met)
             self.add_name_to_metabolite_hash(met.name,met)
             for anno in met.annotation:
-                self.add_name_to_metabolite_hash(met.annotation[anno],met)
+                if isinstance(met.annotation[anno], list):
+                    for item in met.annotation[anno]:
+                        self.add_name_to_metabolite_hash(item,met)
+                else:
+                    self.add_name_to_metabolite_hash(met.annotation[anno],met)
 
     def add_name_to_metabolite_hash(self,name,met):
         if name not in self.metabolite_hash:
@@ -100,6 +104,7 @@ def add_exchanges_for_metabolites(self,cpds,uptake=0,excretion=0,prefix='EX_', p
             drain_reaction.annotation["sbo"] = 'SBO:0000627'    
             drains.append(drain_reaction)
         self.model.add_reactions(drains)
+        return drains
 
     def reaction_scores(self):
         return {}

modelseedpy/fbapkg/__init__.py

@@ -20,3 +20,4 @@
 from modelseedpy.fbapkg.problemreplicationpkg import ProblemReplicationPkg
 from modelseedpy.fbapkg.fullthermopkg import FullThermoPkg
 from modelseedpy.fbapkg.objconstpkg import ObjConstPkg
+from modelseedpy.fbapkg.changeoptpkg import ChangeOptPkg

modelseedpy/fbapkg/basefbapkg.py

@@ -100,3 +100,15 @@ def all_variables(self):
 
     def all_constraints(self):
         return self.pkgmgr.all_constraints()
+    
+    def add_variable_type(self,name,type):
+        if name not in self.variables:
+            self.variables[name] = dict()
+        if name not in self.variable_types:
+            self.variable_types[name] = type
+            
+    def add_constraint_type(self,name,type):
+        if name not in self.constraints:
+            self.constraints[name] = dict()
+        if name not in self.constraint_types:
+            self.constraint_types[name] = type

modelseedpy/fbapkg/bilevelpkg.py

@@ -12,21 +12,21 @@ def __init__(self,model):
         BaseFBAPkg.__init__(self,model,"reaction use",{"dualconst":"string","dualub":"string","duallb":"string"},{"dualvar":"string","objective":"string"})
 
     def build_package(self,filter = None,binary_variable_count = 0):
-        self.validate_parameters(parameters, [], {
+        self.validate_parameters({}, [], {
             "binary_variable_count":binary_variable_count
         });
+        print("binary_variable_count:",binary_variable_count)
         coefficients = {}
         obj_coef = {}
         obj = self.model.solver.objective
-        #Using the JSON calls because get_linear_coefficients is REALLY slow
         #Creating new objective coefficient and bound variables
         bound_variables = {}
         reactions = self.model.reactions
-        if binary_variable_count > 0:
+        if self.parameters["binary_variable_count"] > 0:
             for reaction in self.model.reactions: 
-                var = self.build_variable("fluxcomp",reaction,None)
-                const = self.build_constraint("fluxcomp",reaction,None,None,None)
+                var = self.build_variable("flxcmp",reaction,None)
         #Retreiving model data with componenent flux variables
+        #Using the JSON calls because get_linear_coefficients is REALLY slow
         mdldata = self.model.solver.to_json()
         consthash = {}
         for const in mdldata["constraints"]:
@@ -35,10 +35,10 @@ def build_package(self,filter = None,binary_variable_count = 0):
         variables = list(self.model.solver.variables)
         objterms = obj.get_linear_coefficients(variables)
         #Adding binary variables and constraints which should not be included in dual formulation
-        if binary_variable_count > 0:
-            for variable in self.variables["fluxcomp"]: 
-                var = self.build_variable("bfluxcomp",variable,None)
-                const = self.build_constraint("bfluxcomp",variable,None,None,None)                
+        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:
@@ -65,7 +65,6 @@ def build_package(self,filter = None,binary_variable_count = 0):
                     coefficients[var][dvar] = 1
                 self.build_constraint("dualvar",var,obj,objterms,coefficients)
         self.build_constraint("objective",None,obj,objterms,obj_coef)
-        #Add binary variables and constaints which should not be part of the formulation
 
     def build_variable(self,type,object,obj_coef):
         if type == "dualconst":
@@ -95,30 +94,38 @@ def build_variable(self,type,object,obj_coef):
             for i in range(0,self.parameters["binary_variable_count"]):
                 value = 2**i
                 if object.lower_bound < 0:
-                    var = BaseFBAPkg.build_variable(self,"rflxcmp"+str(i),0,-1*value*object.lower_bound/denominator,"continuous",object.id)
+                    self.add_variable_type("rflxcmp"+str(i),"reaction")
+                    var = BaseFBAPkg.build_variable(self,"rflxcmp"+str(i),0,-1*value*object.lower_bound/denominator,"continuous",object)
                     coefs[0][var] = -1
                 if object.upper_bound > 0:
-                    var = BaseFBAPkg.build_variable(self,"fflxcmp"+str(i),0,value*object.upper_bound/denominator,"continuous",object.id)
+                    self.add_variable_type("fflxcmp"+str(i),"reaction")
+                    var = BaseFBAPkg.build_variable(self,"fflxcmp"+str(i),0,value*object.upper_bound/denominator,"continuous",object)
                     coefs[1][var] = -1
             if object.lower_bound < 0:
                 #flux - flux_comp_0 - flux_comp_n = 0 - restriction of reverse fluxes by component fluxes
+                self.add_constraint_type("rflxcmpc","reaction")
                 coefs[0][object.reverse_variable] = 1
-                BaseFBAPkg.build_constraint(self,"rflxcmpc",0,0,coefs[0],object.id)
+                BaseFBAPkg.build_constraint(self,"rflxcmpc",0,0,coefs[0],object)
             if object.upper_bound > 0:
                 #flux - flux_comp_0 - flux_comp_n = 0 - restriction of forward fluxes by component fluxes
+                self.add_constraint_type("fflxcmpc","reaction")
                 coefs[1][object.forward_variable] = 1
-                BaseFBAPkg.build_constraint(self,"fflxcmpc",0,0,coefs[1],object.id)
+                BaseFBAPkg.build_constraint(self,"fflxcmpc",0,0,coefs[1],object)
             return None
         if type == "bflxcmp" and self.parameters["binary_variable_count"] > 0:
             for i in range(0,self.parameters["binary_variable_count"]):
                 if object.lower_bound < 0:
-                    var = BaseFBAPkg.build_variable(self,"brflxcmp"+str(i),0,1,"binary",object.id)
+                    self.add_variable_type("brflxcmp"+str(i),"reaction")
+                    var = BaseFBAPkg.build_variable(self,"brflxcmp"+str(i),0,1,"binary",object)
                     othervar = self.variables["rflxcmp"+str(i)][object.id]
-                    BaseFBAPkg.build_constraint(self,"brflxcmpc"+str(i),None,0,{othervar:1,var:-1000},object.id)
+                    self.add_constraint_type("brflxcmpc"+str(i),"reaction")
+                    BaseFBAPkg.build_constraint(self,"brflxcmpc"+str(i),None,0,{othervar:1,var:-1000},object)
                 if object.upper_bound > 0:
-                    var = BaseFBAPkg.build_variable(self,"bfflxcmp"+str(i),0,1,"binary",object.id)
+                    self.add_variable_type("bfflxcmp"+str(i),"reaction")
+                    var = BaseFBAPkg.build_variable(self,"bfflxcmp"+str(i),0,1,"binary",object)
                     othervar = self.variables["fflxcmp"+str(i)][object.id]
-                    BaseFBAPkg.build_constraint(self,"bfflxcmpc"+str(i),None,0,{othervar:1,var:-1000},object.id)
+                    self.add_constraint_type("bfflxcmpc"+str(i),"reaction")
+                    BaseFBAPkg.build_constraint(self,"bfflxcmpc"+str(i),None,0,{othervar:1,var:-1000},object)
             return None
 
     def build_constraint(self,type,object,objective,objterms,coefficients):

modelseedpy/fbapkg/changeoptpkg.py

@@ -0,0 +1,59 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import absolute_import
+
+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 ChangeOptPkg(BaseFBAPkg):
+    def __init__(self,model):
+        BaseFBAPkg.__init__(self,model,"change optimization",{"bgoal":"string"},{"bgoalc":"reaction"})
+    
+    def build_package(self,target_values = {},build_objective = True):
+        self.validate_parameters({}, [], {
+            "target_values":target_values
+        });
+        #Creating objective but won't set objective unless build_objective is true
+        goal_objective = self.model.problem.Objective(Zero,direction="max")
+        obj_coef = dict()
+        for rxnid in target_values:
+            print(rxnid)
+            if rxnid in self.model.reactions:
+                print("FOUND!")
+                rxn = self.model.reactions.get_by_id(rxnid)
+                var = self.build_variable("bgoal",rxn)
+                obj_coef[var] = target_values[rxnid]["objcoef"]
+                const = self.build_constraint("bgoalc",rxn)    
+        #Now setting the goal objective if build_objective is true
+        if build_objective:
+            self.model.objective = goal_objective
+            goal_objective.set_linear_coefficients(obj_coef)
+    
+    def build_variable(self,type,object):
+        if type == "bgoal":
+            goal = self.parameters["target_values"][object.id]
+            return BaseFBAPkg.build_variable(self,type,0,1,"binary",object.id+goal["direction"])
+                    
+    def build_constraint(self,type,object):
+        if type == "bgoalc":
+            #For lower: goal - vi + 1000 - 1000 gi >= 0
+            #For higher: vi - goal  + 1000 - 1000 gi >= 0
+            lb = -1000
+            goal = self.parameters["target_values"][object.id]
+            var = self.variables["bgoal"][object.id+goal["direction"]]
+            fluxvar = None
+            if goal["flux"] < 0:
+                fluxvar = object.reverse_variable
+            else:
+                fluxvar = object.forward_variable
+            coef = {var:-1000}
+            if goal["direction"] == "low":
+                coef[fluxvar] = -1
+                lb += -1*abs(goal["flux"])
+            else:
+                coef[fluxvar] = 1
+                lb += abs(goal["flux"])
+            return BaseFBAPkg.build_constraint(self,type,lb,None,coef,object)
+