From 90ef52069698469bd5fdf5dfbf29d886687a9319 Mon Sep 17 00:00:00 2001 From: Andrew's WorkMacbook Date: Tue, 30 Jul 2024 16:33:39 -0500 Subject: [PATCH] debugging during CommScores development --- modelseedpy/community/commhelper.py | 28 ++++++++++++---------------- modelseedpy/community/mscommunity.py | 2 +- modelseedpy/core/msminimalmedia.py | 5 +++-- modelseedpy/core/msmodelutl.py | 4 ++++ 4 files changed, 20 insertions(+), 19 deletions(-) diff --git a/modelseedpy/community/commhelper.py b/modelseedpy/community/commhelper.py index 76626b50..49e10048 100644 --- a/modelseedpy/community/commhelper.py +++ b/modelseedpy/community/commhelper.py @@ -4,12 +4,12 @@ from modelseedpy.core.msmodelutl import MSModelUtil from modelseedpy.core.fbahelper import FBAHelper from cobra import Model, Reaction, Metabolite +from optlang import Constraint, Objective from cobra.medium import minimal_medium # from commscores import GEMCompatibility from cobra.flux_analysis import pfba from collections import OrderedDict from optlang.symbolics import Zero -from optlang import Constraint from math import inf, isclose from pandas import DataFrame from pprint import pprint @@ -60,7 +60,6 @@ def build_from_species_models(org_models, model_id=None, name=None, abundances=N # Renaming compartments output = MSModelUtil.parse_id(met) # if printing: print(met, output) - print(met.id, output) if output is None: if printing: print(f"The {met.id} ({output}; {hasattr(met, 'compartment')}) is unpredictable.") met.id = correct_nonMSID(met, (met.id, "c"), model_index) @@ -75,7 +74,6 @@ def build_from_species_models(org_models, model_id=None, name=None, abundances=N else: met.compartment = compartment + str(index) met.id = name + "_" + met.compartment - print(met.id) new_metabolites.add(met) if "cpd11416_c" in met.id or "biomass" in met.name: met.name = f"{met.id}_{model_util.model.id}" @@ -85,7 +83,7 @@ def build_from_species_models(org_models, model_id=None, name=None, abundances=N if rxn.id[0:3] != "EX_": ## biomass reactions if re.search('^(bio)(\d+)$', rxn.id) or "biomass" in rxn.id: - print(rxn.id) + print(rxn.id, "from", model_util.id, "becomes", end=" ") index = int(re.sub(r"(^biomass)", "", rxn.id)) if "biomass" in rxn.id else int(re.sub(r"(^bio)", "", rxn.id)) if biomass_index == 2: while f"bio{biomass_index}" in model_reaction_ids: biomass_index += 1 @@ -134,11 +132,6 @@ def build_from_species_models(org_models, model_id=None, name=None, abundances=N # else: # # TODO develop a method for compartmentalizing models without editing all reaction IDs or assuming their syntax # pass - # adds only unique reactions and metabolites to the community model - newmodel = Model(model_id or "+".join([model.id for model in models]), - name or "+".join([model.name for model in models])) - newmodel.add_reactions(FBAHelper.filter_cobra_set(new_reactions)) - newmodel.add_metabolites(FBAHelper.filter_cobra_set(new_metabolites)) # Create community biomass comm_biomass = Metabolite("cpd11416_c0", None, "Community biomass", 0, "c0") @@ -150,23 +143,26 @@ def build_from_species_models(org_models, model_id=None, name=None, abundances=N else: abundances = {met: -abundances[memberID]["abundance"] for memberID, met in member_biomasses.items()} else: abundances = {met: -1 / len(member_biomasses) for met in member_biomasses.values()} - # TODO - add the biomass reactions instead of thje biomass metabolites for the commKinetics - print(abundances) + ## TODO - add the biomass reactions instead of the biomass metabolites for the commKinetics ## define community biomass components metabolites.update(abundances) comm_biorxn = Reaction(id="bio1", name="bio1", lower_bound=0, upper_bound=1000) comm_biorxn.add_metabolites(metabolites) - print(comm_biorxn) + + # adds only unique reactions and metabolites to the community model + newmodel = Model(model_id or "+".join([model.id for model in models]), + name or "+".join([model.name for model in models])) + newmodel.add_reactions(FBAHelper.filter_cobra_set(new_reactions)) + newmodel.add_metabolites(FBAHelper.filter_cobra_set(new_metabolites)) newmodel.add_reactions([comm_biorxn]) - # update model components + newmodel.objective = Objective(comm_biorxn.flux_expression) newutl = MSModelUtil(newmodel) - newutl.add_objective(comm_biorxn.flux_expression) + # newutl.add_objective(comm_biorxn.flux_expression) newutl.model.add_boundary(comm_biomass, "sink") # Is a sink reaction for reversible cpd11416_c0 consumption necessary? ## proportionally limit the fluxes to their abundances - # print(abundances) # add the metadata of community composition - print(newutl.model.reactions.bio1.reaction) + print("Community objective", newutl.model.objective.expression) if hasattr(newutl.model, "_context"): newutl.model._contents.append(member_biomasses) elif hasattr(newutl.model, "notes"): newutl.model.notes.update(member_biomasses) # print([cons.name for cons in newutl.model.constraints]) diff --git a/modelseedpy/community/mscommunity.py b/modelseedpy/community/mscommunity.py index dc955266..a73ac549 100644 --- a/modelseedpy/community/mscommunity.py +++ b/modelseedpy/community/mscommunity.py @@ -44,7 +44,7 @@ def __init__(self, community, biomass_cpd, ID=None, index=None, abundance=0): # print(rxn.id, rxn.reaction, "\t\t\t", end="\r") rxnComp = FBAHelper.rxn_compartment(rxn) if rxnComp is None: print(f"The reaction {rxn.id} compartment is undefined.") - if rxnComp[1:] == '': print(rxn, rxnComp) + if rxnComp[1:] == '': print("no compartment", rxn, rxnComp) elif int(rxnComp[1:]) == self.index and 'bio' not in rxn.name: self.reactions.append(rxn) if self.biomass_cpd in rxn.metabolites: # print(rxn.reaction) diff --git a/modelseedpy/core/msminimalmedia.py b/modelseedpy/core/msminimalmedia.py index afa54004..1d235b21 100644 --- a/modelseedpy/core/msminimalmedia.py +++ b/modelseedpy/core/msminimalmedia.py @@ -82,7 +82,7 @@ def _influx_objective(model_util, interacting): return influxes @staticmethod - def minimize_flux(org_model, min_growth=None, environment=None, interacting=True, printing=True): + def minimize_flux(org_model, min_growth=None, environment=None, interacting=True, pfba=True, printing=True): """minimize the total in-flux of exchange reactions in the model""" if org_model.slim_optimize() == 0: raise ObjectiveError(f"The model {org_model.id} possesses an objective value of 0 in complete media, " @@ -90,7 +90,8 @@ def minimize_flux(org_model, min_growth=None, environment=None, interacting=True model_util = MSModelUtil(org_model, True) model_util.add_medium(environment or model_util.model.medium) # define the MILP - min_growth = model_util.model.slim_optimize() if min_growth is None else min(min_growth, model_util.model.slim_optimize()) + sol_growth = model_util.run_fba(None, pfba).fluxes[model_util.biomass_objective] + min_growth = sol_growth if min_growth is None else min(min_growth, sol_growth) # min_flux = MSMinimalMedia._min_consumption_objective(model_util, interacting) media_exchanges = MSMinimalMedia._influx_objective(model_util, interacting) # parse the minimal media diff --git a/modelseedpy/core/msmodelutl.py b/modelseedpy/core/msmodelutl.py index 263fdd0d..5f875f43 100644 --- a/modelseedpy/core/msmodelutl.py +++ b/modelseedpy/core/msmodelutl.py @@ -124,6 +124,10 @@ def __init__(self, model, copy=False, environment=None): self.test_objective = None self.reaction_scores = None self.score = None + try: + self.biomass_objective = list(self.model.objective.variables)[0].name + except IndexError: + print(f"The {self.id} has an improperly defined objective function") self.integrated_gapfillings = [] self.attributes = {} if hasattr(self.model, "computed_attributes"):