Commit 4952dbf7 authored by eladnoor's avatar eladnoor

implementing the uncertainty matrix constraints in MDF analysis

parent d7ebdaba
......@@ -31,6 +31,7 @@ from io import IOBase
from typing import Callable, Iterable, List, TextIO, Tuple
import numpy as np
from scipy.linalg import fractional_matrix_power
from equilibrator_cache.reaction import (
create_stoichiometric_matrix_from_reactions)
from sbtab import SBtab
......@@ -97,7 +98,10 @@ class Pathway(object):
# TODO: right now MDF does not work with uncertainties, and we
# need to fix that. In the meantime, we use a 0-matrix to
# indicate there is no uncertainty in the dG'0 estimates.
self.dg_uncertainties *= 0.0
# for i in range(Nr):
# for j in range(Nr):
# if i != j:
# self.dg_uncertainties[i, j] *= 0.0
else:
self.standard_dg_primes = standard_dg_primes
self.dg_uncertainties = dg_uncertainties
......@@ -189,10 +193,14 @@ class Pathway(object):
@property
def sqrt_dg_uncertainty_over_rt(self) -> np.array:
u_over_rt = self.dg_uncertainties**(0.5) / (R * default_T)
return np.array(list(map(float, u_over_rt.flat))).reshape(
# get rid of the kJ/mol units, by dividing the uncertainties by RT^2
u_over_rt = self.dg_uncertainties / (R * default_T)**2
u_over_rt = np.array(list(map(float, u_over_rt.flat))).reshape(
u_over_rt.shape)
# now take the square root of the matrix
return fractional_matrix_power(u_over_rt, 0.5)
@classmethod
def from_csv_file(cls,
f: TextIO,
......@@ -241,8 +249,8 @@ class Pathway(object):
return sorted(compounds)
def calc_mdf(self) -> object:
return PathwayThermoModel(self).FindMDF()
def calc_mdf(self, stdev_factor: float = 1.0) -> object:
return PathwayThermoModel(self, stdev_factor).FindMDF()
def reaction_ids(self) -> Iterable[str]:
return map(lambda rxn: rxn.reaction_id, self.reactions)
......
......@@ -40,28 +40,28 @@ from .settings import (
class PathwayMDFData(object):
@ureg.check(None, None, "kJ/mol", None, "[concentration]", None, None, None)
def __init__(
self,
pathway,
mdf: float,
B: float,
I_dir: np.array,
concentrations: np.array,
dg_cov_eigen: np.array,
l: np.array,
y: np.array,
reaction_prices: np.array,
compound_prices: np.array
) -> object:
self.mdf = mdf
self.dg_cov_eigen = dg_cov_eigen
self.mdf = B * R * default_T
self.bounds = pathway._bounds.Copy()
concentrations = np.exp(l) * standard_concentration
# calculate
dg_prime = PathwayMDFData.CalculateReactionEnergiesUsingConcentrations(
pathway, concentrations)
dg_prime_raw = (dg_prime +
dg_cov_eigen.T @ pathway.dg_uncertainties)
stdevs = pathway.sqrt_dg_uncertainty_over_rt @ y
dg_prime_raw = dg_prime + stdevs * R * default_T
# adjust dG to flux directions
dg_prime_adj = (I_dir @ dg_prime_raw) * Q_("kJ/mol")
......@@ -174,8 +174,8 @@ class PathwayMDFData(object):
cumulative_dgs = np.array([float(x/Q_('kJ/mol')) for x in
cumulative_dgs]) * Q_('kJ/mol')
xticks = np.arange(0, self.reaction_df.shape[0] + 1)
xticklabels = [''] + self.reaction_df.reaction_id.tolist()
xticks = 0.5 + np.arange(self.reaction_df.shape[0])
xticklabels = self.reaction_df.reaction_id.tolist()
reaction_fig, ax = plt.subplots(1, 1, figsize=(9, 6))
ax.xaxis.set_units(Q_("kJ/mol"))
......@@ -199,8 +199,19 @@ class PathwayMDFData(object):
class PathwayThermoModel(object):
"""Container for doing pathway-level thermodynamic analysis."""
def __init__(self, pathway) -> object:
def __init__(
self,
pathway,
stdev_factor: float = 1.0
) -> object:
"""Create a model for running MDF analysis.
:param pathway: a Pathway object defining the reactions and bounds
:param stdev_factor: the factor with which we multiply the
sqrt uncertainty matrix of the dG'0 estimates
"""
self.pathway = pathway
self.stdev_factor = stdev_factor
self.Nc, self.Nr = pathway.S.shape
# Make sure dG0_r' is the right size
......@@ -242,8 +253,7 @@ class PathwayThermoModel(object):
def _MakeDrivingForceConstraints(
self
) -> Tuple[np.array, np.array, np.array]:
"""
Generates the A matrix and b & c vectors that can be used in a
"""Generates the A matrix and b & c vectors that can be used in a
standard form linear problem:
max c'x
subject to Ax <= b
......@@ -253,8 +263,9 @@ class PathwayThermoModel(object):
are the natural log of the concentrations of metabolites, and
B is the max-min driving force variable which is being maximized
by the LP
:return: (A, b, c) - the parameters of the LP
"""
# Define and apply the constraints on the concentrations
inds = np.nonzero(np.diag(self.I_dir))[0].tolist()
A11 = self.I_dir[inds] @ self.pathway.sqrt_dg_uncertainty_over_rt
......@@ -273,7 +284,7 @@ class PathwayThermoModel(object):
# upper bound values
b1 = -self.I_dir[inds] @ self.pathway.standard_dg_prime_over_rt
b2 = np.ones(self.Nr)
b2 = np.ones(self.Nr) * self.stdev_factor
A = np.vstack([np.hstack([ A11, A12, A13]), # driving force
np.hstack([ A21, A22, A23]), # covariance var ub
......@@ -413,7 +424,7 @@ class PathwayThermoModel(object):
A 2 (optimal dGfs, optimal concentrations, optimal mdf).
"""
def get_primal_array(l):
return np.array([v.primal for v in l])
return np.array([v.primal for v in l], ndmin=1)
lp_primal, y, l, B = self._MakeMDFProblem()
......@@ -421,10 +432,9 @@ class PathwayThermoModel(object):
logging.warning('LP status %s', lp_primal.status)
raise Exception("Cannot solve MDF primal optimization problem")
y = get_primal_array(y)
l = get_primal_array(l)
mdf = lp_primal.variables['mdf'].primal * R * default_T
conc = np.exp(l) * standard_concentration
y = get_primal_array(y) # covariance eigenvalue prefactors
l = get_primal_array(l) # log concentrations
B = lp_primal.variables['mdf'].primal
lp_dual, w, g, z, u = self._MakeMDFProblemDual()
if lp_dual.optimize() != 'optimal':
......@@ -444,9 +454,9 @@ class PathwayThermoModel(object):
return PathwayMDFData(
self.pathway,
mdf,
B,
self.I_dir,
conc,
l,
y,
reaction_prices,
compound_prices
......
......@@ -44,8 +44,6 @@ def test_mdf():
sbtab_fname = os.path.join(TEST_DIR, 'pathway.tsv')
pp = Pathway.from_sbtab(sbtab_fname)
pp_mdf_data = pp.calc_mdf()
net_rxn = pp.net_reaction()
reference_net_rxn = Reaction.parse_formula(
......@@ -53,8 +51,12 @@ def test_mdf():
"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)
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(
'reaction_id').shadow_price
......
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