From 9ff8b0fca3a90895c1572c3c6d50f8f2541f4e3e Mon Sep 17 00:00:00 2001 From: yseif Date: Mon, 8 Jul 2019 17:17:29 -0700 Subject: [PATCH 1/2] Gapfill algorithm adapted to the multi-strain workflow --- gapfill_multi.py | 65 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 gapfill_multi.py diff --git a/gapfill_multi.py b/gapfill_multi.py new file mode 100644 index 000000000..4305da677 --- /dev/null +++ b/gapfill_multi.py @@ -0,0 +1,65 @@ +def gapfill_multi(model, missing_genes, **kwargs): + + """ + + Generate a list of gapfilling reactions from a list of missing genes for a strain-specific model. + + :param model: COBRA model for the base strain with the objective coefficient for the reaction of interest (e.g. biomass reaction) set to 1. + + :param missing_genes: list of genes with no homologs in the strain of interest. + + :param lower_bound: minimum allowable yield of gapfilled model. + + :param biomass: override the current model settings and temporarily assign the objective coefficient for a function of interest to 1. + + :return: a list of gapfilling reactions. + + """ + + # model.solver = 'gurobi' + + if 'lower_bound' in kwargs.keys(): + lower_bound = kwargs['lower_bound'] + else: + lower_bound = model.optimize().objective_value*0.5 + + biomass_reactions = [rx.id for rx in model.reactions if rx.objective_coefficient == 1] + if 'biomass' in kwargs.keys(): + biomass = kwargs['biomass'] + if len(biomass_reactions) > 1: + for rx in set(biomass_reactions) - {biomass}: + model.reactions.get_by_id(rx).objective_coefficient = 0 + + else: + if len(biomass_reactions) > 1: + raise Exception("This model has more than one objective. \n Please adjust the objective coefficient to 1 for the chosen objective reaction (e.g. biomass or ATP) and 0 for the rest of the reactions, \n or specify the reaction ID to use as an objective.") + if len(biomass_reactions) > 1: + raise Exception("The model doesn't have an objective function. Please set the appropriate objective coefficient to 1, or specify the reaction ID to use as an objective.") + biomass = biomass_reactions[0] + + + model.solver.configuration.tolerances.feasibility = 1e-9 + constraints = [] + indicators = [] + + for rx in cobra.manipulation.find_gene_knockout_reactions(model, missing_genes): + + indicator = model.problem.Variable('%s_i'%rx.id , type = 'binary') + indicators.append(indicator) + + new_cstr1 = model.problem.Constraint( rx.flux_expression - rx.upper_bound*indicator ,ub = 0) + new_cstr2 = model.problem.Constraint(-rx.flux_expression + rx.lower_bound*indicator ,ub = 0) + constraints += [new_cstr1, new_cstr2] + model.add_cons_vars([new_cstr1, new_cstr2, indicator]) + + model.reactions.get_by_id(biomass).lower_bound = lower_bound + model.objective = model.problem.Objective(-sum(indicators)) + sol = model.optimize() + indicator_results = [ind.name[:-2] for ind in indicators if ind.primal != 0.0] + + # removing changes to model + model.remove_cons_vars(constraints+indicators) + for rx in set(biomass_reactions): + model.reactions.get_by_id(rx).objective_coefficient = 1 + + return indicator_results \ No newline at end of file From 12336e57c42005ce4206ffb807643f92c0b38042 Mon Sep 17 00:00:00 2001 From: yseif Date: Mon, 8 Jul 2019 17:29:49 -0700 Subject: [PATCH 2/2] gapfilling adapted to multi-strain workflow --- cobra/flux_analysis/gapfill_multi.py | 65 ++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 cobra/flux_analysis/gapfill_multi.py diff --git a/cobra/flux_analysis/gapfill_multi.py b/cobra/flux_analysis/gapfill_multi.py new file mode 100644 index 000000000..4305da677 --- /dev/null +++ b/cobra/flux_analysis/gapfill_multi.py @@ -0,0 +1,65 @@ +def gapfill_multi(model, missing_genes, **kwargs): + + """ + + Generate a list of gapfilling reactions from a list of missing genes for a strain-specific model. + + :param model: COBRA model for the base strain with the objective coefficient for the reaction of interest (e.g. biomass reaction) set to 1. + + :param missing_genes: list of genes with no homologs in the strain of interest. + + :param lower_bound: minimum allowable yield of gapfilled model. + + :param biomass: override the current model settings and temporarily assign the objective coefficient for a function of interest to 1. + + :return: a list of gapfilling reactions. + + """ + + # model.solver = 'gurobi' + + if 'lower_bound' in kwargs.keys(): + lower_bound = kwargs['lower_bound'] + else: + lower_bound = model.optimize().objective_value*0.5 + + biomass_reactions = [rx.id for rx in model.reactions if rx.objective_coefficient == 1] + if 'biomass' in kwargs.keys(): + biomass = kwargs['biomass'] + if len(biomass_reactions) > 1: + for rx in set(biomass_reactions) - {biomass}: + model.reactions.get_by_id(rx).objective_coefficient = 0 + + else: + if len(biomass_reactions) > 1: + raise Exception("This model has more than one objective. \n Please adjust the objective coefficient to 1 for the chosen objective reaction (e.g. biomass or ATP) and 0 for the rest of the reactions, \n or specify the reaction ID to use as an objective.") + if len(biomass_reactions) > 1: + raise Exception("The model doesn't have an objective function. Please set the appropriate objective coefficient to 1, or specify the reaction ID to use as an objective.") + biomass = biomass_reactions[0] + + + model.solver.configuration.tolerances.feasibility = 1e-9 + constraints = [] + indicators = [] + + for rx in cobra.manipulation.find_gene_knockout_reactions(model, missing_genes): + + indicator = model.problem.Variable('%s_i'%rx.id , type = 'binary') + indicators.append(indicator) + + new_cstr1 = model.problem.Constraint( rx.flux_expression - rx.upper_bound*indicator ,ub = 0) + new_cstr2 = model.problem.Constraint(-rx.flux_expression + rx.lower_bound*indicator ,ub = 0) + constraints += [new_cstr1, new_cstr2] + model.add_cons_vars([new_cstr1, new_cstr2, indicator]) + + model.reactions.get_by_id(biomass).lower_bound = lower_bound + model.objective = model.problem.Objective(-sum(indicators)) + sol = model.optimize() + indicator_results = [ind.name[:-2] for ind in indicators if ind.primal != 0.0] + + # removing changes to model + model.remove_cons_vars(constraints+indicators) + for rx in set(biomass_reactions): + model.reactions.get_by_id(rx).objective_coefficient = 1 + + return indicator_results \ No newline at end of file