...
  View open merge request
Commits (5)
...@@ -35,19 +35,18 @@ The easiest way to get eQuilibrator-API up and running is using virtualenv, PyPI ...@@ -35,19 +35,18 @@ The easiest way to get eQuilibrator-API up and running is using virtualenv, PyPI
``` ```
virtualenv -p python3 equilibrator virtualenv -p python3 equilibrator
source equilibrator/bin/activate source equilibrator/bin/activate
pip install equilibrator-api jupyter pip install equilibrator-api
curl https://gitlab.com/elad.noor/equilibrator-api/raw/develop/scripts/equilibrator_cmd.ipynb > equilibrator_cmd.ipynb
jupyter notebook
``` ```
Then select the notebook called `equilibrator_cmd.ipynb` and follow the examples in it. You can then either follow the examples below and write your own code, or start by opening the
provided jupyter notebook:
Alternatively, you could install from source. Make sure you have [git-lfs](https://git-lfs.github.com/)
installed before cloning the repository:
```
git clone https://gitlab.com/elad.noor/equilibrator-api.git
cd equilibrator-api
python setup.py install
``` ```
pip install jupyter
curl https://gitlab.com/elad.noor/equilibrator-api/raw/develop/scripts/equilibrator_cmd.ipynb > equilibrator_cmd.ipynb
jupyter notebook
```
Then select the notebook called `equilibrator_cmd.ipynb` and follow the examples in it.
## Example Usage ## Example Usage
...@@ -63,7 +62,7 @@ You can parse a reaction from a KEGG-style reaction string. The example given ...@@ -63,7 +62,7 @@ You can parse a reaction from a KEGG-style reaction string. The example given
is ATP hydrolysis to ADP and inorganic phosphate. is ATP hydrolysis to ADP and inorganic phosphate.
```python ```python
rxn_str = "KEGG:C00002 + KEGG:C00001 = KEGG:C00008 + KEGG:C00009" rxn_str = "kegg:C00002 + kegg:C00001 = kegg:C00008 + kegg:C00009"
rxn = Reaction.parse_formula(rxn_str) rxn = Reaction.parse_formula(rxn_str)
``` ```
...@@ -84,14 +83,14 @@ calculate the standard change in Gibbs potential due to this reaction. ...@@ -84,14 +83,14 @@ calculate the standard change in Gibbs potential due to this reaction.
# You control the pH and ionic strength! # You control the pH and ionic strength!
# ionic strength is in Molar units. # ionic strength is in Molar units.
standard_dg_prime, uncertainty = eq_api.standard_dg_prime(rxn) standard_dg_prime, uncertainty = eq_api.standard_dg_prime(rxn)
print("dG0' = %s \u00B1 %s\n" % (standard_dg_prime, uncertainty)) print(f"dG0' = {standard_dg_prime:.1f} \u00B1 {uncertainty:1.f}")
``` ```
You can also calculate the [reversibility index](https://doi.org/10.1093/bioinformatics/bts317) for this reaction. You can also calculate the [reversibility index](https://doi.org/10.1093/bioinformatics/bts317) for this reaction.
```python ```python
ln_RI = eq_api.ln_reversibility_index(rxn) ln_RI = eq_api.ln_reversibility_index(rxn)
print('ln(Reversibility Index) = %s\n' % ln_RI) print(f"ln(Reversibility Index) = {ln_RI:.1f}")
``` ```
The reversibility index is a measure of the degree of the reversibility of the The reversibility index is a measure of the degree of the reversibility of the
......
...@@ -37,8 +37,7 @@ import pandas as pd ...@@ -37,8 +37,7 @@ import pandas as pd
def main(args): def main(args):
"""Run main script, calculates the reaction Gibbs energy change.""" """Run main script, calculates the reaction Gibbs energy change."""
from equilibrator_api import ( from equilibrator_api import (
Q_, ComponentContribution, Reaction, ureg) # isort:skip Q_, ComponentContribution, Reaction) # isort:skip
ureg.default_format = ".2f"
p_h = Q_(args.ph) p_h = Q_(args.ph)
assert p_h.check(None) assert p_h.check(None)
...@@ -115,3 +114,4 @@ if __name__ == '__main__': ...@@ -115,3 +114,4 @@ if __name__ == '__main__':
logging.getLogger().setLevel(logging.WARNING) logging.getLogger().setLevel(logging.WARNING)
args = parser.parse_args() args = parser.parse_args()
main(args)
This diff is collapsed.
...@@ -62,7 +62,6 @@ def main(args): ...@@ -62,7 +62,6 @@ def main(args):
from equilibrator_api import ( from equilibrator_api import (
Q_, ComponentContribution, Reaction, ureg) # isort:skip Q_, ComponentContribution, Reaction, ureg) # isort:skip
ureg.default_format = ".2f"
p_h = Q_(args.ph) p_h = Q_(args.ph)
assert p_h.check(None) assert p_h.check(None)
...@@ -80,7 +79,7 @@ def main(args): ...@@ -80,7 +79,7 @@ def main(args):
assert temperature.check("[temperature]") assert temperature.check("[temperature]")
sys.stderr.write(f"pH = {p_h}\n") sys.stderr.write(f"pH = {p_h}\n")
sys.stderr.write(f"I = {ionic_strength:.2g}\n") sys.stderr.write(f"I = {ionic_strength}\n")
sys.stderr.write(f"T = {temperature}\n") sys.stderr.write(f"T = {temperature}\n")
sys.stderr.flush() sys.stderr.flush()
......
...@@ -32,58 +32,103 @@ import sys ...@@ -32,58 +32,103 @@ import sys
from numpy import nan, sqrt from numpy import nan, sqrt
from equilibrator_api import ComponentContribution, Reaction
def main(args):
from equilibrator_api import (
Q_, ComponentContribution, Reaction) # isort:skip
if __name__ == '__main__': p_h = Q_(args.ph)
assert p_h.check(None)
parser = argparse.ArgumentParser( ionic_strength = Q_(args.i)
description='Calculate potentials for a number of reactions.') if ionic_strength.check(None):
parser.add_argument( ionic_strength *= Q_("M")
'outfile', type=argparse.FileType('w'), else:
help='path to output file') assert ionic_strength.check("[concentration]")
parser.add_argument('--i', type=float,
help='ionic strength in M',
default=0.1)
parser.add_argument('--ph', type=float, help='pH level', default=7.0)
logging.getLogger().setLevel(logging.WARNING)
args = parser.parse_args() temperature = Q_(args.t)
if temperature.check(None):
temperature *= Q_("K")
else:
assert temperature.check("[temperature]")
sys.stderr.write(f"pH = {p_h}\n")
sys.stderr.write(f"I = {ionic_strength}\n")
sys.stderr.write(f"T = {temperature}\n")
sys.stderr.write('pH = %.1f\n' % args.ph) sys.stderr.write("Initializing the Component Contribution estimator.\n")
sys.stderr.write('I = %.1f M\n' % args.i) eq_api = ComponentContribution(p_h=p_h, ionic_strength=ionic_strength,
temperature=temperature)
ids = [] sys.stderr.write("Parsing model reactions.\n")
reactions = []
with open('data/iJO1366_reactions.csv', 'r') as f: all_ids = []
comments = []
balanced_ids = [] # only those of real and balanced reactions
balanced_reactions = [] # a list of the reaction objects that are balanced
with open('iJO1366_reactions.csv', 'r') as f:
for row in csv.DictReader(f): for row in csv.DictReader(f):
ids.append(row['bigg.reaction']) all_ids.append(row['bigg.reaction'])
try: try:
reactions.append(Reaction.parse_formula(row['formula'])) rxn = Reaction.parse_formula(row['formula'])
sys.stderr.write(f"Reaction {row['bigg.reaction']}: {rxn}\n")
if rxn.is_empty():
comments.append("reaction is empty")
elif not rxn.is_balanced():
comments.append("reaction is not chemically balanced")
else:
balanced_ids.append(row['bigg.reaction'])
balanced_reactions.append(rxn)
comments.append("")
except ValueError as e: except ValueError as e:
print('warning: cannot parse reaction %s because of %s' % print('warning: cannot parse reaction %s because of %s' %
(row['bigg.reaction'], str(e))) (row['bigg.reaction'], str(e)))
reactions.append(Reaction({}))
continue
cc = ComponentContribution(pH=args.ph, ionic_strength=args.i) sys.stderr.write("Calculating the ΔG'0 of all reactions in the model.\n")
dG0_prime, U = eq_api.standard_dg_prime_multi(balanced_reactions)
dG0_prime, U = cc.standard_dg_prime_multi(reactions) rid_to_reaction = dict(zip(balanced_ids, balanced_reactions))
rid_to_dg0_prime = dict(zip(balanced_ids, dG0_prime))
rid_to_uncertainty = dict(zip(balanced_ids, U.diagonal().flat))
sys.stderr.write("Writing the results to .csv file.\n")
writer = csv.writer(args.outfile) writer = csv.writer(args.outfile)
header = ['reaction', 'pH', 'ionic strength [M]', 'dG\'0 [kJ/mol]', header = ["reaction", "pH", "ionic strength [M]", "ΔG'0",
'uncertainty [kJ/mol]', 'ln(Reversibility Index)', 'comment'] "uncertainty", "ln(Reversibility Index)", "comment"]
writer.writerow(header) writer.writerow(header)
for s, r, dg0, u in zip(ids, reactions,
dG0_prime.flat, U.diagonal().flat): for rid, comment in zip(all_ids, comments):
row = [s, args.ph, args.i] dg0 = nan
if r.is_empty(): sigma = nan
row += [nan, nan, nan, 'reaction is empty'] ln_RI = nan
elif r.check_full_reaction_balancing(): if rid in balanced_ids:
ln_RI = cc.reversibility_index(r) sigma = sqrt(rid_to_uncertainty[rid])
row += ['%.2f' % dg0, '%.2f' % sqrt(u), '%.2f' % ln_RI, ''] if sigma < Q_(100, "kJ/mol"):
else: dg0 = rid_to_dg0_prime[rid]
row += [nan, nan, nan, 'reaction is not chemically balanced'] ln_RI = eq_api.ln_reversibility_index(rid_to_reaction[rid])
writer.writerow(row)
writer.writerow([rid, args.ph, args.i, dg0, sigma, ln_RI, comment])
args.outfile.flush() args.outfile.flush()
if __name__ == '__main__':
parser = argparse.ArgumentParser(
description='Calculate potentials for a number of reactions.')
parser.add_argument(
'outfile', type=argparse.FileType('w'),
help='path to output file')
parser.add_argument('--ph', type=str,
help='pH level',
default="7.0")
parser.add_argument('--i', type=str,
help='ionic strength (in molar, default 0.25 M)',
default="0.25M")
parser.add_argument('--t', type=str,
help='temperature (in kalvin, default 298.15 K)',
default="298.15K")
logging.getLogger().setLevel(logging.WARNING)
args = parser.parse_args()
main(args)
This diff is collapsed.
...@@ -30,13 +30,15 @@ del get_versions ...@@ -30,13 +30,15 @@ del get_versions
import numpy as np import numpy as np
from equilibrator_cache import Compound from equilibrator_cache import CompoundCache, BaseCompound, BaseReaction
from equilibrator_cache import CompoundCache
from equilibrator_cache.thermodynamic_constants import ( from equilibrator_cache.thermodynamic_constants import (
FARADAY, POSSIBLE_REACTION_ARROWS, Q_, R, default_I, default_pH, FARADAY, POSSIBLE_REACTION_ARROWS, Q_, R, default_I, default_pH,
default_pMg, default_T, physiological_concentration, default_pMg, default_T, physiological_concentration,
standard_concentration, ureg) standard_concentration, ureg)
ureg.default_format = ".2f"
ureg.setup_matplotlib(True)
ccache = CompoundCache() ccache = CompoundCache()
default_phase = 'aqueous' default_phase = 'aqueous'
...@@ -49,6 +51,7 @@ def strip_units(v: np.array) -> np.array: ...@@ -49,6 +51,7 @@ def strip_units(v: np.array) -> np.array:
from equilibrator_api.bounds import Bounds from equilibrator_api.bounds import Bounds
from equilibrator_api.pathway import Pathway from equilibrator_api.pathway import Pathway
from equilibrator_api.phased_reaction import PhasedReaction as Reaction from equilibrator_api.phased_compound import Compound
from equilibrator_api.phased_reaction import Reaction
from equilibrator_api.component_contribution import ComponentContribution from equilibrator_api.component_contribution import ComponentContribution
from equilibrator_api.thermo_models import PathwayThermoModel from equilibrator_api.thermo_models import PathwayThermoModel
...@@ -33,7 +33,7 @@ import pkg_resources ...@@ -33,7 +33,7 @@ import pkg_resources
from numpy import log from numpy import log
from . import ( from . import (
Q_, Compound, ccache, default_conc_lb, default_conc_ub, Q_, BaseCompound, ccache, default_conc_lb, default_conc_ub,
standard_concentration, ureg) standard_concentration, ureg)
...@@ -44,7 +44,7 @@ class BaseBounds(object): ...@@ -44,7 +44,7 @@ class BaseBounds(object):
"""Return a (deep) copy of self.""" """Return a (deep) copy of self."""
raise NotImplementedError raise NotImplementedError
def GetLowerBound(self, compound: Compound): def GetLowerBound(self, compound: BaseCompound):
"""Get the lower bound for this key. """Get the lower bound for this key.
Args: Args:
...@@ -52,7 +52,7 @@ class BaseBounds(object): ...@@ -52,7 +52,7 @@ class BaseBounds(object):
""" """
raise NotImplementedError raise NotImplementedError
def GetUpperBound(self, compound: Compound): def GetUpperBound(self, compound: BaseCompound):
"""Get the upper bound for this key. """Get the upper bound for this key.
Args: Args:
...@@ -62,7 +62,7 @@ class BaseBounds(object): ...@@ -62,7 +62,7 @@ class BaseBounds(object):
def GetLowerBounds( def GetLowerBounds(
self, self,
compounds: Iterable[Compound] compounds: Iterable[BaseCompound]
) -> Iterable[float]: ) -> Iterable[float]:
"""Get the bounds for a set of keys in order. """Get the bounds for a set of keys in order.
...@@ -73,7 +73,7 @@ class BaseBounds(object): ...@@ -73,7 +73,7 @@ class BaseBounds(object):
def GetUpperBounds( def GetUpperBounds(
self, self,
compounds: Iterable[Compound] compounds: Iterable[BaseCompound]
) -> Iterable[float]: ) -> Iterable[float]:
"""Get the bounds for a set of keys in order. """Get the bounds for a set of keys in order.
...@@ -82,17 +82,17 @@ class BaseBounds(object): ...@@ -82,17 +82,17 @@ class BaseBounds(object):
""" """
return map(self.GetUpperBound, compounds) return map(self.GetUpperBound, compounds)
def GetBoundTuple(self, compound: Compound) -> Tuple[float, float]: def GetBoundTuple(self, compound: BaseCompound) -> Tuple[float, float]:
"""Get both upper and lower bounds for this key. """Get both upper and lower bounds for this key.
:param compound: a Compound object :param compound: a BaseCompound object
:return: a 2-tuple (lower bound, upper bound) :return: a 2-tuple (lower bound, upper bound)
""" """
return self.GetLowerBound(compound), self.GetUpperBound(compound) return self.GetLowerBound(compound), self.GetUpperBound(compound)
def GetBounds( def GetBounds(
self, self,
compounds: Iterable[Compound] compounds: Iterable[BaseCompound]
) -> Tuple[Iterable[float], Iterable[float]]: ) -> Tuple[Iterable[float], Iterable[float]]:
"""Get the bounds for a set of compounds. """Get the bounds for a set of compounds.
...@@ -115,7 +115,7 @@ class BaseBounds(object): ...@@ -115,7 +115,7 @@ class BaseBounds(object):
def GetLnBounds( def GetLnBounds(
self, self,
compounds: Iterable[Compound] compounds: Iterable[BaseCompound]
) -> Tuple[Iterable[float], Iterable[float]]: ) -> Tuple[Iterable[float], Iterable[float]]:
"""Get the log-bounds for a set of compounds. """Get the log-bounds for a set of compounds.
...@@ -128,7 +128,7 @@ class BaseBounds(object): ...@@ -128,7 +128,7 @@ class BaseBounds(object):
def GetLnLowerBounds( def GetLnLowerBounds(
self, self,
compounds: Iterable[Compound] compounds: Iterable[BaseCompound]
) -> Iterable[float]: ) -> Iterable[float]:
"""Get the log lower bounds for a set of compounds. """Get the log lower bounds for a set of compounds.
...@@ -140,7 +140,7 @@ class BaseBounds(object): ...@@ -140,7 +140,7 @@ class BaseBounds(object):
def GetLnUpperBounds( def GetLnUpperBounds(
self, self,
compounds: Iterable[Compound] compounds: Iterable[BaseCompound]
) -> Iterable[float]: ) -> Iterable[float]:
"""Get the log upper bounds for a set of compounds. """Get the log upper bounds for a set of compounds.
...@@ -151,7 +151,7 @@ class BaseBounds(object): ...@@ -151,7 +151,7 @@ class BaseBounds(object):
return map(self.conc2ln_conc, ubs) return map(self.conc2ln_conc, ubs)
@ureg.check(None, None, '[concentration]', '[concentration]') @ureg.check(None, None, '[concentration]', '[concentration]')
def SetBounds(self, compound: Compound, lb: float, ub: float): def SetBounds(self, compound: BaseCompound, lb: float, ub: float):
"""Set bounds for a specific key. """Set bounds for a specific key.
Args: Args:
...@@ -173,8 +173,8 @@ class Bounds(BaseBounds): ...@@ -173,8 +173,8 @@ class Bounds(BaseBounds):
@ureg.check(None, None, None, '[concentration]', '[concentration]') @ureg.check(None, None, None, '[concentration]', '[concentration]')
def __init__(self, def __init__(self,
lower_bounds: Dict[Compound, float] = None, lower_bounds: Dict[BaseCompound, float] = None,
upper_bounds: Dict[Compound, float] = None, upper_bounds: Dict[BaseCompound, float] = None,
default_lb: float = default_conc_lb, default_lb: float = default_conc_lb,
default_ub: float = default_conc_ub) -> object: default_ub: float = default_conc_ub) -> object:
"""Initialize the bounds object. """Initialize the bounds object.
...@@ -212,8 +212,8 @@ class Bounds(BaseBounds): ...@@ -212,8 +212,8 @@ class Bounds(BaseBounds):
provided provided
:return: a Bounds object :return: a Bounds object
""" """
lbs: Dict[Compound, float] = dict() lbs: Dict[BaseCompound, float] = dict()
ubs: Dict[Compound, float] = dict() ubs: Dict[BaseCompound, float] = dict()
bounds_df = pd.read_csv(f) bounds_df = pd.read_csv(f)
for row in bounds_df.itertuples(): for row in bounds_df.itertuples():
...@@ -244,7 +244,7 @@ class Bounds(BaseBounds): ...@@ -244,7 +244,7 @@ class Bounds(BaseBounds):
bounds_unit: str = None, bounds_unit: str = None,
default_lb: float = default_conc_lb, default_lb: float = default_conc_lb,
default_ub: float = default_conc_ub default_ub: float = default_conc_ub
) -> Tuple[object, Dict[Compound, str]]: ) -> Tuple[object, Dict[BaseCompound, str]]:
"""Read bounds from a Pandas DataFrame. """Read bounds from a Pandas DataFrame.
:param df: a pandas.DataFrame with the bounds :param df: a pandas.DataFrame with the bounds
...@@ -318,7 +318,7 @@ class Bounds(BaseBounds): ...@@ -318,7 +318,7 @@ class Bounds(BaseBounds):
self.default_lb, self.default_lb,
self.default_ub) self.default_ub)
def GetLowerBound(self, compound: Compound) -> float: def GetLowerBound(self, compound: BaseCompound) -> float:
"""Get the lower bound for this key. """Get the lower bound for this key.
:param compound: a compound :param compound: a compound
...@@ -326,7 +326,7 @@ class Bounds(BaseBounds): ...@@ -326,7 +326,7 @@ class Bounds(BaseBounds):
""" """
return self.lower_bounds.get(compound, self.default_lb) return self.lower_bounds.get(compound, self.default_lb)
def GetUpperBound(self, compound: Compound): def GetUpperBound(self, compound: BaseCompound):
"""Get the upper bound for this key. """Get the upper bound for this key.
:param compound: a compound :param compound: a compound
......
...@@ -32,7 +32,7 @@ import numpy as np ...@@ -32,7 +32,7 @@ import numpy as np
from component_contribution import GibbsEnergyPredictor from component_contribution import GibbsEnergyPredictor
from . import FARADAY, R, default_I, default_pH, default_pMg, default_T, ureg from . import FARADAY, R, default_I, default_pH, default_pMg, default_T, ureg
from .phased_reaction import PhasedReaction from .phased_reaction import Reaction
class ComponentContribution(object): class ComponentContribution(object):
...@@ -69,7 +69,7 @@ class ComponentContribution(object): ...@@ -69,7 +69,7 @@ class ComponentContribution(object):
return R * self.temperature return R * self.temperature
@staticmethod @staticmethod
def standard_dg(reaction: PhasedReaction) -> Tuple[float, float]: def standard_dg(reaction: Reaction) -> Tuple[float, float]:
"""Calculate the chemical reaction energies of a reaction. """Calculate the chemical reaction energies of a reaction.
:param reaction: the input Reaction object :param reaction: the input Reaction object
...@@ -87,7 +87,7 @@ class ComponentContribution(object): ...@@ -87,7 +87,7 @@ class ComponentContribution(object):
@staticmethod @staticmethod
def standard_dg_multi( def standard_dg_multi(
reactions: List[PhasedReaction]) -> Tuple[np.array, np.array]: reactions: List[Reaction]) -> Tuple[np.array, np.array]:
"""Calculate the chemical reaction energies of a list of reactions. """Calculate the chemical reaction energies of a list of reactions.
Using the major microspecies of each of the reactants. Using the major microspecies of each of the reactants.
...@@ -106,7 +106,7 @@ class ComponentContribution(object): ...@@ -106,7 +106,7 @@ class ComponentContribution(object):
def standard_dg_prime( def standard_dg_prime(
self, self,
reaction: PhasedReaction reaction: Reaction
) -> Tuple[float, float]: ) -> Tuple[float, float]:
"""Calculate the transformed reaction energies of a reaction. """Calculate the transformed reaction energies of a reaction.
...@@ -134,7 +134,7 @@ class ComponentContribution(object): ...@@ -134,7 +134,7 @@ class ComponentContribution(object):
def dg_prime( def dg_prime(
self, self,
reaction: PhasedReaction reaction: Reaction
) -> Tuple[float, float]: ) -> Tuple[float, float]:
"""Calculate the dG'0 of a single reaction. """Calculate the dG'0 of a single reaction.
...@@ -150,7 +150,7 @@ class ComponentContribution(object): ...@@ -150,7 +150,7 @@ class ComponentContribution(object):
def standard_dg_prime_multi( def standard_dg_prime_multi(
self, self,
reactions: List[PhasedReaction] reactions: List[Reaction]
) -> Tuple[np.array, np.array]: ) -> Tuple[np.array, np.array]:
"""Calculate the transformed reaction energies of a list of reactions. """Calculate the transformed reaction energies of a list of reactions.
...@@ -179,14 +179,14 @@ class ComponentContribution(object): ...@@ -179,14 +179,14 @@ class ComponentContribution(object):
def physiological_dg_prime( def physiological_dg_prime(
self, self,
reaction: PhasedReaction reaction: Reaction
) -> Tuple[float, float]: ) -> Tuple[float, float]:
"""Calculate the dG'm of a single reaction. """Calculate the dG'm of a single reaction.
Assume all aqueous reactants are at 1 mM, gas reactants at 1 mbar and Assume all aqueous reactants are at 1 mM, gas reactants at 1 mbar and
the rest at their standard concentration. the rest at their standard concentration.
:param reaction: an object of type PhasedReaction :param reaction: an object of type Reaction
:return: a tuple (dG_r_prime, dG_uncertainty) where dG_r_prime is :return: a tuple (dG_r_prime, dG_uncertainty) where dG_r_prime is
the estimated Gibbs free energy of reaction and dG_uncertainty is the the estimated Gibbs free energy of reaction and dG_uncertainty is the
standard deviation of estimation. Multiply it by 1.96 to get a 95% standard deviation of estimation. Multiply it by 1.96 to get a 95%
...@@ -198,7 +198,7 @@ class ComponentContribution(object): ...@@ -198,7 +198,7 @@ class ComponentContribution(object):
self.RT * reaction.physiological_dg_correction() self.RT * reaction.physiological_dg_correction()
return physiological_dg_prime, dg_uncertainty return physiological_dg_prime, dg_uncertainty
def ln_reversibility_index(self, reaction: PhasedReaction) -> float: def ln_reversibility_index(self, reaction: Reaction) -> float:
"""Calculate the reversibility index (ln Gamma) of a single reaction. """Calculate the reversibility index (ln Gamma) of a single reaction.
:return: ln_RI - The reversibility index (in natural log scale). :return: ln_RI - The reversibility index (in natural log scale).
...@@ -211,10 +211,10 @@ class ComponentContribution(object): ...@@ -211,10 +211,10 @@ class ComponentContribution(object):
ln_RI = (2.0 / abs_sum_coeff) * physiological_dg_prime / self.RT ln_RI = (2.0 / abs_sum_coeff) * physiological_dg_prime / self.RT
return ln_RI return ln_RI
def standard_e_prime(self, reaction: PhasedReaction) -> Tuple[float, float]: def standard_e_prime(self, reaction: Reaction) -> Tuple[float, float]:
"""Calculate the E'0 of a single half-reaction. """Calculate the E'0 of a single half-reaction.
:param reaction: a PhasedReaction object :param reaction: a Reaction object
:return: a tuple (E0_prime, E0_uncertainty) where E0_prime is :return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
the estimated standard electrostatic potential of reaction and the estimated standard electrostatic potential of reaction and
E0_uncertainty is the standard deviation of estimation. Multiply it E0_uncertainty is the standard deviation of estimation. Multiply it
...@@ -237,11 +237,11 @@ class ComponentContribution(object): ...@@ -237,11 +237,11 @@ class ComponentContribution(object):
def physiological_e_prime( def physiological_e_prime(
self, self,
reaction: PhasedReaction reaction: Reaction
) -> Tuple[float, float]: ) -> Tuple[float, float]:
"""Calculate the E'0 of a single half-reaction. """Calculate the E'0 of a single half-reaction.
:param reaction: a PhasedReaction object :param reaction: a Reaction object
:return: a tuple (E0_prime, E0_uncertainty) where E0_prime is :return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
the estimated standard electrostatic potential of reaction and the estimated standard electrostatic potential of reaction and
E0_uncertainty is the standard deviation of estimation. Multiply it E0_uncertainty is the standard deviation of estimation. Multiply it
...@@ -264,11 +264,11 @@ class ComponentContribution(object): ...@@ -264,11 +264,11 @@ class ComponentContribution(object):
def e_prime( def e_prime(
self, self,
reaction: PhasedReaction reaction: Reaction
) -> Tuple[float, float]: ) -> Tuple[float, float]:
"""Calculate the E'0 of a single half-reaction. """Calculate the E'0 of a single half-reaction.
:param reaction: a PhasedReaction object :param reaction: a Reaction object
:return: a tuple (E0_prime, E0_uncertainty) where E0_prime is :return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
the estimated standard electrostatic potential of reaction and the estimated standard electrostatic potential of reaction and
E0_uncertainty is the standard deviation of estimation. Multiply it E0_uncertainty is the standard deviation of estimation. Multiply it
...@@ -291,19 +291,19 @@ class ComponentContribution(object): ...@@ -291,19 +291,19 @@ class ComponentContribution(object):
def dg_analysis( def dg_analysis(
self, self,
reaction: PhasedReaction reaction: Reaction
) -> List[Dict[str, object]]: ) -> List[Dict[str, object]]:
"""Get the analysis of the component contribution estimation process. """Get the analysis of the component contribution estimation process.
:param reaction: a PhasedReaction object. :param reaction: a Reaction object.
:return: the analysis results as a list of dictionaries :return: the analysis results as a list of dictionaries
""" """
return ComponentContribution.predictor.get_dg_analysis(reaction) return ComponentContribution.predictor.get_dg_analysis(reaction)
def is_using_group_contribution(self, reaction: PhasedReaction) -> bool: def is_using_group_contribution(self, reaction: Reaction) -> bool:
"""Check whether group contribution is needed to get this reactions' dG. """Check whether group contribution is needed to get this reactions' dG.
:param reaction: a PhasedReaction object. :param reaction: a Reaction object.
:return: true iff group contribution is needed :return: true iff group contribution is needed
""" """
return ComponentContribution.predictor.is_using_group_contribution( return ComponentContribution.predictor.is_using_group_contribution(
......
This diff is collapsed.
...@@ -33,16 +33,15 @@ from io import IOBase ...@@ -33,16 +33,15 @@ from io import IOBase
from typing import Callable, Iterable, List, TextIO, Tuple from typing import Callable, Iterable, List, TextIO, Tuple
import numpy as np import numpy as np
from equilibrator_cache.reaction import (
create_stoichiometric_matrix_from_reactions)
from sbtab import SBtab from sbtab import SBtab
from scipy.linalg import fractional_matrix_power from scipy.linalg import fractional_matrix_power
from . import ( from . import (
Q_, Compound, R, default_I, default_pH, default_pMg, default_T, strip_units) Q_, BaseCompound, BaseReaction, R, default_I, default_pH, default_pMg,
default_T, strip_units)
from .bounds import Bounds from .bounds import Bounds
from .component_contribution import ComponentContribution from .component_contribution import ComponentContribution
from .phased_reaction import PhasedReaction from .phased_reaction import Reaction
from .thermo_models import PathwayMDFData, PathwayThermoModel from .thermo_models import PathwayMDFData, PathwayThermoModel
...@@ -53,7 +52,7 @@ class Pathway(object): ...@@ -53,7 +52,7 @@ class Pathway(object):
""" """
def __init__(self, def __init__(self,
reactions: List[PhasedReaction], reactions: List[Reaction],
fluxes: np.array, fluxes: np.array,
standard_dg_primes: np.array = None, standard_dg_primes: np.array = None,
dg_sigma: np.array = None, dg_sigma: np.array = None,
...@@ -66,7 +65,7 @@ class Pathway(object): ...@@ -66,7 +65,7 @@ class Pathway(object):
"""Initialize. """Initialize.
Args: Args:
reactions: a list of gibbs.reaction.Reaction objects. reactions: a list of Reaction objects.
fluxes: numpy.array of relative fluxes in same order as reactions. fluxes: numpy.array of relative fluxes in same order as reactions.
standard_dg_primes: reaction energies (in kJ/mol) standard_dg_primes: reaction energies (in kJ/mol)
dg_sigma: square root of the uncertainty covariance matrix dg_sigma: square root of the uncertainty covariance matrix
...@@ -82,7 +81,8 @@ class Pathway(object): ...@@ -82,7 +81,8 @@ class Pathway(object):
else: else:
self._bounds = bounds.Copy() self._bounds = bounds.Copy()
self.S = create_stoichiometric_matrix_from_reactions(reactions) self.S = BaseReaction.create_stoichiometric_matrix_from_reactions(
reactions)
self.fluxes = fluxes self.fluxes = fluxes
assert self.fluxes.shape == (Nr, ) assert self.fluxes.shape == (Nr, )
...@@ -167,7 +167,7 @@ class Pathway(object): ...@@ -167,7 +167,7 @@ class Pathway(object):
def set_compound_names( def set_compound_names(
self, self,
mapping: Callable[[Compound], str] mapping: Callable[[BaseCompound], str]
) -> None: ) -> None:
"""Use alternative compound names for outputs such as plots. """Use alternative compound names for outputs such as plots.
...@@ -206,7 +206,7 @@ class Pathway(object): ...@@ -206,7 +206,7 @@ class Pathway(object):
self._bounds.GetLnUpperBounds(self.S.index)))) self._bounds.GetLnUpperBounds(self.S.index))))
@staticmethod @staticmethod
def get_compounds(reactions: Iterable[PhasedReaction]) -> List[Compound]: def get_compounds(reactions: Iterable[Reaction]) -> List[BaseCompound]:
"""Get a unique list of all compounds in all reactions. """Get a unique list of all compounds in all reactions.
:param reactions: an iterator of reactions :param reactions: an iterator of reactions
...@@ -243,7 +243,7 @@ class Pathway(object): ...@@ -243,7 +243,7 @@ class Pathway(object):
return map(str, self.reactions) return map(str, self.reactions)
@property @property
def net_reaction(self) -> PhasedReaction: def net_reaction(self) -> Reaction:
"""Calculate the sum of all the reactions in the pathway. """Calculate the sum of all the reactions in the pathway.
:return: the net reaction :return: the net reaction
...@@ -252,7 +252,7 @@ class Pathway(object): ...@@ -252,7 +252,7 @@ class Pathway(object):
net_rxn_stoich = self.S @ v net_rxn_stoich = self.S @ v
net_rxn_stoich = net_rxn_stoich[net_rxn_stoich != 0] net_rxn_stoich = net_rxn_stoich[net_rxn_stoich != 0]
sparse = net_rxn_stoich.to_dict() sparse = net_rxn_stoich.to_dict()
return PhasedReaction(sparse) return Reaction(sparse)
@classmethod @classmethod
def from_csv_file(cls, def from_csv_file(cls,
...@@ -283,7 +283,7 @@ class Pathway(object): ...@@ -283,7 +283,7 @@ class Pathway(object):
flux = float(row.get('Flux', 0.0)) flux = float(row.get('Flux', 0.0))
logging.debug('formula = %f x (%s)', flux, rxn_formula) logging.debug('formula = %f x (%s)', flux, rxn_formula)
rxn = PhasedReaction.parse_formula(rxn_formula) rxn = Reaction.parse_formula(rxn_formula)
rxn.check_full_reaction_balancing() rxn.check_full_reaction_balancing()
reactions.append(rxn) reactions.append(rxn)
...@@ -353,9 +353,9 @@ class Pathway(object): ...@@ -353,9 +353,9 @@ class Pathway(object):
reactions = [] reactions = []
reaction_ids = [] reaction_ids = []
for row in reaction_df.itertuples(): for row in reaction_df.itertuples():
rxn = PhasedReaction.parse_formula(row.ReactionFormula, rxn = Reaction.parse_formula(row.ReactionFormula,
row.ID, row.ID,
name_to_compound.get) name_to_compound.get)
if not rxn.is_balanced(ignore_atoms=("H")): if not rxn.is_balanced(ignore_atoms=("H")):
warnings.warn(f"Reaction {rxn.rid} is not balanced") warnings.warn(f"Reaction {rxn.rid} is not balanced")
reactions.append(rxn) reactions.append(rxn)
......
"""inherit from equilibrator_cache.models.compound.Compound an add phases.""" """inherit from equilibrator_cache.BaseCompound an add phases."""
# The MIT License (MIT) # The MIT License (MIT)
# #
# Copyright (c) 2013 Weizmann Institute of Science # Copyright (c) 2013 Weizmann Institute of Science
...@@ -30,11 +30,9 @@ from collections import namedtuple ...@@ -30,11 +30,9 @@ from collections import namedtuple
from typing import Set, Tuple from typing import Set, Tuple
import numpy as np import numpy as np
from equilibrator_cache import Compound
from equilibrator_cache.reaction import Reaction
from equilibrator_cache.thermodynamic_constants import legendre_transform from equilibrator_cache.thermodynamic_constants import legendre_transform
from . import Q_, R, ccache, ureg from . import Q_, BaseCompound, R, ccache, ureg
AQUEOUS_PHASE_NAME = 'aqueous' AQUEOUS_PHASE_NAME = 'aqueous'
...@@ -189,53 +187,40 @@ class Condition(object): ...@@ -189,53 +187,40 @@ class Condition(object):
if self.physiolical_abundance is None: if self.physiolical_abundance is None:
return True # when the phase does not allow relative abundances return True # when the phase does not allow relative abundances
else: else:
return self._abundance == self.standard_abundance return self._abundance == self.physiolical_abundance
class PhasedCompound(object):
"""A class that combines a equilibrator_api Compound and a Condition."""
def __init__(self, compound: Compound, condition: Condition = None):
"""Create a new PhasedCompound object."""
self.compound = 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 @property
def mnx_id(self) -> str: def is_standard(self) -> bool:
"""Get the compound's MetaNetX ID.""" """Return True if the abundance is standard for this phase."""
return self.compound.mnx_id return (
self.standard_abundance is None or
self._abundance == self.standard_abundance
)
@property class Compound(BaseCompound):
def id(self) -> int: """A class that combines a equilibrator_api BaseCompound and a Condition."""
"""Get the compound's equilibrator internal ID."""
return self.compound.id def __init__(self, compound: BaseCompound, condition: Condition = None):
"""Create a new Compound object."""
super(Compound, self).__init__()
self.id = compound.id
self.mnx_id = compound.mnx_id
self.inchi_key = compound.inchi_key
self.inchi = compound.inchi
self.smiles = compound.smiles
self.mass = compound.mass
self.atom_bag = compound.atom_bag
self.dissociation_constants = compound.dissociation_constants
self.group_vector = compound.group_vector
self.microspecies = compound.microspecies
self.identifiers = compound.identifiers
@property self.condition = condition or Condition(self.possible_phases[0])
def formula(self) -> str:
"""Get the chemical formula."""
return self.compound.formula
@property @property
def html_formula(self) -> str: def html_formula(self) -> str:
"""Get the chemical formula.""" """Get the chemical formula."""
return re.sub(r'(\d+)', r'<sub>\1</sub>', self.compound.formula) return re.sub(r'(\d+)', r'<sub>\1</sub>', self.formula)
@property
def mass(self) -> float:
"""Get the chemical molecular mass."""
return self.compound.mass
@property @property
def phase(self) -> str: def phase(self) -> str:
...@@ -245,7 +230,7 @@ class PhasedCompound(object): ...@@ -245,7 +230,7 @@ class PhasedCompound(object):
@property @property
def phase_shorthand(self) -> str: def phase_shorthand(self) -> str:
"""Get the phase shorthand (i.e. 'l' for liquid).""" """Get the phase shorthand (i.e. 'l' for liquid)."""
if self.compound.mnx_id == "MNXM60": if self == ccache.co2_total:
# for CO2(total) there should be no indication of the phase, # for CO2(total) there should be no indication of the phase,
# since it is actually a combination of gas and aqueous. # since it is actually a combination of gas and aqueous.
return "" return ""
...@@ -260,7 +245,7 @@ class PhasedCompound(object): ...@@ -260,7 +245,7 @@ class PhasedCompound(object):
:return: :return:
""" """
assert phase in self.possible_phases, ( assert phase in self.possible_phases, (
f"The phase of {self.compound} must be one of the following: " f"The phase of {self} must be one of the following: "
f"{str(self.possible_phases)}." f"{str(self.possible_phases)}."
) )
self.condition = Condition(phase) self.condition = Condition(phase)
...@@ -284,7 +269,7 @@ class PhasedCompound(object): ...@@ -284,7 +269,7 @@ class PhasedCompound(object):
:param abundance: the new abundance in the correct units :param abundance: the new abundance in the correct units
:return: :return:
""" """
assert self.compound != Reaction.PROTON_COMPOUND, ( assert self != ccache.proton, (
f"Cannot directly change the relative abundance of H+. Use the pH " f"Cannot directly change the relative abundance of H+. Use the pH "
f"parameter instead for adjusting this value." f"parameter instead for adjusting this value."
) )
...@@ -294,7 +279,7 @@ class PhasedCompound(object): ...@@ -294,7 +279,7 @@ class PhasedCompound(object):
@property @property
def ln_abundance(self) -> float: def ln_abundance(self) -> float:
"""Return the log of the abundance (for thermodynamic calculations).""" """Return the log of the abundance (for thermodynamic calculations)."""
if self.compound == Reaction.PROTON_COMPOUND: if self == ccache.proton:
# Since pH is an environmental parameter, we do not use the # Since pH is an environmental parameter, we do not use the
# concentration of protons for energy calculations # concentration of protons for energy calculations
return 0.0 return 0.0
...@@ -303,7 +288,7 @@ class PhasedCompound(object): ...@@ -303,7 +288,7 @@ class PhasedCompound(object):
@property @property
def ln_physiological_abundance(self) -> float: def ln_physiological_abundance(self) -> float:
"""Return the log of the default physiological abundance.""" """Return the log of the default physiological abundance."""
if self.compound == Reaction.PROTON_COMPOUND: if self == ccache.proton:
# Since pH is an environmental parameter, we do not use the # Since pH is an environmental parameter, we do not use the
# concentration of protons for energy calculations # concentration of protons for energy calculations
return 0.0 return 0.0
...@@ -314,6 +299,11 @@ class PhasedCompound(object): ...@@ -314,6 +299,11 @@ class PhasedCompound(object):
"""Check if the abundance is physiological.""" """Check if the abundance is physiological."""
return self.condition.is_physiological return self.condition.is_physiological
@property
def is_standard(self) -> bool:
"""Check if the abundance is standard."""
return self.condition.is_standard
@ureg.check(None, None, "[concentration]", "[temperature]") @ureg.check(None, None, "[concentration]", "[temperature]")
def get_stored_standard_dgf_prime( def get_stored_standard_dgf_prime(
self, self,
...@@ -381,7 +371,7 @@ class PhasedCompound(object): ...@@ -381,7 +371,7 @@ class PhasedCompound(object):
if ms is not None: if ms is not None:
ms_list.append(ms.__dict__()) ms_list.append(ms.__dict__())
else: else:
for ms in self.compound.microspecies: for ms in self.microspecies:
d = {"num_protons": ms.number_protons, d = {"num_protons": ms.number_protons,
"charge": ms.charge, "charge": ms.charge,
"ddg_over_rt": ms.ddg_over_rt, "ddg_over_rt": ms.ddg_over_rt,
...@@ -398,9 +388,21 @@ class PhasedCompound(object): ...@@ -398,9 +388,21 @@ class PhasedCompound(object):
def first_name(self): def first_name(self):
"""Get the first name that appears in the compound-cache.""" """Get the first name that appears in the compound-cache."""
# TODO: use MetaNetX database to get the first name # TODO: use MetaNetX database to get the first name
return ccache.get_compound_first_name(self.compound) return ccache.get_compound_first_name(self)
@property @property
def names(self) -> Set[str]: def names(self) -> Set[str]:
"""Get the set of all associated names in the compound cache.""" """Get the set of all associated names in the compound cache."""
return ccache.get_compound_names(self.compound) return ccache.get_compound_names(self)
def __str__(self) -> str:
s = super(Compound, self).__str__()
if self.condition.is_standard:
return f"{s}{self.phase_shorthand}"
else:
return f"{s}{self.phase_shorthand} [{self.abundance:.2g}]"
def __repr__(self):
return (
f"{type(self).__name__}(mnx_id={self.mnx_id}, phase={self.phase}, abundance={self.abundance})"
)
This diff is collapsed.
...@@ -159,8 +159,6 @@ class PathwayMDFData(object): ...@@ -159,8 +159,6 @@ class PathwayMDFData(object):
:return: matplotlib Figure :return: matplotlib Figure
""" """
ureg.setup_matplotlib(True)
data_df = self.compound_df.copy() data_df = self.compound_df.copy()
data_df['y'] = np.arange(0, data_df.shape[0]) data_df['y'] = np.arange(0, data_df.shape[0])
......
...@@ -82,6 +82,9 @@ def test_atp_hydrolysis_physiological_dg(comp_contribution, atp_hydrolysis): ...@@ -82,6 +82,9 @@ def test_atp_hydrolysis_physiological_dg(comp_contribution, atp_hydrolysis):
"""Test the dG adjustments for physiological conditions (with H2O).""" """Test the dG adjustments for physiological conditions (with H2O)."""
warnings.simplefilter('ignore', ResourceWarning) warnings.simplefilter('ignore', ResourceWarning)
assert atp_hydrolysis.is_balanced() assert atp_hydrolysis.is_balanced()
assert atp_hydrolysis.get_phase(ccache.water) == 'liquid'
assert float(atp_hydrolysis.physiological_dg_correction()) == \ assert float(atp_hydrolysis.physiological_dg_correction()) == \
pytest.approx(numpy.log(1e-3), rel=1e-3) pytest.approx(numpy.log(1e-3), rel=1e-3)
......