core: add the API surface MSCommunity needs from ModelSEEDpy

modelseedpy/core/__init__.py

@@ -3,7 +3,11 @@
 
 from modelseedpy.core.rast_client import RastClient
 from modelseedpy.core.msgenome import MSGenome
-from modelseedpy.core.fbahelper import FBAHelper
+from modelseedpy.core.fbahelper import (
+    FBAHelper,
+    bioFlux_check,
+    minimizeFlux_withGrowth,
+)
 from modelseedpy.core.msbuilder import MSBuilder
 from modelseedpy.core.msmedia import MSMedia
 from modelseedpy.core.mseditorapi import MSEditorAPI, MSEquation

modelseedpy/core/exceptions.py

@@ -1,5 +1,6 @@
 # -*- coding: utf-8 -*-
 
+
 # Adding a few exception classes to handle different types of errors in a central file
 class ModelSEEDError(Exception):
     """Error in ModelSEED execution logic"""
@@ -24,3 +25,33 @@ class GapfillingError(Exception):
     """Error in model gapfilling"""
 
     pass
+
+
+class ParameterError(Exception):
+    """Error in a parameterization"""
+
+    pass
+
+
+class ObjectAlreadyDefinedError(Exception):
+    """Error from defining an object that is already defined"""
+
+    pass
+
+
+class NoFluxError(Exception):
+    """Error for FBA solutions that carry no flux"""
+
+    pass
+
+
+class ObjectiveError(Exception):
+    """Erroneous assignment of a secondary objective via a constraint"""
+
+    pass
+
+
+class ModelError(Exception):
+    """Errors in a model that corrupt the simulation"""
+
+    pass

modelseedpy/core/fbahelper.py

@@ -16,6 +16,7 @@
 from scipy.odr.odrpack import Output  # !!! Output is never used
 from chemw import ChemMW
 from warnings import warn
+from modelseedpy.core.exceptions import ObjectiveError
 
 # from Carbon.Aliases import false
 
@@ -353,3 +354,76 @@ def parse_df(df):
         return array(
             dtype=object, object=[array(df.index), array(df.columns), df.to_numpy()]
         )
+
+    # --- helpers required by MSCommunity -------------------------------------
+    @staticmethod
+    def isnumber(string):
+        if str(string) in ["nan", "inf", "True", "False"]:
+            return False
+        try:
+            float(string)
+            return True
+        except:
+            return False
+
+    @staticmethod
+    def mediaName(media):
+        if media == None:
+            return "Complete"
+        return media.id
+
+    @staticmethod
+    def rxn_mets_list(rxn):
+        return [met for met in rxn.reactants + rxn.products]
+
+    @staticmethod
+    def solution_to_variables_dict(solution, model):
+        return {model.variables.get(key): flux for key, flux in solution.fluxes.items()}
+
+    @staticmethod
+    def convert_kbase_media(kbase_media, uniform_uptake=1000):
+        if uniform_uptake is None:
+            return {
+                "EX_" + exID: -bound[0]
+                for exID, bound in kbase_media.get_media_constraints().items()
+            }
+        return {
+            "EX_" + exID: uniform_uptake
+            for exID in kbase_media.get_media_constraints().keys()
+        }
+
+
+# --- FBA solution helpers required by MSCommunity ---------------------------
+
+
+def bioFlux_check(model, sol=None, sol_dict=None, min_growth=0.1):
+    """Verify that a solution maintains the requested minimal biomass flux.
+
+    :param model: the cobra model that produced the solution
+    :param sol: a cobra Solution, used when ``sol_dict`` is not supplied
+    :param sol_dict: a {variable: flux} mapping, computed from ``sol`` if absent
+    :param min_growth: the biomass flux the solution was required to maintain
+    :raises ObjectiveError: when the simulated growth falls below ``min_growth``
+    :return: the {variable: flux} mapping
+    """
+    sol_dict = sol_dict or FBAHelper.solution_to_variables_dict(sol, model)
+    simulated_growth = sum(
+        [flux for var, flux in sol_dict.items() if re.search(r"(^bio\d+$)", var.name)]
+    )
+    if simulated_growth < min_growth * 0.9999 and simulated_growth + min_growth > 1e-8:
+        raise ObjectiveError(
+            f"The assigned minimal_growth of {min_growth} was not maintained during "
+            f"the simulation, where the observed growth value was {simulated_growth}."
+        )
+    if sol is not None and sol.status != "optimal":
+        logger.warning(f"The solution is {sol.status}, not optimal.")
+    return sol_dict
+
+
+def minimizeFlux_withGrowth(model_util, min_growth, obj):
+    """Minimize ``obj`` subject to maintaining at least ``min_growth`` biomass flux."""
+    model_util.add_minimal_objective_cons(min_growth, name="min_growth")
+    model_util.add_objective(obj, "min")
+    sol = model_util.model.optimize()
+    sol_dict = bioFlux_check(model_util.model, sol, None, min_growth)
+    return sol, sol_dict

modelseedpy/core/msmodelutl.py

@@ -1,5 +1,6 @@
 # -*- coding: utf-8 -*-
 import logging
+from math import isclose
 import re
 import time
 import json
@@ -8,6 +9,7 @@
 from modelseedpy.fbapkg.mspackagemanager import MSPackageManager
 from modelseedpy.biochem.modelseed_biochem import ModelSEEDBiochem
 from modelseedpy.core.fbahelper import FBAHelper
+from modelseedpy.core.exceptions import ModelError
 
 logger = logging.getLogger(__name__)
 logger.setLevel(logging.DEBUG)
@@ -93,9 +95,36 @@ def get(model, create_if_missing=True):
         else:
             return None
 
-    def __init__(self, model):
+    def __init__(self, model, copy=False, environment=None, climit=None, o2limit=None):
+        """
+        :param model: the cobra model to wrap
+        :param copy: work on an objective-preserving copy rather than the model itself
+        :param environment: a media dict applied to the model before anything else
+        :param climit: total carbon uptake cap, in mmol/gDW/h; False disables
+        :param o2limit: oxygen uptake cap, in mmol/gDW/h; False disables
+
+        The four keyword arguments all default to the previous no-op behaviour,
+        so existing callers of ``MSModelUtil(model)`` are unaffected.
+        """
         self.model = model
-        self.pkgmgr = MSPackageManager.get_pkg_mgr(model)
+        if environment is not None:
+            self.add_medium(environment)
+        self.id = model.id
+        if copy:
+            org_obj_val = model.slim_optimize()
+            self.model = model.copy()
+            self.model.objective = model.objective
+            new_obj_val = self.model.slim_optimize()
+            if (
+                not isclose(org_obj_val, new_obj_val, rel_tol=1e-2)
+                and org_obj_val > 1e-2
+            ):
+                raise ModelError(
+                    f"The {model.id} objective value is corrupted by being copied, "
+                    f"where the original objective value is {org_obj_val} and the "
+                    f"new objective value is {new_obj_val}."
+                )
+        self.pkgmgr = MSPackageManager.get_pkg_mgr(self.model)
         self.atputl = None
         self.gfutl = None
         self.metabolite_hash = None
@@ -104,6 +133,28 @@ def __init__(self, model):
         self.reaction_scores = None
         self.score = None
         self.integrated_gapfillings = []
+        self.apply_uptake_limits(climit, o2limit)
+
+    def apply_uptake_limits(self, climit=None, o2limit=None):
+        """Cap total carbon and oxygen uptake.
+
+        Passing ``False`` for both is an explicit opt-out and returns immediately;
+        passing ``None`` for both leaves the model untouched, which is the
+        behaviour callers that predate these arguments rely on.
+        """
+        if climit is False and o2limit is False:
+            return
+        if climit is None and o2limit is None:
+            return
+        if not FBAHelper.isnumber(climit) and FBAHelper.isnumber(o2limit):
+            climit = 3 * o2limit
+        elif not FBAHelper.isnumber(climit):
+            climit = 60
+        if not FBAHelper.isnumber(o2limit):
+            o2limit = climit / 3
+        self.pkgmgr.getpkg("ElementUptakePkg").build_package({"C": climit})
+        if "EX_cpd00007_e0" in [rxn.id for rxn in self.model.reactions]:
+            self.model.reactions.get_by_id("EX_cpd00007_e0").lower_bound = -o2limit
 
     def compute_automated_reaction_scores(self):
         """
@@ -938,3 +989,113 @@ def parse_id(object):
             m = re.search("(.+)_([a-z]+)(\d*)$", object.id)
             return (m[1], m[2], m[3])
         return None
+
+    # --- helpers required by MSCommunity -------------------------------------
+    def add_medium(self, media, uniform_uptake=None):
+        # add the new media and its flux constraints
+        exIDs = [exRXN.id for exRXN in self.exchange_list()]
+        if not hasattr(media, "items"):
+            media = FBAHelper.convert_kbase_media(media)
+        elif not any(["EX_" in x for x in list(media.keys())]):
+            media = {"EX_" + k + "_e0": v for k, v in media.items()}
+        self.model.medium = {ex: uptake for ex, uptake in media.items() if ex in exIDs}
+        if uniform_uptake is not None:
+            self.model.medium = dict(
+                zip(
+                    list(self.model.medium.keys()),
+                    [uniform_uptake] * len(self.model.medium),
+                )
+            )
+        return self.model.medium
+
+    def add_minimal_objective_cons(
+        self, min_value=0.1, objective_expr=None, name="min_value"
+    ):
+        if name not in self.model.constraints:
+            objective_expr = objective_expr or self.model.objective.expression
+            self.create_constraint(
+                self.model.problem.Constraint(
+                    objective_expr, lb=min_value, ub=None, name=name
+                )
+            )
+            # print(self.model.constraints["min_value"])
+        else:
+            print(
+                f"The {name} constraint already exists in {self.model.id}, "
+                f"hence the lb is simply updated from"
+                f" {self.model.constraints[name].lb} to {min_value}.\n"
+            )
+            self.model.constraints[name].lb = min_value
+
+    def add_objective(self, objective, direction="max", coef=None):
+        self.model.objective = self.model.problem.Objective(
+            objective, direction=direction
+        )
+        self.model.solver.update()
+        if coef:
+            self.model.objective.set_linear_coefficients(coef)
+            self.model.solver.update()
+
+    def carbon_exchange_list(self, include_unknown=True):
+        if not include_unknown:
+            return [
+                ex for ex in self.exchange_list() if "C" in ex.reactants[0].elements
+            ]
+        return [
+            ex
+            for ex in self.exchange_list()
+            if not ex.reactants[0].elements or "C" in ex.reactants[0].elements
+        ]
+
+    def carbon_exchange_mets_list(self, include_unknown=True):
+        return self.metabolites_set(self.carbon_exchange_list(include_unknown))
+
+    def create_constraint(self, constraint, coef=None, sloppy=False, printing=False):
+        # if printing:   print(coef)
+        self.model.add_cons_vars(constraint, sloppy=sloppy)
+        self.model.solver.update()
+        if coef:
+            constraint.set_linear_coefficients(coef)
+        self.model.solver.update()
+
+    def exchange_mets_list(self):
+        return self.metabolites_set(self.exchange_list())
+
+    def remove_constraint(self, consName):
+        for cons in self.model.constraints:
+            if consName not in cons.name:
+                continue
+            # if self.printing:   print(f"Removing {consName} from {self.model.id}")
+            self.model.remove_cons_vars(cons)
+
+    def run_fba(self, media=None, pfba=False, fva_reactions=None):
+        from cobra import flux_analysis
+
+        if media:
+            self.pkgmgr.getpkg("KBaseMediaPkg").build_package(media)
+        if pfba:
+            return flux_analysis.pfba(self.model)
+        if fva_reactions is not None:
+            return flux_analysis.variability.flux_variability_analysis(
+                self.model, fva_reactions
+            )
+        return self.model.optimize()
+
+    def standard_exchanges(self):
+        for ex in self.exchange_list():
+            if len(ex.reactants) != 1 and len(ex.products) != 0:
+                raise ModelError(
+                    f"The ex {ex.id} possesses {len(ex.reactants)} reactants and "
+                    f"{len(ex.products)} products, which are non-standard and are incompatible"
+                    f" with various ModelSEED operations."
+                )
+
+    def metabolites_set(self, reactions_set=None, ids=False):
+        rxns = reactions_set or self.model.reactions
+        if ids:
+            return {met.id for rxn in rxns for met in rxn.metabolites}
+        return {met for rxn in rxns for met in rxn.metabolites}
+
+    def remove_cons_vars(self, vars_cons):
+        self.model.remove_cons_vars(vars_cons)
+        self.model.solver.update()