From b0129609d6840cc4861a02313a8730e3807acd3c Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Tue, 28 Jun 2022 20:56:03 -0400 Subject: [PATCH 1/8] gapfillpkg edits and error corrections --- modelseedpy/fbapkg/gapfillingpkg.py | 245 +++++++++++----------------- 1 file changed, 96 insertions(+), 149 deletions(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index ba1cda6b..0a67ad8d 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -3,15 +3,15 @@ from __future__ import absolute_import import logging +logger = logging.getLogger(__name__) + 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 +from cobra import Model, Reaction, Metabolite # !!! Model is never used from modelseedpy.fbapkg.basefbapkg import BaseFBAPkg from modelseedpy.core.fbahelper import FBAHelper -logger = logging.getLogger(__name__) - default_blacklist = ["rxn12985", "rxn00238", "rxn07058", "rxn05305", "rxn00154", "rxn09037", "rxn10643", "rxn11317", "rxn05254", "rxn05257", "rxn05258", "rxn05259", "rxn05264", "rxn05268", "rxn05269", "rxn05270", "rxn05271", "rxn05272", "rxn05273", "rxn05274", "rxn05275", @@ -65,9 +65,8 @@ class GapfillingPkg(BaseFBAPkg): """ - + """ - def __init__(self, model): BaseFBAPkg.__init__(self, model, "gapfilling", {}, {}) self.gapfilling_penalties = None @@ -82,10 +81,7 @@ def build(self, template, minimum_objective=0.01): self.build_package(parameters) def get_model_index_hash(self): - """ - Determine all indices that should be gap filled - :return: - """ + """Determine all indices that should be gap filled""" index_hash = {"none": 0} for metabolite in self.model.metabolites: if re.search('_[a-z]\d+$', metabolite.id) is not None: @@ -98,7 +94,7 @@ def get_model_index_hash(self): index_hash["none":0] # Iterating over all indecies with more than 10 intracellular compounds: return index_hash - + def build_package(self, parameters): self.validate_parameters(parameters, [], { "auto_sink": ["cpd02701", "cpd11416", "cpd15302"], @@ -125,31 +121,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] @@ -166,9 +156,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: @@ -182,22 +170,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" @@ -205,71 +188,66 @@ 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: new_penalties[cobra_reaction.id]["forward"] = self.parameters["model_penalty"] # Updating metabolites in reaction to new model - metabolites = cobra_reaction.metabolites; new_stoichiometry = {} - for metabolite in metabolites: + for met in cobra_reaction.metabolites: # Adding new coefficient: - new_stoichiometry[local_remap[metabolite.id]] = metabolites[metabolite] + new_stoichiometry[local_remap[met.id]] = cobra_reaction.metabolites[met] # Zeroing out current coefficients - if local_remap[metabolite.id] != metabolite: - new_stoichiometry[metabolite] = 0 + if local_remap[met.id] != met: + new_stoichiometry[met] = 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 for cpd in new_exchange: - drain_reaction = FBAHelper.add_drain_from_metabolite_id(self.model, cpd.id, - self.parameters["default_uptake"], - self.parameters["default_excretion"]) - if drain_reaction.id not in self.cobramodel.reactions and drain_reaction.id not in new_reactions: - new_reactions[drain_reaction.id] = drain_reaction + drain_reaction = FBAHelper.add_drain_from_metabolite_id( + self.model, cpd.id, self.parameters["default_uptake"], self.parameters["default_excretion"]) + if drain_reaction.id not in self.cobramodel.reactions and drain_reaction.id not in self.new_reactions: + self.new_reactions[drain_reaction.id] = drain_reaction - # Only run this on new demands so we don't readd for all exchanges + # Only run this on new demands so we don't read for all exchanges for cpd_id in new_demand: - drain_reaction = FBAHelper.add_drain_from_metabolite_id(self.model, cpd.id, - self.parameters["default_uptake"], - self.parameters["default_excretion"]) - if drain_reaction.id not in self.cobramodel.reactions and drain_reaction.id not in new_reactions: - new_reactions[drain_reaction.id] = drain_reaction + drain_reaction = FBAHelper.add_drain_from_metabolite_id( + self.model, cpd.id, self.parameters["default_uptake"], self.parameters["default_excretion"]) + if drain_reaction.id not in self.cobramodel.reactions and drain_reaction.id not in self.new_reactions: + self.new_reactions[drain_reaction.id] = drain_reaction # 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,33 +257,30 @@ 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 + 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.id) if compartment == "e": new_exchange.append(cobra_metabolite.id) # 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() - # Adding all metabolites to model prior to adding reactions + self.new_reactions, new_penalties = {}, {} new_exchange, new_demand = self.extend_model_with_template_metabolites(template, index) for template_reaction in template.reactions: @@ -313,20 +288,20 @@ def extend_model_with_template_for_gapfilling(self, template, index): continue cobra_reaction = self.convert_template_reaction(template_reaction, index, template, 1) # 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() @@ -339,76 +314,56 @@ def extend_model_with_template_for_gapfilling(self, template, index): "forward"] = template_reaction.base_cost + template_reaction.forward_penalty # Only run this on new exchanges so we don't read for all exchanges for cpd_id in new_exchange: - drain_reaction = FBAHelper.add_drain_from_metabolite_id(self.model, cpd_id, - self.parameters["default_uptake"], - self.parameters["default_excretion"]) - if drain_reaction is not None and drain_reaction.id not in new_reactions: + drain_reaction = FBAHelper.add_drain_from_metabolite_id( + self.model, cpd_id, self.parameters["default_uptake"], self.parameters["default_excretion"]) + if drain_reaction and drain_reaction.id not in self.new_reactions: new_penalties[drain_reaction.id] = { 'added': 1, 'reverse': 1, 'forward': 1 - } - new_reactions[drain_reaction.id] = drain_reaction - # Only run this on new demands so we don't readd for all exchanges + } + self.new_reactions[drain_reaction.id] = drain_reaction + # Only run this on new demands so we don't read for all exchanges for cpd_id in new_demand: - drain_reaction = FBAHelper.add_drain_from_metabolite_id(self.model, cpd_id, - self.parameters["default_uptake"], - self.parameters["default_excretion"], - "DM_") - - if drain_reaction is not None and drain_reaction.id not in new_reactions: + drain_reaction = FBAHelper.add_drain_from_metabolite_id( + self.model, cpd_id, self.parameters["default_uptake"], self.parameters["default_excretion"], "DM_") + if drain_reaction and drain_reaction.id not in self.new_reactions: new_penalties[drain_reaction.id] = { 'added': 1, 'reverse': 1, 'forward': 1 } - new_reactions[drain_reaction.id] = drain_reaction + self.new_reactions[drain_reaction.id] = drain_reaction # 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) - + 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: 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 @@ -421,46 +376,39 @@ 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 = FBAHelper.compute_flux_values_from_variables(self.model) - 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() + new_solution = self.compute_gapfilled_solution() # !!! should flux_values be added as an argument here? return new_solution #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(flux_values) if flux_values == None: flux_values = FBAHelper.compute_flux_values_from_variables(self.model) for rxnobj in self.model.reactions: @@ -471,16 +419,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 @@ -503,7 +450,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: From b71ad9b0855851b7bb0747ca068885f08c6db8ec Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Tue, 28 Jun 2022 21:05:06 -0400 Subject: [PATCH 2/8] gapfillpkg polishing --- modelseedpy/fbapkg/gapfillingpkg.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index 0a67ad8d..3010a031 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -7,7 +7,7 @@ import re import json # !!! import is never used -from optlang.symbolics import Zero +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 @@ -65,7 +65,7 @@ class GapfillingPkg(BaseFBAPkg): """ - + """ def __init__(self, model): BaseFBAPkg.__init__(self, model, "gapfilling", {}, {}) @@ -81,7 +81,10 @@ def build(self, template, minimum_objective=0.01): self.build_package(parameters) def get_model_index_hash(self): - """Determine all indices that should be gap filled""" + """ + Determine all indices that should be gap filled + :return: + """ index_hash = {"none": 0} for metabolite in self.model.metabolites: if re.search('_[a-z]\d+$', metabolite.id) is not None: @@ -94,7 +97,7 @@ def get_model_index_hash(self): index_hash["none":0] # Iterating over all indecies with more than 10 intracellular compounds: return index_hash - + def build_package(self, parameters): self.validate_parameters(parameters, [], { "auto_sink": ["cpd02701", "cpd11416", "cpd15302"], @@ -203,7 +206,7 @@ def extend_model_with_model_for_gapfilling(self, source_model, index): m = re.search('(.+)_([a-z])\d+$', modelreaction.id) if m[1] not in self.parameters["blacklist"]: cobra_reaction = modelreaction.copy() - cobra_reaction.id = m[1] + "_" + m[2] + index + 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] = {} @@ -450,7 +453,7 @@ def run_test_conditions(self, condition_list, solution = None, max_iterations = return solution def filter_database_based_on_tests(self,test_conditions): - filtered_list = [] + filtered_list = [] with self.model: rxnlist = [] for reaction in self.model.reactions: From d75426cb0771b2bc4a3326e2c492d95328cb4dd8 Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Tue, 28 Jun 2022 21:42:23 -0400 Subject: [PATCH 3/8] gapfillpkg polishing --- modelseedpy/fbapkg/gapfillingpkg.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index 3010a031..18137ca3 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -216,13 +216,14 @@ def extend_model_with_model_for_gapfilling(self, source_model, index): if cobra_reaction.upper_bound > 0: new_penalties[cobra_reaction.id]["forward"] = self.parameters["model_penalty"] # Updating metabolites in reaction to new model + metabolites = cobra_reaction.metabolites; new_stoichiometry = {} - for met in cobra_reaction.metabolites: + for metabolite in metabolites: # Adding new coefficient: - new_stoichiometry[local_remap[met.id]] = cobra_reaction.metabolites[met] + new_stoichiometry[local_remap[metabolite.id]] = metabolites[metabolite] # Zeroing out current coefficients - if local_remap[met.id] != met: - new_stoichiometry[met] = 0 + 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: self.model.reactions.get_by_id(cobra_reaction.id).lower_bound = cobra_reaction.lower_bound From d1960a2d8eb11bf92dbdc1f894ff69da166a8092 Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Tue, 28 Jun 2022 21:46:32 -0400 Subject: [PATCH 4/8] gapfillpkg polishing --- modelseedpy/fbapkg/gapfillingpkg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index 18137ca3..f7918e69 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -406,13 +406,13 @@ def binary_check_gapfilling_solution(self, solution=None, flux_values=None): 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() # !!! should flux_values be added as an argument here? + new_solution = self.compute_gapfilled_solution() return new_solution #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(flux_values) + solution = self.compute_gapfilled_solution() # !!! should flux_values be added as an argument here? if flux_values == None: flux_values = FBAHelper.compute_flux_values_from_variables(self.model) for rxnobj in self.model.reactions: From a67148fdca2bd5056da00c4a1b715416b9c95a23 Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Mon, 4 Jul 2022 10:26:21 -0400 Subject: [PATCH 5/8] 1->True conversion in alignment with the merged edits from template.py --- modelseedpy/fbapkg/gapfillingpkg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index f7918e69..b9cd27ce 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -290,7 +290,7 @@ def extend_model_with_template_for_gapfilling(self, 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 self.new_reactions: # Adding any template reactions missing from the present model From a5e0de24d4656dc01a661a9824cf67475ffcd250 Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Mon, 4 Jul 2022 10:34:37 -0400 Subject: [PATCH 6/8] comments for possible integration with template.py --- modelseedpy/fbapkg/gapfillingpkg.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index b9cd27ce..8973c52e 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -11,6 +11,7 @@ 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? default_blacklist = ["rxn12985", "rxn00238", "rxn07058", "rxn05305", "rxn00154", "rxn09037", "rxn10643", "rxn11317", "rxn05254", "rxn05257", "rxn05258", "rxn05259", "rxn05264", "rxn05268", @@ -70,6 +71,7 @@ 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 = { @@ -266,7 +268,7 @@ def extend_model_with_template_metabolites(self, template, index='0'): 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 + 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]) @@ -355,13 +357,13 @@ def convert_template_compound(self, template_compound, index, template): met.annotation["seed.compound"] = base_id return met - def convert_template_reaction(self, template_reaction, index, template, for_gapfilling=1): + 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 From f249bf9acec40d2010c519b2a1912d0213f4850f Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Fri, 22 Jul 2022 22:45:07 -0400 Subject: [PATCH 7/8] gapfillingpkg polishing --- modelseedpy/fbapkg/gapfillingpkg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index 8973c52e..db66c25b 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -285,8 +285,8 @@ def extend_model_with_template_metabolites(self, template, index='0'): # 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): - # Adding all metabolites to model prior to adding reactions 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: From a277496ce439c2346dcecebe053b111f3d9aae0c Mon Sep 17 00:00:00 2001 From: Andrew Freiburger Date: Fri, 22 Jul 2022 22:51:58 -0400 Subject: [PATCH 8/8] gapfillingpkg additional edits --- modelseedpy/fbapkg/gapfillingpkg.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/modelseedpy/fbapkg/gapfillingpkg.py b/modelseedpy/fbapkg/gapfillingpkg.py index db66c25b..90a1210d 100755 --- a/modelseedpy/fbapkg/gapfillingpkg.py +++ b/modelseedpy/fbapkg/gapfillingpkg.py @@ -3,7 +3,6 @@ from __future__ import absolute_import import logging -logger = logging.getLogger(__name__) import re import json # !!! import is never used @@ -13,6 +12,8 @@ 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__) + default_blacklist = ["rxn12985", "rxn00238", "rxn07058", "rxn05305", "rxn00154", "rxn09037", "rxn10643", "rxn11317", "rxn05254", "rxn05257", "rxn05258", "rxn05259", "rxn05264", "rxn05268", "rxn05269", "rxn05270", "rxn05271", "rxn05272", "rxn05273", "rxn05274", "rxn05275",