Commit 0635e72d authored by Elad Noor's avatar Elad Noor

fixing bug in mdf_plot

parent a7f5c1b9
......@@ -17,7 +17,7 @@ if __name__ == '__main__':
parser = argparse.ArgumentParser(
description='Calculate the Max-min Driving Force (MDF) of a pathway.')
parser.add_argument(
'infile', type=argparse.FileType(),
'infile', type=str,
help='path to input file containing reactions')
parser.add_argument(
'outfile', type=str,
......@@ -27,13 +27,13 @@ if __name__ == '__main__':
args = parser.parse_args()
pp = Pathway.from_sbtab(args.infile)
output_pdf = PdfPages(args.outfile)
mdf_res = pp.calc_mdf()
output_pdf.savefig(mdf_res.conc_plot)
output_pdf.savefig(mdf_res.mdf_plot)
output_pdf.close()
rxn_df = pd.DataFrame(mdf_res.report_reactions)
cpd_df = pd.DataFrame(mdf_res.report_compounds)
......@@ -35,6 +35,7 @@ class ReactionMDFData(object):
'dGm': self.dGm, 'dG0': self.dG0,
'shadow_price': self.shadow_price}
class CompoundMDFData(object):
def __init__(self, compound, concentration_bounds,
concentration, shadow_price):
......@@ -55,13 +56,14 @@ class CompoundMDFData(object):
'ub': self.ub,
'shadow_price': self.shadow_price}
class PathwayMDFData(object):
def __init__(self, parsed_pathway, mdf_result):
self.parsed_pathway = parsed_pathway
self.mdf_result = mdf_result
self.model = mdf_result.model
self.reaction_data = []
for i in range(len(self.parsed_pathway.reactions)):
mdf_data = ReactionMDFData(self.parsed_pathway.reactions[i],
......@@ -71,7 +73,7 @@ class PathwayMDFData(object):
self.mdf_result.dG_r_prime_adj[i],
self.mdf_result.reaction_prices[i])
self.reaction_data.append(mdf_data)
compound_ids = parsed_pathway.compound_kegg_ids
compound_names = parsed_pathway.compound_names
cbounds = [self.model.concentration_bounds.GetBoundTuple(cid)
......@@ -79,7 +81,8 @@ class PathwayMDFData(object):
concs = self.mdf_result.concentrations
prices = self.mdf_result.compound_prices
self.compound_data = [CompoundMDFData(*t)
for t in zip(compound_names, cbounds, concs, prices)]
for t in zip(compound_names, cbounds,
concs, prices)]
@property
def mdf(self):
......@@ -106,7 +109,8 @@ class PathwayMDFData(object):
import matplotlib.pyplot as plt
ys = np.arange(0, len(self.compound_data))
concs = np.array([c.concentration for c in self.compound_data], ndmin=2)
concs = np.array([c.concentration for c in self.compound_data],
ndmin=2)
cnames = [str(c.compound) for c in self.compound_data]
default_lb = self.model.concentration_bounds.default_lb
default_ub = self.model.concentration_bounds.default_ub
......@@ -122,13 +126,13 @@ class PathwayMDFData(object):
concs_equal = concs[bounds_equal]
# Special color for metabolites with nonzero shadow prices.
shadow_prices = np.array([c.shadow_price for c in self.compound_data], ndmin=2)
shadow_prices = np.array([c.shadow_price for c in self.compound_data],
ndmin=2)
nz_shadow = np.where(shadow_prices != 0)[0]
ys_nz_shadow = ys[nz_shadow]
concs_nz_shadow = concs[nz_shadow]
conc_figure, ax = plt.subplots(1, 1, figsize=(9, 6))
#ax = plt.axes([0.2, 0.1, 0.7, 0.8])
ax.axvspan(1e-8, default_lb, color='y', alpha=0.5)
ax.axvspan(default_ub, 1e3, color='y', alpha=0.5)
ax.scatter(concs, ys,
......@@ -145,10 +149,10 @@ class PathwayMDFData(object):
ax.set_xlim(default_lb*0.1, 1.5)
ax.set_ylim(-1.5, len(self.compound_data) + 0.5)
conc_figure.tight_layout()
return conc_figure
@property
def mdf_plot(self):
import matplotlib.pyplot as plt
......@@ -160,7 +164,8 @@ class PathwayMDFData(object):
cumulative_dgs = np.cumsum(dgs)
xticks = np.arange(0, len(cumulative_dgs))
xticklabels = [''] + [r.reaction.reaction_id for r in self.reaction_data]
xticklabels = [''] + \
[r.reaction.reaction_id for r in self.reaction_data]
mdf_fig, ax = plt.subplots(1, 1, figsize=(9, 6))
ax.plot(cumulative_dgms, label='Physiological concentrations (1 mM)',
......@@ -172,7 +177,7 @@ class PathwayMDFData(object):
ax.grid(color='grey', linestyle='--', linewidth=1, alpha=0.2)
ax.set_xlabel('Reaction Step')
ax.set_ylabel("Cumulative $\Delta_r G'$ (kJ/mol)")
ax.set_ylabel(r'Cumulative $\Delta_r G^\prime$ (kJ/mol)')
ax.legend(loc='best')
ax.set_title('MDF = %.1f' % self.mdf)
......
......@@ -54,7 +54,8 @@ class Pathway(object):
bounds: bounds on metabolite concentrations.
Uses default bounds if None provided.
pH: (optional) specify the pH at which the dG values are calculated
ionic_strength: (optional) specify the I at which the dG values are calculated
ionic_strength: (optional) specify the I at which the dG values
are calculated
"""
assert len(reactions) == len(fluxes)
assert len(reactions) == len(dG0_r_primes)
......@@ -75,7 +76,7 @@ class Pathway(object):
self.compound_kegg_ids))
else:
self.compound_names = list(self.compound_kegg_ids)
nr, nc = self.S.shape
# dGr should be orthogonal to nullspace of S
......@@ -114,8 +115,10 @@ class Pathway(object):
reactions.append(rxn)
fluxes.append(flux)
equilibrator = ComponentContribution(pH=pH, ionic_strength=ionic_strength)
dG0_r_primes, dG0_uncertainties = zip(*map(equilibrator.dG0_prime, reactions))
equilibrator = ComponentContribution(pH=pH,
ionic_strength=ionic_strength)
dG0_r_primes, dG0_uncertainties = zip(*map(equilibrator.dG0_prime,
reactions))
dG0_r_primes = list(dG0_r_primes)
return Pathway(reactions, fluxes, dG0_r_primes, bounds=bounds)
......@@ -132,14 +135,15 @@ class Pathway(object):
"""Builds a stoichiometric matrix.
Returns:
Two tuple (S, compounds) where compounds is the KEGG IDs of the compounds
in the order defining the column order of the stoichiometric matrix S.
Two tuple (S, compounds) where compounds is the KEGG IDs of
the compounds in the order defining the column order of the
stoichiometric matrix S.
"""
compounds = set()
for r in self.reactions:
compounds.update(r.kegg_ids)
compounds = sorted(compounds)
smat = zeros((len(compounds), len(self.reactions)))
for j, r in enumerate(self.reactions):
for i, c in enumerate(compounds):
......@@ -170,20 +174,20 @@ class Pathway(object):
table_ids = ['ConcentrationConstraint', 'Reaction', 'RelativeFlux',
'Parameter']
dfs = []
for table_id in table_ids:
sbtab = sbtabdoc.get_sbtab_by_id(table_id)
if sbtab is None:
logging.error('%s contains the following TableIDs: %s' %
(file_name,
', '.join(map(lambda s: s.table_id, sbtabdoc.sbtabs))))
tables = ', '.join(map(lambda s: s.table_id, sbtabdoc.sbtabs))
logging.error('%s contains the following TableIDs: %s' %
(file_name, tables))
raise PathwayParseError('The SBtabDocument must have a table '
'with the following TableID: %s'
% table_id)
dfs.append(sbtab.to_data_frame())
bounds_df, reaction_df, flux_df, keqs_df = dfs
bounds_unit = sbtabdoc.get_sbtab_by_id(
'ConcentrationConstraint').get_attribute('Unit')
bounds = Bounds.from_dataframe(bounds_df, bounds_unit)
......
......@@ -109,20 +109,29 @@ def test_mdf():
sbtab_fname = os.path.join(TEST_DIR, 'pathway_ethanol_SBtab.tsv')
pp = Pathway.from_sbtab(sbtab_fname)
mdf_res = pp.calc_mdf().mdf_result
pp_mdf_data = pp.calc_mdf()
assert mdf_res.mdf == pytest.approx(1.69, abs=1e-2)
assert mdf_res.max_total_dG == pytest.approx(-159.05, abs=1e-2)
assert mdf_res.min_total_dG == pytest.approx(-181.05, abs=1e-2)
assert pp_mdf_data.mdf == pytest.approx(1.69, abs=1e-2)
assert pp_mdf_data.max_total_dG == pytest.approx(-159.05, abs=1e-2)
assert pp_mdf_data.min_total_dG == pytest.approx(-181.05, abs=1e-2)
assert mdf_res.reaction_prices[0, 0] == pytest.approx(0.0, abs=1e-2)
assert mdf_res.reaction_prices[3, 0] == pytest.approx(0.25, abs=1e-2)
assert mdf_res.reaction_prices[4, 0] == pytest.approx(0.25, abs=1e-2)
assert mdf_res.reaction_prices[5, 0] == pytest.approx(0.50, abs=1e-2)
shadow_prices = list(map(lambda r: r.shadow_price,
pp_mdf_data.reaction_data))
assert mdf_res.compound_prices[1, 0] == pytest.approx(0.0, abs=1e-2)
assert mdf_res.compound_prices[2, 0] == pytest.approx(1.24, abs=1e-2)
assert mdf_res.compound_prices[3, 0] == pytest.approx(-1.24, abs=1e-2)
assert shadow_prices[0] == pytest.approx(0.0, abs=1e-2)
assert shadow_prices[3] == pytest.approx(0.25, abs=1e-2)
assert shadow_prices[4] == pytest.approx(0.25, abs=1e-2)
assert shadow_prices[5] == pytest.approx(0.50, abs=1e-2)
compound_prices = list(map(lambda r: r.shadow_price,
pp_mdf_data.compound_data))
assert compound_prices[1] == pytest.approx(0.0, abs=1e-2)
assert compound_prices[2] == pytest.approx(1.24, abs=1e-2)
assert compound_prices[3] == pytest.approx(-1.24, abs=1e-2)
# test the mdf_plot() runs without crashing
pp_mdf_data.mdf_plot
def test_reaction_matcher():
......
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