diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index 5a9bd5f5..547a436f 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -3,12 +3,14 @@ from __future__ import absolute_import import logging + import re -import json -from optlang.symbolics import Zero, add -from cobra import Model, Reaction, Metabolite +import json # !!! import is never used +from optlang.symbolics import Zero, add # !!! add is never used +from cobra import Model, Reaction, Metabolite # !!! Model is never used from modelseedpy.fbapkg.basefbapkg import BaseFBAPkg from modelseedpy.core.fbahelper import FBAHelper +from modelseedpy.core.template import Template # !!! possibly use the conversion functions from this class to replace the conversion functions defined herein? logger = logging.getLogger(__name__) @@ -67,10 +69,10 @@ class GapfillingPkg(BaseFBAPkg): """ """ - def __init__(self, model): BaseFBAPkg.__init__(self, model, "gapfilling", {}, {}) self.gapfilling_penalties = None + self.tempalte = Template() def build(self, template, minimum_objective=0.01): parameters = { @@ -137,31 +139,25 @@ def build_package(self, parameters): # Iterating over all indecies with more than 10 intracellular compounds: self.gapfilling_penalties = dict() - for index in indexhash: - if indexhash[index] > 10: + for index, val in indexhash.items(): + if val > 10: if index == "none": for template in self.parameters["default_gapfill_templates"]: - new_penalties = self.extend_model_with_template_for_gapfilling(template, index) - self.gapfilling_penalties.update(new_penalties) + self.gapfilling_penalties.update(self.extend_model_with_template_for_gapfilling(template, index)) for gfmdl in self.parameters["default_gapfill_models"]: - new_penalties = self.extend_model_with_model_for_gapfilling(gfmdl, index) - self.gapfilling_penalties.update(new_penalties) + self.gapfilling_penalties.update(self.extend_model_with_model_for_gapfilling(gfmdl, index)) if index in self.parameters["gapfill_templates_by_index"]: for template in self.parameters["gapfill_templates_by_index"][index]: - new_penalties = self.extend_model_with_template_for_gapfilling(template, index) - self.gapfilling_penalties.update(new_penalties) + self.gapfilling_penalties.update(self.extend_model_with_template_for_gapfilling(template, index)) if index in self.parameters["gapfill_models_by_index"]: for gfmdl in self.parameters["gapfill_models_by_index"]: - new_penalties = self.extend_model_with_model_for_gapfilling(gfmdl, index) - self.gapfilling_penalties.update(new_penalties) + self.gapfilling_penalties.update(self.extend_model_with_model_for_gapfilling(gfmdl, index)) if self.parameters["gapfill_all_indecies_with_default_templates"]: for template in self.parameters["default_gapfill_templates"]: - new_penalties = self.extend_model_with_template_for_gapfilling(template, index) - self.gapfilling_penalties.update(new_penalties) + self.gapfilling_penalties.update(self.extend_model_with_template_for_gapfilling(template, index)) if self.parameters["gapfill_all_indecies_with_default_models"]: for gfmdl in self.parameters["default_gapfill_models"]: - new_penalties = self.extend_model_with_model_for_gapfilling(gfmdl, index) - self.gapfilling_penalties.update(new_penalties) + self.gapfilling_penalties.update(self.extend_model_with_model_for_gapfilling(gfmdl, index)) # Rescaling penalties by reaction scores and saving genes for reaction in self.gapfilling_penalties: rxnid = reaction.split("_")[0] @@ -178,9 +174,7 @@ def build_package(self, parameters): self.model.solver.update() if self.parameters["set_objective"] == 1: - reaction_objective = self.model.problem.Objective( - Zero, - direction="min") + reaction_objective = self.model.problem.Objective(Zero, direction="min") obj_coef = dict() for reaction in self.model.reactions: if reaction.id in self.gapfilling_penalties: @@ -194,22 +188,17 @@ def build_package(self, parameters): # elif default_penalty != 0: # obj_coef[reaction.forward_variable] = 0 else: - obj_coef[reaction.forward_variable] = 0 - obj_coef[reaction.reverse_variable] = 0 + obj_coef[reaction.forward_variable] = obj_coef[reaction.reverse_variable] = 0 self.model.objective = reaction_objective reaction_objective.set_linear_coefficients(obj_coef) def extend_model_with_model_for_gapfilling(self, source_model, index): - new_metabolites = {} - new_reactions = {} - new_exchange = [] - new_demand = [] - new_penalties = dict() - local_remap = {} + self.new_metabolites, self.new_reactions, local_remap, new_penalties = {}, {}, {}, {} + new_exchange, new_demand = [], [] # Adding metabolites from source model to gapfill model for cobra_metabolite in source_model.metabolites: original_id = cobra_metabolite.id - if re.search('(.+)_([a-z])\d+$', cobra_metabolite.id) != None: + if re.search('(.+)_([a-z])\d+$', cobra_metabolite.id): m = re.search('(.+)_([a-z])\d+$', cobra_metabolite.id) if m[2] == "e": cobra_metabolite.compartment = "e0" @@ -217,26 +206,26 @@ def extend_model_with_model_for_gapfilling(self, source_model, index): else: cobra_metabolite.compartment = m[2] + index cobra_metabolite.id = m[1] + "_" + m[2] + index - if cobra_metabolite.id not in self.model.metabolites and cobra_metabolite.id not in new_metabolites: - new_metabolites[cobra_metabolite.id] = cobra_metabolite + if cobra_metabolite.id not in self.model.metabolites and cobra_metabolite.id not in self.new_metabolites: + self.new_metabolites[cobra_metabolite.id] = cobra_metabolite local_remap[original_id] = cobra_metabolite if m[1] + "_" + m[2] in self.parameters["auto_sink"]: new_demand.append(cobra_metabolite) if m[2] == "e": new_exchange.append(cobra_metabolite) # Adding all metabolites to model prior to adding reactions - self.model.add_metabolites(new_metabolites.values()) + self.model.add_metabolites(self.new_metabolites.values()) # Adding reactions from source model to gapfill model for modelreaction in source_model.reactions: if re.search('(.+)_([a-z])\d+$', modelreaction.id) != None: m = re.search('(.+)_([a-z])\d+$', modelreaction.id) if m[1] not in self.parameters["blacklist"]: cobra_reaction = modelreaction.copy() - cobra_reaction.id = groups[1] + "_" + groups[2] + index - if cobra_reaction.id not in self.model.reactions and cobra_reaction.id not in new_reactions: - new_reactions[cobra_reaction.id] = cobra_reaction - new_penalties[cobra_reaction.id] = dict(); - new_penalties[cobra_reaction.id]["added"] = 1 + cobra_reaction.id = m[1] + "_" + m[2] + index + if cobra_reaction.id not in self.model.reactions and cobra_reaction.id not in self.new_reactions: + self.new_reactions[cobra_reaction.id] = cobra_reaction + new_penalties[cobra_reaction.id] = {} + new_penalties[cobra_reaction.id]["added"] = True if cobra_reaction.lower_bound < 0: new_penalties[cobra_reaction.id]["reverse"] = self.parameters["model_penalty"] if cobra_reaction.upper_bound > 0: @@ -251,25 +240,23 @@ def extend_model_with_model_for_gapfilling(self, source_model, index): if local_remap[metabolite.id] != metabolite: new_stoichiometry[metabolite] = 0 cobra_reaction.add_metabolites(new_stoichiometry, combine=False) - elif cobra_reaction.lower_bound < 0 and self.model.reactions.get_by_id( - cobra_reaction.id).lower_bound == 0: + elif cobra_reaction.lower_bound < 0 and self.model.reactions.get_by_id(cobra_reaction.id).lower_bound == 0: self.model.reactions.get_by_id(cobra_reaction.id).lower_bound = cobra_reaction.lower_bound self.model.reactions.get_by_id(cobra_reaction.id).update_variable_bounds() new_penalties[cobra_reaction.id]["reverse"] = self.parameters["model_penalty"] - new_penalties[cobra_reaction.id]["reversed"] = 1 - elif cobra_reaction.upper_bound > 0 and self.model.reactions.get_by_id( - cobra_reaction.id).upper_bound == 0: + new_penalties[cobra_reaction.id]["reversed"] = True + elif cobra_reaction.upper_bound > 0 and self.model.reactions.get_by_id(cobra_reaction.id).upper_bound == 0: self.model.reactions.get_by_id(cobra_reaction.id).upper_bound = cobra_reaction.upper_bound self.model.reactions.get_by_id(cobra_reaction.id).update_variable_bounds() - new_penalties[cobra_reaction.id]["forward"] = model_penalty - new_penalties[cobra_reaction.id]["reversed"] = 1 + new_penalties[cobra_reaction.id]["forward"] = self.parameters["model_penalty"] + new_penalties[cobra_reaction.id]["reversed"] = True - # Only run this on new exchanges so we don't readd for all exchanges + # Only run this on new exchanges so we don't readd for all exchanges self.modelutl.add_exchanges_for_metabolites(new_exchange,self.parameters["default_uptake"],self.parameters["default_excretion"]) # Only run this on new demands so we don't readd for all exchanges self.modelutl.add_exchanges_for_metabolites(new_demand,self.parameters["default_uptake"],self.parameters["default_excretion"],"DM_") # Adding all new reactions to the model at once (much faster than one at a time) - self.model.add_reactions(new_reactions.values()) + self.model.add_reactions(self.new_reactions.values()) return new_penalties def extend_model_with_template_metabolites(self, template, index='0'): @@ -279,54 +266,51 @@ def extend_model_with_template_metabolites(self, template, index='0'): :param index: :return: """ - new_metabolites = {} - new_exchange = [] - new_demand = [] + self.new_metabolites = {} + new_exchange, new_demand = [], [] for template_compound in template.compcompounds: compartment = template_compound.compartment compartment_index = "0" if compartment == 'e' else index - cobra_metabolite = self.convert_template_compound(template_compound, compartment_index, template) # TODO: move function out - if cobra_metabolite.id not in self.model.metabolites and cobra_metabolite.id not in new_metabolites: - new_metabolites[cobra_metabolite.id] = cobra_metabolite + cobra_metabolite = self.convert_template_compound(template_compound, compartment_index, template) # TODO: move function out # !!! core.template has the same function? + if cobra_metabolite.id not in self.model.metabolites and cobra_metabolite.id not in self.new_metabolites: + self.new_metabolites[cobra_metabolite.id] = cobra_metabolite #self.model.add_metabolites([cobra_metabolite]) - msid = FBAHelper.modelseed_id_from_cobra_metabolite(cobra_metabolite) + msid = FBAHelper.modelseed_id_from_cobra_metaboliteabolite(cobra_metabolite) if msid in self.parameters["auto_sink"]: if msid != "cpd11416" or cobra_metabolite.compartment == "c0": new_demand.append(cobra_metabolite) if compartment == "e": new_exchange.append(cobra_metabolite) # Adding all metabolites to model prior to adding reactions - self.model.add_metabolites(new_metabolites.values()) + self.model.add_metabolites(self.new_metabolites.values()) return new_exchange, new_demand # Possible new function to add to the KBaseFBAModelToCobraBuilder to extend a model with a template for gapfilling for a specific index def extend_model_with_template_for_gapfilling(self, template, index): - new_reactions = {} - new_penalties = dict() - + self.new_reactions, new_penalties = {}, {} # Adding all metabolites to model prior to adding reactions new_exchange, new_demand = self.extend_model_with_template_metabolites(template, index) for template_reaction in template.reactions: if template_reaction.reference_id in self.parameters["blacklist"]: continue - cobra_reaction = self.convert_template_reaction(template_reaction, index, template, 1) # TODO: move function out + cobra_reaction = self.convert_template_reaction(template_reaction, index, template, True) # TODO: move function out new_penalties[cobra_reaction.id] = dict() - if cobra_reaction.id not in self.model.reactions and cobra_reaction.id not in new_reactions: + if cobra_reaction.id not in self.model.reactions and cobra_reaction.id not in self.new_reactions: # Adding any template reactions missing from the present model - new_reactions[cobra_reaction.id] = cobra_reaction + self.new_reactions[cobra_reaction.id] = cobra_reaction if cobra_reaction.lower_bound < 0: new_penalties[cobra_reaction.id][ "reverse"] = template_reaction.base_cost + template_reaction.reverse_penalty if cobra_reaction.upper_bound > 0: new_penalties[cobra_reaction.id][ "forward"] = template_reaction.base_cost + template_reaction.forward_penalty - new_penalties[cobra_reaction.id]["added"] = 1 + new_penalties[cobra_reaction.id]["added"] = True elif template_reaction.GapfillDirection == "=": # Adjusting directionality as needed for existing reactions model_reaction = self.model.reactions.get_by_id(cobra_reaction.id) - new_penalties[cobra_reaction.id]["reversed"] = 1 + new_penalties[cobra_reaction.id]["reversed"] = True if model_reaction.lower_bound == 0: model_reaction.lower_bound = template_reaction.lower_bound model_reaction.update_variable_bounds() @@ -341,65 +325,49 @@ def extend_model_with_template_for_gapfilling(self, template, index): exchanges = self.modelutl.add_exchanges_for_metabolites(new_exchange,self.parameters["default_uptake"],self.parameters["default_excretion"]) for ex in exchanges: new_penalties[ex.id] = { - 'added': 1, - 'reverse': 1, - 'forward': 1 - } + 'added': True, + 'reverse': True, + 'forward': True + } # Only run this on new demands so we don't readd for all exchanges exchanges = self.modelutl.add_exchanges_for_metabolites(new_demand,self.parameters["default_uptake"],self.parameters["default_excretion"],"DM_") for ex in exchanges: new_penalties[ex.id] = { - 'added': 1, - 'reverse': 1, - 'forward': 1 - } + 'added': True, + 'reverse': True, + 'forward': True + } # Adding all new reactions to the model at once (much faster than one at a time) - self.model.add_reactions(new_reactions.values()) + self.model.add_reactions(self.new_reactions.values()) return new_penalties def convert_template_compound(self, template_compound, index, template): base_id = template_compound.id.split("_")[0] base_compound = template.compounds.get_by_id(base_id) - new_id = template_compound.id - new_id += str(index) - compartment = template_compound.compartment - compartment += str(index) - - met = Metabolite(new_id, - formula=base_compound.formula, - name=base_compound.name, - charge=template_compound.charge, - compartment=compartment) + new_id = template_compound.id+str(index) + compartment = template_compound.compartment+str(index) + met = Metabolite(new_id, formula=base_compound.formula, name=base_compound.name, + charge=template_compound.charge, compartment=compartment) met.annotation["sbo"] = "SBO:0000247" # simple chemical - Simple, non-repetitive chemical entity. met.annotation["seed.compound"] = base_id return met - def convert_template_reaction(self, template_reaction, index, template, for_gapfilling=1): - array = template_reaction.id.split("_") - base_id = array[0] - new_id = template_reaction.id - new_id += str(index) - + def convert_template_reaction(self, template_reaction, index, template, for_gapfilling=True): # !!! this function is also in core.template?? + base_id = template_reaction.id.split("_")[0] # !!! base_id is never used + new_id = template_reaction.id+str(index) lower_bound = template_reaction.lower_bound upper_bound = template_reaction.upper_bound - direction = template_reaction.GapfillDirection - if for_gapfilling == 0: + if not for_gapfilling: direction = template_reaction.direction - if direction == ">": lower_bound = 0 elif direction == "<": upper_bound = 0 - cobra_reaction = Reaction(new_id, - name=template_reaction.name, - lower_bound=lower_bound, - upper_bound=upper_bound) - object_stoichiometry = {} for m, value in template_reaction.metabolites.items(): metabolite_id = m.id @@ -412,38 +380,31 @@ def convert_template_reaction(self, template_reaction, index, template, for_gapf metabolite = self.model.metabolites.get_by_id(metabolite_id) object_stoichiometry[metabolite] = value + cobra_reaction = Reaction(new_id, name=template_reaction.name, lower_bound=lower_bound, upper_bound=upper_bound) cobra_reaction.add_metabolites(object_stoichiometry) - cobra_reaction.annotation["sbo"] = "SBO:0000176" # biochemical reaction cobra_reaction.annotation["seed.reaction"] = template_reaction.reference_id - return cobra_reaction def binary_check_gapfilling_solution(self, solution=None, flux_values=None): - if solution is None: - solution = self.compute_gapfilled_solution() - if flux_values is None: - flux_values = self.modelutl.compute_flux_values_from_variables() - filter = {} - for rxn_id in solution["reversed"]: - filter[rxn_id] = solution["reversed"][rxn_id] - for rxn_id in solution["new"]: - filter[rxn_id] = solution["new"][rxn_id] - self.pkgmgr.getpkg("ReactionUsePkg").build_package(filter) + solution = solution or self.compute_gapfilled_solution(flux_values) + flux_values = flux_values or FBAHelper.compute_flux_values_from_variables(self.model) + rxn_filter = {rxn_id:solution["reversed"][rxn_id] for rxn_id in solution["reversed"]} + rxn_filter.update({rxn_id:solution["new"][rxn_id] for rxn_id in solution["new"]}) + self.pkgmgr.getpkg("ReactionUsePkg").build_package(rxn_filter) objcoef = {} - for rxnid in filter: - if filter[rxnid] == ">": + for rxnid in rxn_filter: + if rxn_filter[rxnid] == ">": objcoef[self.pkgmgr.getpkg("ReactionUsePkg").variables["fu"][rxnid]] = 1 - if filter[rxnid] == "<": + if rxn_filter[rxnid] == "<": objcoef[self.pkgmgr.getpkg("ReactionUsePkg").variables["ru"][rxnid]] = 1 new_solution = {} - with self.model: + with self.model: # to prevent the model for permanently assuming the zeroed reactions # Setting all gapfilled reactions not in the solution to zero self.knockout_gf_reactions_outside_solution(solution,flux_values) - # Setting the objective to be minimization of sum of binary variables - min_reaction_objective = self.model.problem.Objective(Zero, direction="min") - self.model.objective = min_reaction_objective - min_reaction_objective.set_linear_coefficients(objcoef) + # Setting the objective to the minimum sum of binary variables + self.model.objective = self.model.problem.Objective(Zero, direction="min") + self.model.objective.set_linear_coefficients(objcoef) self.model.optimize() new_solution = self.compute_gapfilled_solution() return new_solution @@ -451,7 +412,7 @@ def binary_check_gapfilling_solution(self, solution=None, flux_values=None): #This function is designed to KO all gapfilled reactions not included in the solution def knockout_gf_reactions_outside_solution(self,solution = None,flux_values = None): if solution == None: - solution = self.compute_gapfilled_solution() + solution = self.compute_gapfilled_solution() # !!! should flux_values be added as an argument here? if flux_values == None: flux_values = self.modelutl.compute_flux_values_from_variables() for rxnobj in self.model.reactions: @@ -462,16 +423,15 @@ def knockout_gf_reactions_outside_solution(self,solution = None,flux_values = No rxnobj.upper_bound = 0 rxnobj.update_variable_bounds() - def run_test_conditions(self,condition_list,solution = None,max_iterations = 10): + def run_test_conditions(self, condition_list, solution = None, max_iterations = 10): + reaction_list, filtered_list = [], [] if solution == None: solution = self.compute_gapfilled_solution() - reaction_list = [] for rxnid in solution["reversed"]: reaction_list.append([self.model.reactions.get_by_id(rxnid),solution["reversed"][rxnid]]) for rxnid in solution["new"]: reaction_list.append([self.model.reactions.get_by_id(rxnid),solution["new"][rxnid]]) - filtered_list = [] - with self.model: + with self.model: # to prevent the model for permanently assuming the zeroed reactions #Setting all gapfilled reactions not in the solution to zero self.knockout_gf_reactions_outside_solution(solution) self.pkgmgr.getpkg("ObjConstPkg").constraints["objc"]["1"].lb = 0 @@ -498,7 +458,7 @@ def run_test_conditions(self,condition_list,solution = None,max_iterations = 10) return solution def filter_database_based_on_tests(self,test_conditions): - filetered_list = [] + filtered_list = [] with self.model: rxnlist = [] for reaction in self.model.reactions: