Commit 2d4f192b authored by eladnoor's avatar eladnoor

adding serialize function to Reaction

parent c602faa2
Pipeline #48438902 failed with stage
in 18 minutes and 50 seconds
......@@ -39,10 +39,10 @@ from equilibrator_cache.thermodynamic_constants import (
ccache = CompoundCache()
default_pMg = Q_("14.0")
default_phase = 'aqueous'
default_conc_lb = Q_("1e-6 M")
default_conc_ub = Q_("1e-2 M")
default_e_potential = Q_("0 V")
def strip_units(v: np.array) -> np.array:
return np.array(list(map(float, v.flat))).reshape(v.shape)
......
......@@ -68,34 +68,41 @@ class ComponentContribution(object):
"""Get the value of RT."""
return R * self.temperature
def standard_dg_prime_multi(
self,
reactions: List[PhasedReaction]
) -> Tuple[np.array, np.array]:
"""Calculate the transformed reaction energies of a list of reactions.
@staticmethod
def standard_dg(reaction: PhasedReaction) -> Tuple[float, float]:
"""Calculate the chemical reaction energies of a reaction.
:return: a tuple of two arrays. the first is a 1D NumPy array
containing the CC estimates for the reactions' transformed dG0
the second is a 2D numpy array containing the covariance matrix
of the standard errors of the estimates. one can use the eigenvectors
of the matrix to define a confidence high-dimensional space, or use
U as the covariance of a Gaussian used for sampling
(where dG0_cc is the mean of that Gaussian).
:param reaction: the input Reaction object
:return: a tuple of the dG0 in kJ/mol and standard error. to
calculate the confidence interval, use the range -1.96 to 1.96 times
this value
"""
rxn_dg_pairs = map(lambda r: r.separate_stored_dgf_values(
p_h=self.p_h, ionic_strength=self.ionic_strength,
temperature=self.temperature), reactions)
residual_reactions, stored_dg_primes = zip(*rxn_dg_pairs)
stored_dg_primes = np.array(stored_dg_primes)
residual_reaction, stored_dg = reaction.separate_stored_dg()
standard_dg_prime, dg_sigma = \
ComponentContribution.predictor.standard_dg_prime_multi(
residual_reactions,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
return standard_dg_prime + stored_dg_primes, dg_sigma
standard_dg, dg_sigma = \
ComponentContribution.predictor.standard_dg(
residual_reaction)
return standard_dg + stored_dg, dg_sigma
@staticmethod
def standard_dg_multi(
reactions: List[PhasedReaction]) -> Tuple[np.array, np.array]:
"""Calculate the chemical reaction energies of a list of reactions.
Using the major microspecies of each of the reactants.
:return: a tuple with the CC estimation of the major microspecies'
standard formation energy, and the uncertainty matrix.
"""
rxn_dg_pairs = map(lambda r: r.separate_stored_dg(), reactions)
residual_reactions, stored_dg = zip(*rxn_dg_pairs)
stored_dg = np.array(stored_dg)
standard_dg, dg_sigma = \
ComponentContribution.predictor.standard_dg_multi(
residual_reactions)
return standard_dg + stored_dg, dg_sigma
def standard_dg_prime(
self,
......@@ -109,7 +116,7 @@ class ComponentContribution(object):
this value
"""
residual_reaction, stored_dg_prime = \
reaction.separate_stored_dgf_values(
reaction.separate_stored_dg_prime(
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
......@@ -141,6 +148,35 @@ class ComponentContribution(object):
dg_prime += self.RT * reaction.dg_correction()
return dg_prime, dg_uncertainty
def standard_dg_prime_multi(
self,
reactions: List[PhasedReaction]
) -> Tuple[np.array, np.array]:
"""Calculate the transformed reaction energies of a list of reactions.
:return: a tuple of two arrays. the first is a 1D NumPy array
containing the CC estimates for the reactions' transformed dG0
the second is a 2D numpy array containing the covariance matrix
of the standard errors of the estimates. one can use the eigenvectors
of the matrix to define a confidence high-dimensional space, or use
U as the covariance of a Gaussian used for sampling
(where dG0_cc is the mean of that Gaussian).
"""
rxn_dg_pairs = map(lambda r: r.separate_stored_dg_prime(
p_h=self.p_h, ionic_strength=self.ionic_strength,
temperature=self.temperature), reactions)
residual_reactions, stored_dg_primes = zip(*rxn_dg_pairs)
stored_dg_primes = np.array(stored_dg_primes)
standard_dg_prime, dg_sigma = \
ComponentContribution.predictor.standard_dg_prime_multi(
residual_reactions,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
return standard_dg_prime + stored_dg_primes, dg_sigma
def physiological_dg_prime(
self,
reaction: PhasedReaction
......
......@@ -26,6 +26,7 @@
# THE SOFTWARE.
from collections import namedtuple
from typing import Tuple
import numpy as np
from equilibrator_cache import Compound
......@@ -40,21 +41,17 @@ GAS_PHASE_NAME = 'gas'
LIQUID_PHASE_NAME = 'liquid'
SOLID_PHASE_NAME = 'solid'
PHASE_SUBSCRIPT_TO_NAME = {'aq': AQUEOUS_PHASE_NAME, 's': SOLID_PHASE_NAME,
'l': LIQUID_PHASE_NAME, 'g': GAS_PHASE_NAME}
PhaseInfo = namedtuple("phase_info", "subscript "
PhaseInfo = namedtuple("phase_info", "shorthand "
"standard_abundance "
"physiolical_abundance "
"dimensionality")
PHASE_INFO_DICT = {
AQUEOUS_PHASE_NAME: PhaseInfo("aq", Q_("1 M"), Q_("1 mM"),
AQUEOUS_PHASE_NAME: PhaseInfo("(aq)", Q_("1 M"), Q_("1 mM"),
"[concentration]"),
GAS_PHASE_NAME: PhaseInfo("g", Q_("1 bar"), Q_("1 mbar"), "[pressure]"),
SOLID_PHASE_NAME: PhaseInfo("s", None, None, None),
LIQUID_PHASE_NAME: PhaseInfo("l", None, None, None)
GAS_PHASE_NAME: PhaseInfo("(g)", Q_("1 bar"), Q_("1 mbar"), "[pressure]"),
SOLID_PHASE_NAME: PhaseInfo("(s)", None, None, None),
LIQUID_PHASE_NAME: PhaseInfo("(l)", None, None, None)
}
# dictionary of compounds that can be non aqueous. The values are tuples of
......@@ -133,18 +130,6 @@ class Condition(object):
"""
return PHASE_INFO_DICT[self._phase].dimensionality
@staticmethod
def get_default(compound: Compound):
"""Get the default phase of a compound.
:param compound: a Compound
:return: the default phase
"""
if compound.mnx_id in NON_AQUEOUS_COMPOUND_DICT:
return Condition(NON_AQUEOUS_COMPOUND_DICT[compound.mnx_id][0])
else:
return Condition(AQUEOUS_PHASE_NAME)
@property
def ln_abundance(self) -> float:
"""Return the log of the ratio between given and std abundances."""
......@@ -212,18 +197,45 @@ class PhasedCompound(object):
def __init__(self, compound: Compound, condition: Condition = None):
"""Create a new PhasedCompound object."""
self.compound = compound
self.condition = condition or Condition.get_default(compound)
self.condition = condition or Condition(self.possible_phases[0])
@staticmethod
def get_default(compound: Compound) -> Condition:
"""Get the default phase of a compound.
:param compound: a Compound
:return: the default phase
"""
if compound.mnx_id in NON_AQUEOUS_COMPOUND_DICT:
return Condition(NON_AQUEOUS_COMPOUND_DICT[compound.mnx_id][0])
else:
return Condition(AQUEOUS_PHASE_NAME)
@property
def mnx_id(self):
def mnx_id(self) -> str:
"""Get the compound's MetaNetX ID."""
return self.compound.mnx_id
@property
def phase(self):
def formula(self) -> str:
"""Get the chemical formula."""
return self.compound.formula
@property
def phase(self) -> str:
"""Get the phase."""
return self.condition.phase
@property
def phase_shorthand(self) -> str:
"""Get the phase shorthand (i.e. 'l' for liquid)."""
if self.compound.mnx_id == "MNXM60":
# for CO2(total) there should be no indication of the phase,
# since it is actually a combination of gas and aqueous.
return ""
else:
return PHASE_INFO_DICT[self.condition.phase].shorthand
@phase.setter
def phase(self, phase: str) -> None:
"""Change the phase of a specific compound.
......@@ -231,19 +243,19 @@ class PhasedCompound(object):
:param phase: the new phase
:return:
"""
if self.mnx_id in NON_AQUEOUS_COMPOUND_DICT:
possible_phases = NON_AQUEOUS_COMPOUND_DICT[self.mnx_id]
assert phase in possible_phases, (
f"The phase of {self.compound} must be one of the following: "
f"{str(possible_phases)}."
)
else:
assert phase == AQUEOUS_PHASE_NAME, (
f"The phase of {self.compound} can only by aqueous, "
f"not {phase}."
)
assert phase in self.possible_phases, (
f"The phase of {self.compound} must be one of the following: "
f"{str(self.possible_phases)}."
)
self.condition = Condition(phase)
@property
def possible_phases(self) -> Tuple[str]:
"""Get the possible phases for this compound."""
return NON_AQUEOUS_COMPOUND_DICT.get(
self.mnx_id, (AQUEOUS_PHASE_NAME, )
)
@property
def abundance(self) -> ureg.Quantity:
"""Get the abundance."""
......@@ -288,10 +300,10 @@ class PhasedCompound(object):
@ureg.check(None, None, "[concentration]", "[temperature]")
def get_stored_standard_dgf_prime(
self,
p_h: ureg.Quantity,
ionic_strength: ureg.Quantity,
temperature: ureg.Quantity
self,
p_h: ureg.Quantity,
ionic_strength: ureg.Quantity,
temperature: ureg.Quantity
) -> ureg.Quantity:
"""Return the stored formation energy of this phased compound.
......@@ -303,7 +315,7 @@ class PhasedCompound(object):
:param temperature:
:return: standard_dgf_prime (in kJ/mol)
"""
ms = PHASED_COMPOUND_DICT.get((self.mnx_id, self.phase), None)
ms = self.get_stored_microspecie()
if ms is not None:
# for compound-phase pairs in this dictionary use the stored
# value and transform directly using the Legendre transform
......@@ -315,3 +327,53 @@ class PhasedCompound(object):
else:
# for all other compounds, we will use component_contribution
return None
def get_stored_standard_dgf(
self
) -> ureg.Quantity:
"""Return the stored formation energy of this phased compound.
Only if it exists, otherwise return None (and we will use
component-contribution later to get the reaction energy).
:return: standard_dgf (in kJ/mol)
"""
ms = self.get_stored_microspecie()
if ms is not None:
# for compound-phase pairs in this dictionary use the stored
# value
return ms.standard_dgf
else:
# for all other compounds, we will use component_contribution
return None
def get_stored_microspecie(self) -> MicroSpecie:
"""Get the stored microspecies (from the PHASED_COMPOUND_DICT).
:return: The MicroSpecie namedtuple with the stored formation energy,
or None if this compound has no stored value at this phase.
"""
return PHASED_COMPOUND_DICT.get((self.mnx_id, self.phase), None)
def serialize(self) -> dict:
"""Return a serialized version of all the compound thermo data.
:return: a list of dictionaries with all the microspecies data
"""
ms_list = []
ms = self.get_stored_microspecie()
if ms is not None:
ms_list.append(ms.__dict__())
else:
for ms in self.compound.microspecies:
d = {"num_protons": ms.number_protons,
"charge": ms.charge,
"ddg_over_rt": ms.ddg_over_rt,
"num_magnesiums": 0
}
ms_list.append(d)
return {"mnx_id": self.mnx_id,
"phase": self.phase,
"ln_abundance": self.ln_abundance,
"microspecies": ms_list}
......@@ -27,7 +27,7 @@
import logging
from typing import Dict, Tuple
from typing import Dict, List, Tuple
import numpy as np
from equilibrator_cache.reaction import Reaction
......@@ -87,6 +87,34 @@ class PhasedReaction(Reaction):
if phased_compound.compound == compound:
phased_compound.phase = phase
def get_phased_compound(
self,
compound: Compound
) -> Tuple[PhasedCompound, float]:
"""Get the PhasedCompound object by the Compound object."""
# TODO: This is not ideal. We should try to not keep to dictionaries (
# one with Compounds and one with PhasedCompounds).
for phased_compound, coeff in self.sparse_with_phases.items():
if phased_compound.compound == compound:
return phased_compound, coeff
return None, 0
def get_phase(self, compound: Compound) -> str:
"""Get the phase of the compound."""
phased_compound, _ = self.get_phased_compound(compound)
if phased_compound is None:
return ""
else:
return phased_compound.phase
def get_abundance(self, compound: Compound) -> ureg.Quantity:
"""Get the abundance of the compound."""
phased_compound, _ = self.get_phased_compound(compound)
if phased_compound is None:
return None
else:
return phased_compound.abundance
@property
def is_physiological(self) -> bool:
"""Check if all concentrations are physiological.
......@@ -102,6 +130,13 @@ class PhasedReaction(Reaction):
return False
return True
def get_stoichiometry(self, compound: Compound) -> float:
"""Get the abundance of the compound."""
for phased_compound, coeff in self.sparse_with_phases.items():
if phased_compound.compound == compound:
return coeff
return 0.0
def add_stoichiometry(self,
cpd: Compound,
coeff: float) -> None:
......@@ -119,7 +154,7 @@ class PhasedReaction(Reaction):
self.sparse_with_phases[PhasedCompound(cpd)] = coeff
@ureg.check(None, None, "[concentration]", "[temperature]")
def separate_stored_dgf_values(
def separate_stored_dg_prime(
self,
p_h: ureg.Quantity,
ionic_strength: ureg.Quantity,
......@@ -148,11 +183,31 @@ class PhasedReaction(Reaction):
return Reaction(sparse, arrow=self.arrow, rid=self.rid), stored_dg_prime
def separate_stored_dg(
self,
) -> Tuple[Reaction, ureg.Quantity]:
"""Split the PhasedReaction to aqueous phase and all the rest.
:return: a tuple (residual_reaction, stored_dg) where
residual_reaction is a Reaction object (excluding the compounds that
had stored values), and stored_dg is the total dG of the
compounds with stored values (in kJ/mol).
"""
stored_dg_prime = 0
sparse = {} # all compounds without stored dgf values
for phased_compound, coeff in self.sparse_with_phases.items():
dgf_prime = phased_compound.get_stored_standard_dgf()
if dgf_prime is None:
sparse[phased_compound.compound] = coeff
else:
stored_dg_prime += coeff * dgf_prime
return Reaction(sparse, arrow=self.arrow, rid=self.rid), stored_dg_prime
def dg_correction(self) -> float:
"""Calculate the concentration adjustment in the dG' of reaction.
:param compound_to_conc: a dictionary mapping Compound objects to
concentration in M (default is 1M)
:return: the correction for delta G in units of RT
"""
# here we check that each concentration is in suitable units,
......@@ -180,18 +235,23 @@ class PhasedReaction(Reaction):
By adding H2O, O2, Pi, CO2, and NH4+ to both sides.
"""
balancing_ids = ['MNX:MNXM2', # H2O
# We need to make sure that the number of balancing compounds is the
# same as the number of elements we are trying to balance. Here both
# numbers are 6, i.e. the elements are (e-, H, O, P, C, N)
balancing_ids = ['MNX:MNXM1', # H+
'MNX:MNXM2', # H2O
'MNX:MNXM4', # O2
'MNX:MNXM9', # Pi
'MNX:MNXM13', # CO2
'MNX:MNXM15', # NH4+
]
S = ccache.get_element_data_frame(balancing_ids)
S = ccache.get_element_data_frame(
map(ccache.get_compound, balancing_ids)
)
balancing_atoms = S.index
atom_bag = self._get_reaction_atom_balance()
atom_bag = self._get_reaction_atom_bag()
if atom_bag is None:
logging.warning('Cannot balance this reaction due to'
' missing chemical formulas')
......@@ -204,16 +264,17 @@ class PhasedReaction(Reaction):
raise ValueError(f"Cannot oxidize {self} only with these atoms: "
f"{other_atoms}")
imbalance = np.linalg.inv(S) @ atom_vector
# solve the linear equation S * i = a,
# and ignore small coefficients (< 1e-3)
imbalance = (-np.linalg.inv(S.values) @ atom_vector).round(3)
for compound, coeff in zip(map(ccache.get_compound, balancing_ids),
imbalance.flat):
self.sparse[compound] = self.sparse.get(compound, 0) - coeff
for compound, coeff in zip(S.columns, imbalance.flat):
self.add_stoichiometry(compound, coeff)
return self
@classmethod
def get_oxidation_reaction(cls, compound: Compound) -> object:
def get_oxidation_reaction(cls, compound: Compound):
"""Generate an oxidation Reaction for a single compound.
Generate a Reaction object which represents the oxidation reaction
......@@ -223,29 +284,11 @@ class PhasedReaction(Reaction):
"""
return cls({compound: -1}).balance_by_oxidation()
if __name__ == '__main__':
import sys
r = PhasedReaction.get_oxidation_reaction('C00031')
print(r.write_formula())
print('standard oxidation energy of glucose: %.2f kJ/mol' % r.dG0_prime())
r = PhasedReaction.get_oxidation_reaction('C00064')
print(r.write_formula())
print('standard oxidation energy of acetate: %.2f kJ/mol' % r.dG0_prime())
r = PhasedReaction.parse_formula('C00031 = ').balance_by_oxidation()
print(r.write_formula())
print('standard oxidation energy of glucose: %.2f kJ/mol' % r.dG0_prime())
print('oxidation energy of 1 mM glucose: %.2f kJ/mol' %
r.dG_prime({'C00031': 1e-3}))
print('\nNow, trying to use a compound with an unspecific formula:')
sys.stdout.flush()
r = PhasedReaction.get_oxidation_reaction('C04619')
print(r.write_formula())
print('standard oxidation energy of '
'(3R)-3-Hydroxydecanoyl-[acyl-carrier protein]: %s kJ/mol'
% r.dG0_prime())
def serialize(self) -> List[dict]:
"""Return a serialized version of all the reaction thermo data."""
list_of_compounds = []
for phased_compound, coeff in self.sparse_with_phases.items():
compound_dict = phased_compound.serialize()
compound_dict["coefficient"] = coeff
list_of_compounds.append(compound_dict)
return list_of_compounds
"""unit test for balance with oxidation function."""
# The MIT License (MIT)
#
# Copyright (c) 2013 Weizmann Institute of Science
# Copyright (c) 2018 Institute for Molecular Systems Biology,
# ETH Zurich
# Copyright (c) 2018 Novo Nordisk Foundation Center for Biosustainability,
# Technical University of Denmark
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
import pytest
from equilibrator_api import Q_, ComponentContribution, Reaction, ccache
from . import approx_unit
@pytest.fixture(scope="module")
def comp_contribution() -> ComponentContribution:
"""Create a ComponentContribution object."""
return ComponentContribution(p_h=Q_("7"), ionic_strength=Q_("0.25M"))
@pytest.mark.parametrize(
"compound_id, exp_standard_dg_prime",
[
("kegg:C00031", Q_("-2929.5 kJ/mol")),
("kegg:C00033", Q_("-872.6 kJ/mol")),
("kegg:C00064", Q_("-2040.3 kJ/mol")),
])
def test_oxidation(compound_id, exp_standard_dg_prime, comp_contribution):
"""Test how reactions are balanced by oxidation."""
compound = ccache.get_compound(compound_id)
r1 = Reaction.get_oxidation_reaction(compound)
approx_unit(comp_contribution.dg_prime(r1)[0],
exp_standard_dg_prime,
abs=0.1)
r2 = Reaction.parse_formula(f"{compound_id} = ").balance_by_oxidation()
approx_unit(comp_contribution.dg_prime(r2)[0],
exp_standard_dg_prime,
abs=0.1)
"""unit test for balance with oxidation function."""
# The MIT License (MIT)
#
# Copyright (c) 2013 Weizmann Institute of Science
# Copyright (c) 2018 Institute for Molecular Systems Biology,
# ETH Zurich
# Copyright (c) 2018 Novo Nordisk Foundation Center for Biosustainability,
# Technical University of Denmark
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
import pytest
from equilibrator_api import Q_, Reaction, ccache
@pytest.fixture(scope="module")
def atp_hydrolysis() -> Reaction:
"""Create a ATP hydrolysis reaction."""
formula = "KEGG:C00002 + KEGG:C00001 = KEGG:C00008 + KEGG:C00009"
return Reaction.parse_formula(formula)
def test_serialize(atp_hydrolysis):
"""Test the serialize function of PhasedRection."""
atp_compound = ccache.get_compound("KEGG:C00002")
atp_hydrolysis.set_abundance(atp_compound, Q_("10 mM"))
ser = atp_hydrolysis.serialize()
assert ser[0]["mnx_id"] == "MNXM3" # ATP
assert ser[0]["phase"] == "aqueous"
assert ser[0]["coefficient"] == -1
assert len(ser[0]["microspecies"]) == 7
assert ser[1]["mnx_id"] == "MNXM2" # H2O
assert ser[1]["phase"] == "liquid"
assert len(ser[1]["microspecies"]) == 1
assert ser[2]["mnx_id"] == "MNXM7" # ADP
assert len(ser[2]["microspecies"]) == 7
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment