diff --git a/corda/corda.py b/corda/corda.py index 285636a..1380910 100644 --- a/corda/corda.py +++ b/corda/corda.py @@ -17,6 +17,7 @@ TOL = 1e-6 # Tolerance to judge whether a flux is non-zero UPPER = 1e6 # default upper bound +UN = (1e-6, 1e-4) # uniform noise class CORDA(object): @@ -119,7 +120,6 @@ def __init__(self, model, confidence, met_prod=None, n=5, self.tflux = 1 self.impossible = [] self.n = n - self.noise = 0.1 self.support = support self.pf = penalty_factor self.solver = solver_dict[get_solver_name() if solver is None @@ -127,10 +127,10 @@ def __init__(self, model, confidence, met_prod=None, n=5, self.sargs = solver_kwargs def __perturb(self, lp, pen): - noise = np.random.uniform(high=self.noise, size=len(pen)) + noise = np.random.uniform(low=UN[0], high=UN[1], size=len(pen)) for i, p in enumerate(pen): - if p > 0.0: + if p < 2.0: self.solver.change_variable_objective(lp, i, p + noise[i]) def __quiet_solve(self, lp, os): @@ -154,8 +154,8 @@ def __reduce_conf(self, conf): red_conf = dict.fromkeys(rids, -1) for k, v in conf.items(): - k = k.replace("_reverse", "") - red_conf[k] = max(red_conf[k], v) + kr = k.replace("_reverse", "") + red_conf[kr] = max(red_conf[kr], v) return red_conf def associated(self, targets, conf=None, penalize_medium=True): @@ -309,7 +309,7 @@ def info(self, reversible=True): " - high: {}\n".format(old_counts[3]) else: old = np.array([conf_old[k] for k in conf_old]) - new = np.array([conf[k] for k in conf]) + new = np.array([conf[k] for k in conf_old]) med_inc = np.sum(((old == 1) | (old == 2)) & (new == 3)) noc_inc = np.sum((old == -1) & (new == 3)) free_inc = np.sum((old == 0) & (new == 3)) diff --git a/docs/Makefile b/docs/Makefile index 30e95ff..e715749 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -1,2 +1,2 @@ all: - jupyter-nbconvert *.ipynb + jupyter-nbconvert --execute *.ipynb diff --git a/docs/index.html b/docs/index.html index d6d02b0..cae981b 100644 --- a/docs/index.html +++ b/docs/index.html @@ -11951,7 +11951,7 @@

A small example
from corda import CORDA
 
-opt = CORDA(mod, conf)
+opt = CORDA(mod, conf, n=20)
 opt.build()
 print(opt)
 
@@ -11967,9 +11967,9 @@

A small example
@@ -12010,9 +12010,9 @@

A small example
+ +
+
+
In [7]:
+
+
+
print([opt.model.reactions.get_by_id(k).reaction for k, v in opt.conf.items() if v == 3])
+
+ +
+
+
+ +
+
+ + +
+
+
['adp + pi --> atp', ' --> lac_L', 'ru5p_D --> r5p', 'amp --> ', 'lac_L + nad --> nadh + pyr', ' --> adp', 'e4p + f6p --> g3p + s7p', 'f6p + g3p --> e4p + xu5p_D', 'atp + r5p --> amp + prpp', 'g3p + nad + pi --> 13dpg + nadh', 'gln --> glu', 'atp + pyr --> adp + oaa + pi', ' --> g1p', ' --> gln', 'atp + f6p --> adp + fdp', ' --> pi', 'g6p --> f6p', '13dpg --> 3pg + pi', 'xu5p_D --> ru5p_D', '0.39253 3pg + 20.7045 atp + 0.15446 cit + 0.38587 glu + 0.35261 oaa + 0.053446 prpp + 0.50563 pyr --> 20.6508 adp + 20.6508 pi', 'nadh --> nad', 'g3p + s7p --> r5p + xu5p_D', 'icit --> cit', 'fdp --> dhap + g3p', 'dhap --> g3p', 'g1p --> g6p', 'glu + nad --> akg + nadh', 'akg + nadh --> icit + nad']
+
+
+
+ +
+
+ +
+
+
+
+
+
+

This gives a reconstruction that looks like this (orange denotes included reactions): reconstruction

@@ -12047,7 +12084,7 @@

A small example
-
In [7]:
+
In [8]:
opt = CORDA(mod, conf, met_prod="pep")
@@ -12066,9 +12103,9 @@ 

A small example
build status: reconstruction complete
-Inc. reactions: 38/61
+Inc. reactions: 31/61
  - unclear: 0/0
- - exclude: 36/59
+ - exclude: 29/59
  - low and medium: 0/0
  - high: 2/2
 
@@ -12092,11 +12129,10 @@ 

A small example diff --git a/docs/index.ipynb b/docs/index.ipynb index 153a96d..96f666a 100644 --- a/docs/index.ipynb +++ b/docs/index.ipynb @@ -37,7 +37,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 1, "metadata": { "collapsed": false }, @@ -48,7 +48,7 @@ "1" ] }, - "execution_count": 28, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" } @@ -78,7 +78,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 2, "metadata": { "collapsed": false }, @@ -89,7 +89,7 @@ "60" ] }, - "execution_count": 29, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } @@ -112,7 +112,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 3, "metadata": { "collapsed": false }, @@ -123,7 +123,7 @@ "'0.39253 3pg + 20.7045 atp + 0.15446 cit + 0.38587 glu + 0.35261 oaa + 0.053446 prpp + 0.50563 pyr --> 20.6508 adp + 20.6508 pi'" ] }, - "execution_count": 30, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -141,7 +141,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 4, "metadata": { "collapsed": true }, @@ -161,7 +161,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -171,9 +171,9 @@ "output_type": "stream", "text": [ "build status: reconstruction complete\n", - "Inc. reactions: 33/60\n", + "Inc. reactions: 28/60\n", " - unclear: 0/0\n", - " - exclude: 32/59\n", + " - exclude: 27/59\n", " - low and medium: 0/0\n", " - high: 1/1\n", "\n" @@ -183,7 +183,7 @@ "source": [ "from corda import CORDA\n", "\n", - "opt = CORDA(mod, conf)\n", + "opt = CORDA(mod, conf, n=20)\n", "opt.build()\n", "print(opt)" ] @@ -192,12 +192,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The metric you see are for reversible reactions. We can obtain the irreversible reconstruction metrics by using:" + "The metric you sare for reversible reactions. We can obtain the irreversible reconstruction metrics by using:" ] }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 6, "metadata": { "collapsed": false }, @@ -207,9 +207,9 @@ "output_type": "stream", "text": [ "build status: reconstruction complete\n", - "Inc. reactions: 33/101\n", + "Inc. reactions: 28/101\n", " - unclear: 0/0\n", - " - exclude: 32/100\n", + " - exclude: 27/100\n", " - low and medium: 0/0\n", " - high: 1/1\n", "\n" @@ -224,7 +224,34 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "This gives a reconstruction that looks like this (red denotes included reactions):\n", + "In order to see the list of reactions in the reconstructed model without specifically creating the new model you can use the following snippet." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[' --> adp', 'atp + r5p --> amp + prpp', 'lac_L + nad --> nadh + pyr', 'nadh --> nad', ' --> pi', '13dpg --> 3pg + pi', 'g1p --> g6p', 'f6p + g3p --> e4p + xu5p_D', '0.39253 3pg + 20.7045 atp + 0.15446 cit + 0.38587 glu + 0.35261 oaa + 0.053446 prpp + 0.50563 pyr --> 20.6508 adp + 20.6508 pi', 'icit --> cit', 'akg + nadh --> icit + nad', 'g6p --> f6p', 'adp + pi --> atp', ' --> gln', 'glu + nad --> akg + nadh', ' --> lac_L', 'g3p + nad + pi --> 13dpg + nadh', 'ru5p_D --> r5p', 'xu5p_D --> ru5p_D', 'atp + pyr --> adp + oaa + pi', 'atp + f6p --> adp + fdp', 'gln --> glu', 'amp --> ', 'dhap --> g3p', 'e4p + f6p --> g3p + s7p', 'fdp --> dhap + g3p', ' --> g1p', 'g3p + s7p --> r5p + xu5p_D']\n" + ] + } + ], + "source": [ + "print([opt.model.reactions.get_by_id(k).reaction for k, v in opt.conf.items() if v == 3])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This gives a reconstruction that looks like this (orange denotes included reactions):\n", "![reconstruction](reconstruction.png)" ] }, @@ -237,7 +264,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -247,9 +274,9 @@ "output_type": "stream", "text": [ "build status: reconstruction complete\n", - "Inc. reactions: 39/61\n", + "Inc. reactions: 31/61\n", " - unclear: 0/0\n", - " - exclude: 37/59\n", + " - exclude: 29/59\n", " - low and medium: 0/0\n", " - high: 2/2\n", "\n" @@ -271,7 +298,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -280,15 +307,12 @@ "name": "stdout", "output_type": "stream", "text": [ - "38\n", - "gtp + oaa <=> gdp + pep\n", "2pg <=> pep\n" ] } ], "source": [ "rec = opt.cobra_model(\"plus_pep\")\n", - "print(len(rec.reactions))\n", "use = rec.metabolites.pep.reactions\n", "for r in use: print(r.reaction)" ] @@ -297,12 +321,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "By default CORDA uses redundancy. This means, in case there are several minimal pathways to reach your objective, CORDA will include several of those (which is good since it gives your model some robustness). If we do not want that feature we can modify the parameter n in the CORDA initializer which denotes the maximum number of redundant pathways to include." + "This is the unique solution for that additional constraint because it only includes 2 new reactions (3pg --> 2pg --> pep). The alternative would be using oxaloacetate, however that would require an additional import og GTP and a sink for GDP and would thus be more expensive (3 reactions).\n", + "\n", + "By default CORDA uses redundancy. This means, in case there are several minimal pathways to reach your objective, CORDA will include several of those (which is good since it gives your model some robustness). In this case this did not happen since the best solution was unique. If we do not want that feature we can modify the parameter n in the CORDA initializer which denotes the maximum number of redundant pathways to include." ] }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 10, "metadata": { "collapsed": false, "scrolled": true @@ -312,8 +338,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "33\n", - "gtp + oaa <=> gdp + pep\n", + "30\n", "2pg <=> pep\n" ] } @@ -332,8 +357,17 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "As we can see we can now produce pep via oaa alone. However, the model will now be less robust to deletions." + "Which gives the same solution in that case." ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/docs/model.svg b/docs/model.svg new file mode 100644 index 0000000..fdc9c63 --- /dev/null +++ b/docs/model.svg @@ -0,0 +1,10368 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + ATP + + + Pi + + + NADP + + + NADPH + + + AMP + + + GTP + + GDP + + + + g3p + + + + + + + + + + + + + + + ADP + + + ADP + + + ADP + + + + + + + + lac-L + + + + CO2 + + + + + + + + + + H+ + + + + + H+ + + + + ATP + + + ATP + + + ATP + + + ATP + + + Pi + + + Pi + + + Pi + + + Pi + + + NAD + + + NADP + + + NADH + + + AMP + + + 13dpg + + + 3pg + + + + + + + 23dpg + + + + + + + + + + + + g3p + + + + + + + + + + + ADP + + + ADP + + + ADP + + + + + + 2pg + + + fdp + + + dhap + + + f6p + + + + + pep + + + + g6p + + + + g1p + + + + + + pyr + + + + H+ + + + + NADH + + + + NAD + + + lac-L + + + + H+ + + + + + NADPH + + + NADP + + + 6pgl + + + + + + + 6pgc + + + + + H+ + + + + + NADPH + + + ru5p-D + + + + r1p + + + + r5p + + + + + + + + prpp + + + + + xu5p-D + + + + f6p + + + g3p + + + + + + e4p + + + s7p + + + + + + e4p + + + + + + + + + NAD + + + + + CO2 + + + + H+ + + + + + + NADH + + + oaa + + + cit + + + icit + + + + akg + + + + + A/GDP + + + Pi + + + NAD + + + + + + + + + A/GTP + + + + + H+ + + + + NADH + + + + + + + + + succ + + + + + + + + + + + + fum + + + + mal-L + + + + + + TCA cycle Glycolysis/Gluconeogenesis Pentose phosphate + + + + + + + + + + + + + + + + + gsh + + + gssg + + + + + H+ + + + + + NADPH + + + + NADP + oxidative stress + + + + + + H202 + + + + H20 + + + + + + ATP + + + ADP + + + Pi + + + exchange reactions + + g1p + + + lac-L + + + + + + X + + + + + gln + + + + + + H+ + + + + NAD(P)H + + + NAD(P) + + + + + + glu + + + + ATP + + + + + H+ + + + + ATP + + + Pi + + + ADP + + + + + + + + H+ + + + + NAD + + + NADH + + + O2 + + + + H+ + + + + + + + + + H+ + + oxidativephosphorylation + + + + H+ + + + + + NAD(P)H + + + + NAD(P) + + + + + + + + + NAD(P) + + + + + + H+ + + + + + NAD(P)H + + + + pep + + + GTP + + + GDP + + + + + + + + + + + + + + + + + + + + + ATP + + + cit + proliferation + + + prpp + + + + glu + + + + oaa + + + 3pg + + + pyr + opt. + diff --git a/docs/reconstruction.png b/docs/reconstruction.png index 6632def..a81d7fe 100644 Binary files a/docs/reconstruction.png and b/docs/reconstruction.png differ diff --git a/setup.py b/setup.py index a9fdb04..21d80ba 100644 --- a/setup.py +++ b/setup.py @@ -13,7 +13,7 @@ setup( name='corda', - version='0.1.4', + version='0.2.0', description='Genome-scale model construction with CORDA', long_description=long_description, url='https://github.com/cdiener/corda', diff --git a/tests/tests.py b/tests/tests.py index 082d302..016bf26 100644 --- a/tests/tests.py +++ b/tests/tests.py @@ -8,14 +8,14 @@ from corda import * from cobra import Model, Reaction, Metabolite from cobra.manipulation import convert_to_irreversible, revert_to_reversible -from cobra.io import read_sbml_model + class TestConf(unittest.TestCase): def test_confidence(self): vals = {"g1": -1, "g2": 1, "g3": 2, "g4": 3} cases = [("g1 and g2 or g3", 2), ("g1 and (g2 or g3)", -1), - ("g1 or g2 or g4 or g5", 3), ("g3 and g6", 0), ("", 0)] + ("g1 or g2 or g4 or g5", 3), ("g3 and g6", 0), ("", 0)] for rule, res in cases: conf = reaction_confidence(rule, vals) @@ -24,11 +24,13 @@ def test_confidence(self): def test_eval_safe(self): cases = ["print()", "A + B", "A ^ B"] for ca in cases: - with self.assertRaises(TypeError): reaction_confidence(ca, {}) + with self.assertRaises(TypeError): + reaction_confidence(ca, {}) def test_none(self): self.assertEqual(reaction_confidence(" ", {}), 0) + class TestRevert(unittest.TestCase): def setUp(self): @@ -45,12 +47,14 @@ def setUp(self): def test_remove_breaks(self): self.assertRaises(KeyError, revert_to_reversible, self.model) + class TestExamples(unittest.TestCase): def test_cemet(self): model = test_model() self.assertEqual(len(model.reactions), 60) self.assertEqual(len(model.metabolites), 43) + class TestCORDAsimple(unittest.TestCase): def setUp(self): A = Metabolite("A") @@ -110,7 +114,8 @@ def test_impossible_req(self): def test_association_works(self): need = self.opt.associated(["EX_CORDA_0"]) - self.assertEqual(list(need["EX_CORDA_0"]), ["EX_A", "r1"]) + solutions = (["EX_A", "r1"], ["EX_B", "r2"]) + self.assertTrue(list(need["EX_CORDA_0"]) in solutions) def test_redundancy_works(self): conf = self.opt.conf.copy() @@ -121,13 +126,16 @@ def test_redundancy_works(self): need = self.opt.associated(["EX_CORDA_0"], conf) self.assertEqual(len(need["EX_CORDA_0"]), 2) + class TestCORDAlarge(unittest.TestCase): def setUp(self): model = test_model() conf = {} for i, r in enumerate(model.reactions): - if i % 2 == 0: conf[r.id] = -1 - else: conf[r.id] = 2 + if i % 2 == 0: + conf[r.id] = -1 + else: + conf[r.id] = 2 conf["r60"] = 3 self.model = model self.conf = conf @@ -143,7 +151,8 @@ def test_association_work(self): def test_conf_vals(self): conf = self.conf.copy() - for r in self.model.reactions: conf[r.id] = 1 + for r in self.model.reactions: + conf[r.id] = 1 conf["r60"] = 3 conf["r42"] = 0 conf["r12"] = 1