debugging during CommScores development

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])

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)

modelseedpy/core/msminimalmedia.py

@@ -82,15 +82,16 @@ 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, "
                                  "which is incompatible with minimal media computations.")
         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

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"):