Commit c1ebf118 authored by eladnoor's avatar eladnoor

adding unit_tests for MDF in different conditions and stdev factors

parent 56ff92fb
This diff is collapsed.
......@@ -26,7 +26,7 @@
import logging
from typing import Callable, Dict
from typing import Dict
import numpy as np
from equilibrator_cache.reaction import Reaction
......
......@@ -56,28 +56,35 @@ class PathwayMDFData(object):
concentrations = np.exp(l) * standard_concentration
# calculate
dg_prime = PathwayMDFData.CalculateReactionEnergiesUsingConcentrations(
pathway, concentrations)
standard_dg_primes = pathway.standard_dg_primes
physiological_dg_prime = PathwayMDFData.calc_physiological_dg(pathway)
optimized_dg_prime = PathwayMDFData.calc_dg(pathway, concentrations)
stdevs = pathway.sqrt_dg_uncertainty_over_rt @ y
dg_prime_raw = dg_prime + stdevs * R * default_T
# add the calculated error values (due to the dG'0 uncertainty)
optimized_dg_prime += pathway.dg_sigma @ y
# adjust dG to flux directions
dg_prime_adj = (I_dir @ dg_prime_raw) * Q_("kJ/mol")
# adjust dG to flux directions and convert to kJ/mol
standard_dg_primes = (I_dir @ standard_dg_primes) * R * default_T
physiological_dg_prime = (I_dir @ physiological_dg_prime) * R * default_T
optimized_dg_prime = (I_dir @ optimized_dg_prime) * R * default_T
# all dG values are in units of RT
# all dG values are in units of RT, so we convert them to kJ/mol
reaction_data = zip(pathway.reaction_ids(),
pathway.reaction_formulas(),
pathway.fluxes,
pathway.standard_dg_primes,
pathway.physiological_dg_prime,
dg_prime_adj,
standard_dg_primes,
physiological_dg_prime,
optimized_dg_prime,
reaction_prices)
self.reaction_df = pd.DataFrame(
data=list(reaction_data),
columns=['reaction_id', 'reaction_formula', 'flux', 'dG0', 'dGm',
'dGr', 'shadow_price']
columns=['reaction_id',
'reaction_formula',
'flux',
'standard_dg_prime',
'physiological_dg_prime',
'optimized_dg_prime',
'shadow_price']
)
compound_data = zip(pathway.compound_names,
......@@ -86,30 +93,51 @@ class PathwayMDFData(object):
self.compound_df = pd.DataFrame(
data=list(compound_data),
columns=['compound', 'concentration', 'shadow_price']
columns=['compound',
'concentration',
'shadow_price']
)
lbs, ubs = pathway.bounds
self.compound_df['lower_bound'] = list(lbs)
self.compound_df['upper_bound'] = list(ubs)
@staticmethod
def calc_physiological_dg(
pathway: object
) -> np.array:
""" Add the effect of reactant physiological concentrations on the dG'
:param pathway: Pathway object
:return: the reaction energies (in units of RT)
"""
dg_adj = np.array([float(r.physiological_dg_correction()) for
r in pathway.reactions])
return pathway.standard_dg_primes + dg_adj
@staticmethod
@ureg.check(None, "[concentration]")
def CalculateReactionEnergiesUsingConcentrations(
pathway,
def calc_dg(
pathway: object,
concentrations: np.array
) -> np.array:
""" Add the effect of reactant concentrations on the dG'
:param pathway: Pathway object
:param concentrations: a NumPy array of concentrations
:return: the reaction energies (in units of RT)
"""
log_conc = np.log(concentrations / standard_concentration)
if np.isnan(pathway.standard_dg_primes).any():
dg_adj = np.zeros(pathway.S.shape[1]) * Q_("kJ/mol")
dg_adj = np.zeros(pathway.S.shape[1])
for r in range(pathway.S.shape[1]):
reactants = list(pathway.S[:, r].nonzero()[0].flat)
dg_adj[r] = log_conc[reactants] @ pathway.S[reactants, r]
else:
dg_adj = pathway.S.T.values @ log_conc
return pathway.standard_dg_primes + dg_adj * R * default_T
return pathway.standard_dg_primes + dg_adj
@property
def compound_plot(self) -> plt.Figure:
......@@ -162,9 +190,9 @@ class PathwayMDFData(object):
ureg.setup_matplotlib(True)
cumulative_dgms = [Q_("0 kJ/mol")] + np.cumsum(
self.reaction_df.dGm).tolist()
self.reaction_df.physiological_dg_prime).tolist()
cumulative_dgs = [Q_("0 kJ/mol")] + np.cumsum(
self.reaction_df.dGr).tolist()
self.reaction_df.optimized_dg_prime).tolist()
# This is an ugly hack. For some reason, Pandas messes up the numpy
# arrays with pint units, and they need to be rebuilt from scratch to
......@@ -217,7 +245,7 @@ class PathwayThermoModel(object):
# Make sure dG0_r' is the right size
assert pathway.standard_dg_primes.shape == (self.Nr, ), \
'standard dG required for all reactions'
assert pathway.dg_uncertainties.shape == (self.Nr, self.Nr), \
assert pathway.dg_sigma.shape == (self.Nr, self.Nr), \
'uncertainty in dG required for all reactions'
self.I_dir = np.diag(np.sign(self.pathway.fluxes))
......@@ -268,7 +296,7 @@ class PathwayThermoModel(object):
"""
inds = np.nonzero(np.diag(self.I_dir))[0].tolist()
A11 = self.I_dir[inds] @ self.pathway.sqrt_dg_uncertainty_over_rt
A11 = self.I_dir[inds] @ self.pathway.dg_sigma
A12 = self.I_dir[inds] @ self.pathway.S.T
A13 = np.ones((len(inds), 1))
......@@ -283,7 +311,7 @@ class PathwayThermoModel(object):
A33 = np.zeros((self.Nc, 1))
# upper bound values
b1 = -self.I_dir[inds] @ self.pathway.standard_dg_prime_over_rt
b1 = -self.I_dir[inds] @ self.pathway.standard_dg_primes
b2 = np.ones(self.Nr) * self.stdev_factor
A = np.vstack([np.hstack([ A11, A12, A13]), # driving force
......
......@@ -32,32 +32,31 @@ import warnings
import pytest
from equilibrator_api import Q_, Pathway, Reaction
from sbtab import SBtab
from . import approx_unit
TEST_DIR = os.path.dirname(
@pytest.fixture(scope="module")
def toy_pathway() -> Pathway:
test_dir = os.path.dirname(
os.path.abspath(inspect.getfile(inspect.currentframe())))
sbtab = SBtab.read_csv(os.path.join(test_dir, 'pathway.tsv'), "pathway")
return Pathway.from_sbtab(sbtab)
def test_mdf():
def test_shadow_prices(toy_pathway):
warnings.simplefilter('ignore', ResourceWarning)
sbtab_fname = os.path.join(TEST_DIR, 'pathway.tsv')
pp = Pathway.from_sbtab(sbtab_fname)
net_rxn = pp.net_reaction()
net_rxn = toy_pathway.net_reaction()
reference_net_rxn = Reaction.parse_formula(
"2 KEGG:C00008 + 2 KEGG:C00009 + KEGG:C00031 = "
"2 KEGG:C00002 + 2 KEGG:C00001 + 2 KEGG:C00469 + 2 KEGG:C00011")
assert net_rxn == reference_net_rxn
pp_mdf_data = pp.calc_mdf(stdev_factor=0.0)
approx_unit(pp_mdf_data.mdf, Q_("2.04 kJ/mol"), abs=0.1)
mdf_data = toy_pathway.calc_mdf(stdev_factor=1.0)
pp_mdf_data = pp.calc_mdf(stdev_factor=1.0)
approx_unit(pp_mdf_data.mdf, Q_("2.75 kJ/mol"), abs=0.1)
shadow_prices = pp_mdf_data.reaction_df.set_index(
shadow_prices = mdf_data.reaction_df.set_index(
'reaction_id').shadow_price
assert shadow_prices['glucokinase'] == pytest.approx(0.0, abs=1e-2)
......@@ -65,7 +64,7 @@ def test_mdf():
assert shadow_prices['tim'] == pytest.approx(0.25, abs=1e-2)
assert shadow_prices['gapdh'] == pytest.approx(0.50, abs=1e-2)
compound_prices = pp_mdf_data.compound_df.set_index('compound').shadow_price
compound_prices = mdf_data.compound_df.set_index('compound').shadow_price
assert compound_prices['h2o'] == pytest.approx(0.0, abs=1e-2)
assert compound_prices['nad'] == pytest.approx(0.5, abs=1e-2)
......@@ -73,7 +72,23 @@ def test_mdf():
# test that the plotting functions run without crashing
try:
pp_mdf_data.reaction_plot
pp_mdf_data.compound_plot
mdf_data.reaction_plot
mdf_data.compound_plot
except:
pytest.fail("MDF plot function raises exception")
@pytest.mark.parametrize(
"p_h, ionic_strength, stdev_factor, expected_mdf",
[
(Q_("6"), Q_("0.25M"), 1.0, Q_("-2.08 kJ/mol")),
(Q_("7"), Q_("0.1M"), 1.0, Q_("1.88 kJ/mol")),
(Q_("8"), Q_("0.1M"), 1.0, Q_("4.77 kJ/mol")),
(Q_("7"), Q_("0.25M"), 0.0, Q_("2.04 kJ/mol")),
(Q_("7"), Q_("0.25M"), 1.0, Q_("2.75 kJ/mol")),
(Q_("7"), Q_("0.25M"), 2.0, Q_("3.46 kJ/mol")),
])
def test_mdf(p_h, ionic_strength, stdev_factor, expected_mdf, toy_pathway):
toy_pathway.set_aqueous_params(p_h=p_h, ionic_strength=ionic_strength)
mdf_data = toy_pathway.calc_mdf(stdev_factor=stdev_factor)
approx_unit(mdf_data.mdf, expected_mdf, abs=0.1)
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