diff --git a/pytfa/core/model.py b/pytfa/core/model.py index e46f33d..0593c59 100644 --- a/pytfa/core/model.py +++ b/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 diff --git a/pytfa/redgem/lumpgem.py b/pytfa/redgem/lumpgem.py index a07d8ec..cd1818a 100644 --- a/pytfa/redgem/lumpgem.py +++ b/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: diff --git a/pytfa/thermo/metabolite.py b/pytfa/thermo/metabolite.py index 2cbd82d..b01df77 100644 --- a/pytfa/thermo/metabolite.py +++ b/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) \ No newline at end of file + return (deltaGspA, sp_charge, sp_nH) diff --git a/pytfa/thermo/reaction.py b/pytfa/thermo/reaction.py index 18af129..2c899e2 100644 --- a/pytfa/thermo/reaction.py +++ b/pytfa/thermo/reaction.py @@ -14,8 +14,7 @@ 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 # @@ -23,7 +22,8 @@ -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,13 +69,13 @@ 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 @@ -83,22 +83,22 @@ def calcDGtpt_rhs(reaction, compartmentsData, thermo_units): 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 @@ -106,22 +106,22 @@ def calcDGtpt_rhs(reaction, compartmentsData, thermo_units): # 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 \ No newline at end of file + return std.DEBYE_HUCKEL_B_0 diff --git a/pytfa/thermo/tmodel.py b/pytfa/thermo/tmodel.py index d72f319..8aa9d94 100644 --- a/pytfa/thermo/tmodel.py +++ b/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( diff --git a/pytfa/thermo/utils.py b/pytfa/thermo/utils.py index 96e1044..a8fde5a 100644 --- a/pytfa/thermo/utils.py +++ b/pytfa/thermo/utils.py @@ -13,6 +13,34 @@ Formula_regex = re.compile("([A-Z][a-z]*)([0-9]*)") +# idenfifiers of water metabolite in different databases +PROTON = { + "h", "15378", "10744", "13357", "5584", "24636", "29233", "29234", + "MNXM104313", "MNXM113751", "MNXM145872", "MNXM89553", "HMDB59597", + "C00080", "PROTON", "1132304", "113529", "1470067", "156540", "163953", + "193465", "194688", "2000349", "2872447", "351626", "372511", "374900", + "425969", "425978", "425999", "427899", "428040", "428548", "5668577", + "70106", "74722", "39", "cpd00067", "MNXM1" +} + +# idenfifiers of water metabolite in different databases +WATER = { + "h2o", "oh1", "15377", "10743", "13352", "27313", "42043", "42857", + "43228", "44292", "44701", "44819", "5585", "16234", "13365", "13419", + "44641", "5594", "29356", "29374", "29375", "29412", "30490", "33806", + "33811", "33813", "41981", "29373", "41979", "MNXM114710", "MNXM114753", + "MNXM11838", "MNXM124004", "MNXM124324", "MNXM124831", "MNXM125045", + "MNXM126600", "MNXM128935", "MNXM131091", "MNXM145357", "MNXM49218", + "MNXM56889", "MNXM89551", "5882df9c-dae1-4d80-a40e-db4724271456\/compound\/969d0227-3069-4e44-9525-7ae7bad84170", + "650babc9-9d68-4b73-9332-11972ca26f7b\/compound\/799908db-b8c9-4982-86cb-1f225e2ad08c", + "650babc9-9d68-4b73-9332-11972ca26f7b\/compound\/e7f34a8e-cded-4793-b6d5-792335b38636", + "HMDB02111", "C00001", "D00001", "C01328", "C18714", "D03703", "D06322", + "CPD-15815", "HYDROXYL-GROUP", "OH", "OXONIUM", "WATER", "109276", + "1130930", "113518", "113519", "113521", "141343", "1605715", "189422", + "2022884", "29356", "351603", "5278291", "5668574", "5693747", "8851517", + "40", "cpd00001", "cpd15275", "cpd27222", "MNXM2" +} + def check_reaction_balance(reaction, proton = None): """ @@ -105,7 +133,7 @@ def check_reaction_balance(reaction, proton = None): return 'missing atoms' -def find_transported_mets(reaction): +def find_transported_mets(reaction, annotation_key="seed_id"): """ Get a list of the transported metabolites of the reaction. @@ -135,7 +163,7 @@ def find_transported_mets(reaction): reactants_coeffs = {} for met in reaction.reactants: - reactants_coeffs[met.annotation['seed_id']] = { + reactants_coeffs[met.annotation[annotation_key]] = { 'coeff':reaction.metabolites[met], 'met':met } @@ -144,7 +172,7 @@ def find_transported_mets(reaction): trans_coeffs = {} for met in reaction.products: - seed_id = met.annotation['seed_id'] + seed_id = met.annotation[annotation_key] # If the seed_id also corresponds to a reactant, we add it to the result if seed_id in reactants_coeffs: @@ -160,7 +188,7 @@ def find_transported_mets(reaction): return trans_coeffs -def check_transport_reaction(reaction): +def check_transport_reaction(reaction, annotation_key="seed_id"): """ Check if a reaction is a transport reaction @@ -179,10 +207,10 @@ def check_transport_reaction(reaction): seed_ids = {} try: for reactant in reaction.reactants: - seed_ids[reactant.annotation['seed_id']] = True + seed_ids[reactant.annotation[annotation_key]] = True for product in reaction.products: - if product.annotation['seed_id'] in seed_ids: + if product.annotation[annotation_key] in seed_ids: return True except KeyError: reactants_ids = [x.id.replace(x.compartment,'') for x in reaction.reactants] diff --git a/tests/test_core.py b/tests/test_core.py index f786be1..bf16f8b 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -18,8 +18,11 @@ import pytest import os from pytfa.utils import numerics -from settings import tmodel, this_directory, objective_value - +from settings import ( + tmodel, this_directory, objective_value, thermo_data, small_tmodel, + small_model, lexicon, compartment_data, annotate_from_lexicon, + apply_compartment_data +) # Minimal relative difference between two values to make a test fail test_precision = 1 * 10 ** -5 @@ -134,6 +137,33 @@ def test_reactions_values(reaction): else: assert(relative_error(model_rxn.thermo[thermoval], refval) < test_precision) + +def test_alternative_annotation(): + """Check that an alternative annotation key can be used.""" + alt_small_model = small_model.copy() + other_small_tmodel = pytfa.ThermoModel( + thermo_data, alt_small_model, annotation_key="alt_annot" + ) + annotate_from_lexicon(other_small_tmodel, lexicon) + apply_compartment_data(other_small_tmodel, compartment_data) + + other_small_tmodel.solver = 'optlang-glpk' + for met in other_small_tmodel.metabolites: + if "seed_id" in met.annotation: + annot = met.annotation["seed_id"] + met.annotation["alt_annot"] = annot + del met.annotation["seed_id"] + other_small_tmodel.prepare() + other_small_tmodel.convert() + sol = other_small_tmodel.optimize() + expected_sol = small_tmodel.optimize() + assert pytest.approx(sol.objective_value) == pytest.approx( + expected_sol.objective_value + ) + assert pytest.approx(sol.fluxes.sum()) == pytest.approx( + expected_sol.fluxes.sum() + ) + ############ # LP FILES # ############