ENH: allow different annotation keys

pytfa/core/model.py

@@ -121,7 +121,7 @@ def add_variable(self, kind, hook, queue=False, **kwargs):
                    **kwargs)
 
         self._var_dict[var.name] = var
-        self.logger.debug('Added variable: {}'.format(var.name))
+        # self.logger.debug('Added variable: {}'.format(var.name))
         # self.add_cons_vars(var.variable)
 
         return var
@@ -151,7 +151,7 @@ def add_constraint(self, kind, hook, expr, queue=False,**kwargs):
                     queue=queue,
                     **kwargs)
         self._cons_dict[cons.name] = cons
-        self.logger.debug('Added constraint: {}'.format(cons.name))
+        # self.logger.debug('Added constraint: {}'.format(cons.name))
         # self.add_cons_vars(cons.constraint)
 
         return cons

pytfa/redgem/lumpgem.py

@@ -103,7 +103,7 @@ def __init__(self, tfa_model, additional_core_reactions, params):
             elif is_exchange(rxn):
                 self._exchanges.append(rxn)
             # If it is a transport reaction
-            elif check_transport_reaction(rxn):
+            elif check_transport_reaction(rxn, tfa_model.annotation_key):
                 self._transports.append(rxn)
             # If it's a core reaction
             elif rxn.subsystem in self.core_subsystems:

pytfa/thermo/metabolite.py

@@ -16,8 +16,7 @@
 
 from . import std
 from ..utils.numerics import BIGM_THERMO
-
-CPD_PROTON = 'cpd00067'
+from .utils import PROTON
 
 DEFAULT_VAL = BIGM_THERMO
 
@@ -212,7 +211,7 @@ def calcDGis(self):
             print("Computing DGis...")
 
         # Special case for protons...
-        if self.id == CPD_PROTON:
+        if self.id in PROTON:
             if self.debug:
                 print("Found proton")
             return -self.RT * log(10 ** -self.pH)
@@ -363,7 +362,7 @@ def calcDGspA(self):
             print('Computing DGspA()...')
 
         # Case of the proton
-        if self.id == CPD_PROTON:
+        if self.id in PROTON:
             if self.debug:
                 print('Proton found, returning standard values')
             # we do not adjust for proton so just return the values
@@ -441,4 +440,4 @@ def calcDGspA(self):
         if self.debug:
             print("Found deltaGspA : " + str(deltaGspA))
 
-        return (deltaGspA, sp_charge, sp_nH)
+        return (deltaGspA, sp_charge, sp_nH)

pytfa/thermo/reaction.py

@@ -14,16 +14,16 @@
 from math import log, sqrt
 
 from . import std
-from .utils import find_transported_mets
-from .metabolite import CPD_PROTON
+from .utils import find_transported_mets, PROTON, WATER
 
 ###################
 # REACTIONS TOOLS #
 ###################
 
 
 
-def calcDGtpt_rhs(reaction, compartmentsData, thermo_units):
+def calcDGtpt_rhs(
+        reaction, compartmentsData, thermo_units, annotation_key="seed_id"):
     """ Calculates the RHS of the deltaG constraint, i.e. the sum of the
     non-concentration terms
 
@@ -69,59 +69,59 @@ def calcDGtpt_rhs(reaction, compartmentsData, thermo_units):
     RT_sum_H_LC_tpt = 0  # to include the differential proton concentration
     # effects if protons are transported
 
-    transportedMets = find_transported_mets(reaction)
+    transportedMets = find_transported_mets(reaction, annotation_key)
     compartments = {'reactant': [], 'product': []}
 
-    for seed_id in transportedMets:
+    for met_id in transportedMets:
         for metType in ['reactant', 'product']:
-            if seed_id != 'cpd00001':
-                met = transportedMets[seed_id][metType]
+            if met_id not in WATER:
+                met = transportedMets[met_id][metType]
                 pH_comp = met.thermo.pH
                 ionicStr_comp = met.thermo.ionicStr
 
                 deltaGfsp = met.thermo.deltaGf_tr
 
                 compartments[metType].append(met.compartment)
                 sum_stoich_NH += ((1 if metType == 'product' else -1)
-                                  * transportedMets[seed_id]['coeff']
+                                  * transportedMets[met_id]['coeff']
                                   * met.thermo.nH_std
                                   * RT
                                   * log(10 ** -pH_comp))
                 sum_deltaGFis_trans += ((1 if metType == 'product' else -1)
-                                        * transportedMets[seed_id]['coeff']
+                                        * transportedMets[met_id]['coeff']
                                         * deltaGfsp)
             else:
                 compartments[metType].append('')
 
-            if seed_id == CPD_PROTON:
-                met = transportedMets[seed_id][metType]
+            if met_id in PROTON:
+                met = transportedMets[met_id][metType]
                 pH_comp = met.thermo.pH
                 RT_sum_H_LC_tpt += ((1 if metType == 'product' else -1)
                                     * RT
-                                    * transportedMets[seed_id]['coeff']
+                                    * transportedMets[met_id]['coeff']
                                     * log(10 ** -pH_comp))
 
     # calculate the transport of any ions
     # membrane potential is always defined as inside - outside
     # we should take the larger stoich of the transported compound
     sum_F_memP_charge = 0
 
-    for seed_id in transportedMets:
-        if seed_id != 'cpd00001':
-            out_comp = transportedMets[seed_id]['reactant'].compartment
-            in_comp = transportedMets[seed_id]['product'].compartment
+    for met_id in transportedMets:
+        if met_id not in WATER:
+            out_comp = transportedMets[met_id]['reactant'].compartment
+            in_comp = transportedMets[met_id]['product'].compartment
             mem_pot = compartmentsData[out_comp]['membranePot'][in_comp]
-            charge = transportedMets[seed_id]['reactant'].thermo.charge_std
+            charge = transportedMets[met_id]['reactant'].thermo.charge_std
             # Equal to the product's one
             sum_F_memP_charge += (faraday_const
                                   * (mem_pot / 1000.)
-                                  * transportedMets[seed_id]['coeff']
+                                  * transportedMets[met_id]['coeff']
                                   * charge)
 
     deltaG = 0
 
     for met in reaction.metabolites:
-        if CPD_PROTON != met.annotation['seed_id']:
+        if met.annotation[annotation_key] not in PROTON:
             deltaG += reaction.metabolites[met] * met.thermo.deltaGf_tr
 
     sum_deltaGFis = 0
@@ -133,14 +133,15 @@ def calcDGtpt_rhs(reaction, compartmentsData, thermo_units):
 
     final_coeffs = reaction.metabolites.copy()
 
-    for seed_id in transportedMets:
+    for met_id in transportedMets:
         for metType in ['reactant', 'product']:
-            final_coeffs[transportedMets[seed_id][metType]] -= (
+            final_coeffs[transportedMets[met_id][metType]] -= (
                 (1 if metType == 'product' else -1)
-                * transportedMets[seed_id]['coeff'])
+                * transportedMets[met_id]['coeff'])
 
     for met in final_coeffs:
-        if final_coeffs[met] != 0 and met.annotation['seed_id'] != CPD_PROTON:
+        if (final_coeffs[met] != 0 and
+                met.annotation[annotation_key] not in PROTON):
 
             met_deltaGis = met.thermo.deltaGf_tr
             sum_deltaGFis += final_coeffs[met] * met_deltaGis
@@ -246,4 +247,4 @@ def get_debye_huckel_b(T):
     :param T: Temperature in Kelvin
     :return: Debye_Huckel_B
     """
-    return std.DEBYE_HUCKEL_B_0
+    return std.DEBYE_HUCKEL_B_0

pytfa/thermo/tmodel.py

@@ -27,6 +27,7 @@
     check_transport_reaction,
     find_transported_mets,
     get_reaction_compartment,
+    PROTON
 )
 from ..optim.constraints import (
     SimultaneousUse,
@@ -69,7 +70,8 @@ class ThermoModel(LCSBModel, Model):
     def __init__(self, thermo_data=None, model=Model(), name=None,
                  temperature=std.TEMPERATURE_0,
                  min_ph=std.MIN_PH,
-                 max_ph=std.MAX_PH):
+                 max_ph=std.MAX_PH,
+                 annotation_key="seed_id"):
 
         """
         :param float temperature: the temperature (K) at which to perform the calculations
@@ -84,6 +86,7 @@ def __init__(self, thermo_data=None, model=Model(), name=None,
         self.TEMPERATURE = temperature
         self.thermo_data = thermo_data
         self.parent = model
+        self.annotation_key = annotation_key
 
         # CONSTANTS
         self.MAX_pH = max_ph
@@ -150,21 +153,23 @@ def _prepare_metabolite(self, met):
         CompartmentionicStr = self.compartments[met.compartment]["ionicStr"]
 
         # Which index of the reaction DB do you correspond to ?
-        if not "seed_id" in met.annotation:
+        if not self.annotation_key in met.annotation:
             # raise Exception("seed_id missing for " + met.name)
             self.logger.debug(
-                "Metabolite {} ({}) has no seed_id".format(met.id, met.name)
+                "Metabolite {} ({}) has no {}".format(
+                    met.id, met.name, self.annotation_key
+                )
             )
             metData = None
-        elif not met.annotation["seed_id"] in self.compounds_data:
+        elif not met.annotation[self.annotation_key] in self.compounds_data:
             self.logger.debug(
                 "Metabolite {} ({}) not present in thermoDB".format(
-                    met.annotation["seed_id"], met.name
+                    met.annotation[self.annotation_key], met.name
                 )
             )
             metData = None
         else:
-            metData = self.compounds_data[met.annotation["seed_id"]]
+            metData = self.compounds_data[met.annotation[self.annotation_key]]
             # Override the formula
             met.formula = metData["formula"]
 
@@ -216,7 +221,10 @@ def _prepare_reaction(self, reaction, null_error_override=2):
         )
 
         # Also test if this is a transport reaction
-        reaction.thermo["isTrans"] = check_transport_reaction(reaction)
+        reaction.thermo["isTrans"] = check_transport_reaction(
+            reaction,
+            self.annotation_key
+        )
         # Make sure we have correct thermo values for each metabolites
         correctThermoValues = True
 
@@ -247,7 +255,10 @@ def _prepare_reaction(self, reaction, null_error_override=2):
 
             if reaction.thermo["isTrans"]:
                 (rhs, breakdown) = calcDGtpt_rhs(
-                    reaction, self.compartments, self.thermo_unit
+                    reaction,
+                    self.compartments,
+                    self.thermo_unit,
+                    self.annotation_key
                 )
 
                 reaction.thermo["deltaGR"] = rhs
@@ -256,9 +267,9 @@ def _prepare_reaction(self, reaction, null_error_override=2):
             else:
                 for met in reaction.metabolites:
                     if met.formula != "H" or (
-                        "seed_id" in met.annotation
+                        self.annotation_key in met.annotation
                         # That's H+
-                        and met.annotation["seed_id"] != "cpd00067"
+                        and met.annotation[self.annotation_key] != "cpd00067"
                     ):
                         DeltaGrxn += (
                             reaction.metabolites[met] * met.thermo.deltaGf_tr
@@ -311,8 +322,8 @@ def prepare(self, null_error_override=2):
         proton = {}
         for i in range(num_mets):
             if self.metabolites[i].formula == "H" or (
-                "seed_id" in self.metabolites[i].annotation
-                and self.metabolites[i].annotation["seed_id"] == "cpd00067"
+                self.annotation_key in self.metabolites[i].annotation
+                and self.metabolites[i].annotation[self.annotation_key] in PROTON
             ):
                 proton[self.metabolites[i].compartment] = self.metabolites[i]
 
@@ -363,8 +374,8 @@ def _convert_metabolite(self, met, add_potentials, verbose):
             )
 
         elif (
-            "seed_id" in met.annotation
-            and met.annotation["seed_id"] == "cpd11416"
+            self.annotation_key in met.annotation
+            and met.annotation[self.annotation_key] == "cpd11416"
         ):
             # we do not create the thermo variables for biomass enzyme
             pass
@@ -421,7 +432,7 @@ def _convert_reaction(self, rxn, add_potentials, add_displacement, verbose):
         H2OtRxns = False
         if rxn.thermo["isTrans"] and len(rxn.reactants) == 1:
             try:
-                if rxn.reactants[0].annotation["seed_id"] == "cpd00001":
+                if rxn.reactants[0].annotation[self.annotation_key] == "cpd00001":
                     H2OtRxns = True
             except KeyError:
                 pass
@@ -460,7 +471,7 @@ def _convert_reaction(self, rxn, add_potentials, add_displacement, verbose):
                 # enzyme. This will be added to the constraint on the Right
                 # Hand Side (RHS)
 
-                transportedMets = find_transported_mets(rxn)
+                transportedMets = find_transported_mets(rxn, self.annotation_key)
 
                 # Chemical coefficient, it is the enzyme's coefficient...
                 # + transport coeff for reactants
@@ -564,6 +575,7 @@ def _convert_reaction(self, rxn, add_potentials, add_displacement, verbose):
                     "generating only use constraints for drain reaction"
                     + rxn.id
                 )
+                pass
             else:
                 self.logger.debug(
                     "generating only use constraints for reaction" + rxn.id
@@ -630,7 +642,7 @@ def convert(
             for reaction in self.reactions:
                 if not "isTrans" in reaction.thermo:
                     reaction.thermo["isTrans"] = check_transport_reaction(
-                        reaction
+                        reaction, self.annotation_key
                     )
         except:
             raise Exception(