Commit 5b4474d4 authored by eladnoor's avatar eladnoor

fixing command line script for calculating dG0s in iJO1336

parent 5b666cff
Pipeline #58235483 (#146) failed with stage
in 3 minutes and 51 seconds
......@@ -29,61 +29,105 @@ import argparse
import csv
import logging
import sys
from numpy import nan, sqrt
from equilibrator_api import ComponentContribution, Reaction
def main(args):
from equilibrator_api import (
Q_, ComponentContribution, Reaction) # isort:skip
p_h = Q_(args.ph)
assert p_h.check(None)
if __name__ == '__main__':
ionic_strength = Q_(args.i)
if ionic_strength.check(None):
ionic_strength *= Q_("M")
else:
assert ionic_strength.check("[concentration]")
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('--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)
temperature = Q_(args.t)
if temperature.check(None):
temperature *= Q_("K")
else:
assert temperature.check("[temperature]")
args = parser.parse_args()
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("Initializing the Component Contribution estimator.\n")
eq_api = ComponentContribution(p_h=p_h, ionic_strength=ionic_strength,
temperature=temperature)
sys.stderr.write('pH = %.1f\n' % args.ph)
sys.stderr.write('I = %.1f M\n' % args.i)
sys.stderr.write("Parsing model reactions.\n")
ids = []
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):
ids.append(row['bigg.reaction'])
all_ids.append(row['bigg.reaction'])
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:
print('warning: cannot parse reaction %s because of %s' %
(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)
header = ['reaction', 'pH', 'ionic strength [M]', 'dG\'0 [kJ/mol]',
'uncertainty [kJ/mol]', 'ln(Reversibility Index)', 'comment']
header = ["reaction", "pH", "ionic strength [M]", "ΔG'0",
"uncertainty", "ln(Reversibility Index)", "comment"]
writer.writerow(header)
for s, r, dg0, u in zip(ids, reactions,
dG0_prime.flat, U.diagonal().flat):
row = [s, args.ph, args.i]
if r.is_empty():
row += [nan, nan, nan, 'reaction is empty']
elif r.check_full_reaction_balancing():
ln_RI = cc.reversibility_index(r)
row += ['%.2f' % dg0, '%.2f' % sqrt(u), '%.2f' % ln_RI, '']
else:
row += [nan, nan, nan, 'reaction is not chemically balanced']
writer.writerow(row)
for rid, comment in zip(all_ids, comments):
dg0 = nan
sigma = nan
ln_RI = nan
if rid in balanced_ids:
sigma = sqrt(rid_to_uncertainty[rid])
if sigma < Q_(100, "kJ/mol"):
dg0 = rid_to_dg0_prime[rid]
ln_RI = eq_api.ln_reversibility_index(rid_to_reaction[rid])
writer.writerow([rid, args.ph, args.i, dg0, sigma, ln_RI, comment])
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.
This diff is collapsed.
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