diff --git a/README.md b/README.md
index 1dbe7733..7d138bf3 100644
--- a/README.md
+++ b/README.md
@@ -218,6 +218,27 @@ max_biomass_objective = fva_output[3]
The output of FVA method is tuple that contains `numpy` arrays. The vectors `min_fluxes` and `max_fluxes` contains the minimum and the maximum values of each flux. The vector `max_biomass_flux_vector` is the optimal flux vector according to the biomass objective function and `max_biomass_objective` is the value of that optimal solution.
+```python
+fva_output_df = model.fva_to_df()
+```
+
+returns the solution is an easier to query dataframe with the model's reactions as indices:
+```python
+>>> fva_output_df
+ minimum maximum
+PFK 7.477306 7.477485
+PFL 0.000000 0.000066
+PGI 4.860632 4.861110
+PGK -16.023610 -16.023451
+PGL 4.959736 4.960214
+>>>
+>>> df.loc["NH4t"]
+minimum 4.765316
+maximum 4.765324
+Name: NH4t, dtype: float64
+```
+
+
To apply FBA method,
```python
@@ -229,6 +250,16 @@ max_biomass_objective = fba_output[1]
while the output vectors are the same with the previous example.
+Again, one may use the `fba_to_df()` to get the solution as a pandas dataframe with model's reactions as indices.
+```python
+>>> model.fba_to_df()
+ fluxes
+PFK 7.477382
+PFL 0.000000
+PGI 4.860861
+PGK -16.023526
+PGL 4.959985
+```
### Set the restriction in the flux space
@@ -266,6 +297,46 @@ sampler = polytope_sampler(model)
steady_states = sampler.generate_steady_states()
```
+### Change the medium on your model
+
+To do that, you need to first describe the new medium and then assign it on your `model`.
+For example:
+
+```python
+initial_medium = model.medium
+model.medium
+# {'EX_co2_e': 1000.0, 'EX_glc__D_e': 10.0, 'EX_h_e': 1000.0, 'EX_h2o_e': 1000.0, 'EX_nh4_e': 1000.0, 'EX_o2_e': 1000.0, 'EX_pi_e': 1000.0}
+model.fba()[-1]
+# 0.8739215069684305
+new_medium = initial_medium.copy()
+
+# Set anoxygenic conditions
+new_medium['EX_o2_e'] = 0
+model.medium = new_medium
+model.fba()[-1]
+
+# Check the difference in the optimal value
+# 0.21166294973531055
+```
+
+### Who-is-who
+
+Models may use ids for metabolites and reactions hard to interpret.
+You may use the `reactions_map` and the `metabolites_map` that return the reactions/metabolites ids along with their corresponding names.
+For example:
+
+```python
+>>> model.reactions_map
+ reaction_name
+PFK Phosphofructokinase
+PFL Pyruvate formate lyase
+PGI Glucose-6-phosphate isomerase
+PGK Phosphoglycerate kinase
+PGL 6-phosphogluconolactonase
+```
+
+
+
### Plot flux marginals
@@ -282,8 +353,8 @@ steady_states = sampler.generate_steady_states(ess = 3000)
# plot the histogram for the 14th reaction in e-coli (ACONTa)
reactions = model.reactions
plot_histogram(
- steady_states[13],
- reactions[13],
+ steady_states.loc["ACONTa"],
+ "ACONTa",
n_bins = 60,
)
```
@@ -306,8 +377,8 @@ steady_states = sampler.generate_steady_states(ess = 3000)
# plot the copula between the 13th (PPC) and the 14th (ACONTa) reaction in e-coli
reactions = model.reactions
-data_flux2=[steady_states[12],reactions[12]]
-data_flux1=[steady_states[13],reactions[13]]
+data_flux2=[steady_states.loc["ACONTa"], "ACONTa"]
+data_flux1=[steady_states.loc["PPC"], "PPC"]
plot_copula(data_flux1, data_flux2, n=10)
```
diff --git a/dingo/MetabolicNetwork.py b/dingo/MetabolicNetwork.py
index eec1a8b5..37d2ef6a 100644
--- a/dingo/MetabolicNetwork.py
+++ b/dingo/MetabolicNetwork.py
@@ -7,12 +7,15 @@
# Licensed under GNU LGPL.3, see LICENCE file
import numpy as np
+import pandas as pd
import sys
from typing import Dict
import cobra
from dingo.loading_models import read_json_file, read_mat_file, read_sbml_file, parse_cobra_model
from dingo.fva import slow_fva
from dingo.fba import slow_fba
+import logging
+logger = logging.getLogger(__name__)
try:
import gurobipy
@@ -33,10 +36,12 @@ def __init__(self, tuple_args):
import gurobipy
self._parameters["fast_computations"] = True
+ self._parameters["tol"] = 1e-06
except ImportError as e:
self._parameters["fast_computations"] = False
+ self._parameters["tol"] = 1e-03
- if len(tuple_args) != 10:
+ if len(tuple_args) != 12:
raise Exception(
"An unknown input format given to initialize a metabolic network object."
)
@@ -51,6 +56,8 @@ def __init__(self, tuple_args):
self._medium = tuple_args[7]
self._medium_indices = tuple_args[8]
self._exchanges = tuple_args[9]
+ self._reactions_map = tuple_args[10]
+ self._metabolites_map = tuple_args[11]
try:
if self._biomass_index is not None and (
@@ -104,11 +111,11 @@ def from_cobra_model(cls, arg):
return cls(parse_cobra_model(arg))
- def fva(self):
+ def _fva(self):
"""A member function to apply the FVA method on the metabolic network."""
if self._parameters["fast_computations"]:
- return fast_fva(
+ min_fluxes, max_fluxes, max_biomass_flux_vector, max_biomass_objective = fast_fva(
self._lb,
self._ub,
self._S,
@@ -116,22 +123,43 @@ def fva(self):
self._parameters["opt_percentage"],
)
else:
- return slow_fva(
+ min_fluxes, max_fluxes, max_biomass_flux_vector, max_biomass_objective = slow_fva(
self._lb,
self._ub,
self._S,
self._biomass_function,
self._parameters["opt_percentage"],
)
+ self._min_fluxes = min_fluxes
+ self._max_fluxes = max_fluxes
+ self._opt_vector = max_biomass_flux_vector
+ self._opt_value = max_biomass_objective
+ return min_fluxes, max_fluxes, max_biomass_flux_vector, max_biomass_objective
- def fba(self):
+ def _fba(self):
"""A member function to apply the FBA method on the metabolic network."""
if self._parameters["fast_computations"]:
- return fast_fba(self._lb, self._ub, self._S, self._biomass_function)
+ opt_vector, opt_value = fast_fba(self._lb, self._ub, self._S, self._biomass_function)
else:
- return slow_fba(self._lb, self._ub, self._S, self._biomass_function)
+ opt_vector, opt_value = slow_fba(self._lb, self._ub, self._S, self._biomass_function)
+ self._opt_vector = opt_vector
+ self._opt_value = opt_value
+ return opt_vector, opt_value
+
+ def fba(self):
+ if not hasattr(self, '_opt_vector'):
+ self._fba()
+ fba_df = pd.DataFrame({'fluxes': self._opt_vector}, index=self._reactions)
+ return fba_df
+
+ def fva(self):
+ if not hasattr(self, '_min_fluxes'):
+ self._fva()
+ fva_df = pd.DataFrame({'minimum': self._min_fluxes, 'maximum': self._max_fluxes}, index=self._reactions)
+ return fva_df
+ # Descriptors
@property
def lb(self):
return self._lb
@@ -172,6 +200,30 @@ def exchanges(self):
def parameters(self):
return self._parameters
+ @property
+ def reactions_map(self):
+ return self._reactions_map
+
+ @property
+ def metabolites_map(self):
+ return self._metabolites_map
+
+ @property
+ def opt_value(self):
+ return self._opt_value
+
+ @property
+ def opt_vector(self):
+ return self._opt_vector
+
+ @property
+ def min_fluxes(self):
+ return self._min_fluxes
+
+ @property
+ def max_fluxes(self):
+ return self._max_fluxes
+
@property
def get_as_tuple(self):
return (
@@ -183,16 +235,11 @@ def get_as_tuple(self):
self._biomass_index,
self._biomass_function,
self._medium,
- self._inter_medium,
- self._exchanges
+ self._exchanges,
+ self._reactions_map,
+ self._metabolites_map
)
- def num_of_reactions(self):
- return len(self._reactions)
-
- def num_of_metabolites(self):
- return len(self._metabolites)
-
@lb.setter
def lb(self, value):
self._lb = value
@@ -221,7 +268,6 @@ def biomass_index(self, value):
def biomass_function(self, value):
self._biomass_function = value
-
@medium.setter
def medium(self, medium: Dict[str, float]) -> None:
"""Set the constraints on the model exchanges.
@@ -275,17 +321,51 @@ def set_active_bound(reaction: str, reac_index: int, bound: float) -> None:
# Turn off reactions not present in media
for rxn_id in exchange_rxns - frozen_media_rxns:
"""
- is_export for us, needs to check on the S
- order reactions to their lb and ub
+ is_export for us, needs to check on the S
+ order reactions to their lb and ub
"""
# is_export = rxn.reactants and not rxn.products
reac_index = self._reactions.index(rxn_id)
- products = np.any(self._S[:,reac_index] > 0)
+ products = np.any(self._S[:,reac_index] > 0)
reactants_exist = np.any(self._S[:,reac_index] < 0)
is_export = True if not products and reactants_exist else False
set_active_bound(
rxn_id, reac_index, min(0.0, -self._lb[reac_index] if is_export else self._ub[reac_index])
)
+ self._medium = medium
+ self._opt_vector = None
+
+ @reactions_map.setter
+ def reactions_map(self, value):
+ self._reactions_map = value
+
+ @metabolites_map.setter
+ def metabolites_map(self, value):
+ self._metabolites_map = value
+
+ @min_fluxes.setter
+ def min_fluxes(self, value):
+ self._min_fluxes = value
+
+ @max_fluxes.setter
+ def max_fluxes(self, value):
+ self._max_fluxes = value
+
+ @opt_value.setter
+ def opt_value(self, value):
+ self._opt_value = value
+
+ @opt_vector.setter
+ def opt_vector(self, value):
+ self._opt_vector = value
+
+
+ # Attributes
+ def num_of_reactions(self):
+ return len(self._reactions)
+
+ def num_of_metabolites(self):
+ return len(self._metabolites)
def set_fast_mode(self):
diff --git a/dingo/PolytopeSampler.py b/dingo/PolytopeSampler.py
index 4fe9dca7..e3ba537a 100644
--- a/dingo/PolytopeSampler.py
+++ b/dingo/PolytopeSampler.py
@@ -6,6 +6,7 @@
# Licensed under GNU LGPL.3, see LICENCE file
+import copy
import numpy as np
import warnings
import math
@@ -16,6 +17,7 @@
get_matrices_of_low_dim_polytope,
get_matrices_of_full_dim_polytope,
)
+import pandas as pd
try:
import gurobipy
@@ -37,13 +39,13 @@ def __init__(self, metabol_net):
if not isinstance(metabol_net, MetabolicNetwork):
raise Exception("An unknown input object given for initialization.")
- self._metabolic_network = metabol_net
- self._A = []
- self._b = []
- self._N = []
- self._N_shift = []
- self._T = []
- self._T_shift = []
+ self._metabolic_network = copy.deepcopy(metabol_net)
+ self._A = np.empty( shape=(0, 0) )
+ self._b = np.empty( shape=(0, 0) )
+ self._N = np.empty( shape=(0, 0) )
+ self._N_shift = np.empty( shape=(0, 0) )
+ self._T = np.empty( shape=(0, 0) )
+ self._T_shift = np.empty( shape=(0, 0) )
self._parameters = {}
self._parameters["nullspace_method"] = "sparseQR"
self._parameters["opt_percentage"] = self.metabolic_network.parameters[
@@ -68,18 +70,18 @@ def get_polytope(self):
"""
if (
- self._A == []
- or self._b == []
- or self._N == []
- or self._N_shift == []
- or self._T == []
- or self._T_shift == []
+ self._A.size == 0
+ or self._b.size == 0
+ or self._N.size == 0
+ or self._N_shift.size == 0
+ or self._T.size == 0
+ or self._T_shift.size == 0
):
(
max_biomass_flux_vector,
max_biomass_objective,
- ) = self._metabolic_network.fba()
+ ) = self._metabolic_network._fba()
if (
self._parameters["fast_computations"]
@@ -106,7 +108,7 @@ def get_polytope(self):
max_fluxes,
max_biomass_flux_vector,
max_biomass_objective,
- ) = self._metabolic_network.fva()
+ ) = self._metabolic_network._fva()
A, b, Aeq, beq = get_matrices_of_low_dim_polytope(
self._metabolic_network.S,
@@ -166,6 +168,7 @@ def generate_steady_states(
self._A, self._b, Tr, Tr_shift, samples = P.fast_mmcs(
ess, psrf, parallel_mmcs, num_threads
)
+
else:
self._A, self._b, Tr, Tr_shift, samples = P.slow_mmcs(
ess, psrf, parallel_mmcs, num_threads
@@ -184,8 +187,10 @@ def generate_steady_states(
self._T = np.dot(self._T, Tr)
self._T_shift = np.add(self._T_shift, Tr_shift)
- return steady_states
-
+ steady_states_df = pd.DataFrame(steady_states, index = self._metabolic_network.reactions)
+
+ return steady_states_df
+
def generate_steady_states_no_multiphase(
self, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None
):
@@ -197,11 +202,11 @@ def generate_steady_states_no_multiphase(
burn_in -- the number of points to burn before sampling
thinning -- the walk length of the chain
"""
-
+
self.get_polytope()
P = HPolytope(self._A, self._b)
-
+
if bias_vector is None:
bias_vector = np.ones(self._A.shape[1], dtype=np.float64)
else:
@@ -213,8 +218,9 @@ def generate_steady_states_no_multiphase(
steady_states = map_samples_to_steady_states(
samples_T, self._N, self._N_shift
)
+ steady_states_df = pd.DataFrame(steady_states, index = self._metabolic_network.reactions)
- return steady_states
+ return steady_states_df
@staticmethod
def sample_from_polytope(
@@ -245,7 +251,7 @@ def sample_from_polytope(
)
return samples
-
+
@staticmethod
def sample_from_polytope_no_multiphase(
A, b, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None
@@ -264,7 +270,7 @@ def sample_from_polytope_no_multiphase(
bias_vector = np.ones(A.shape[1], dtype=np.float64)
else:
bias_vector = bias_vector.astype('float64')
-
+
P = HPolytope(A, b)
try:
@@ -286,16 +292,13 @@ def round_polytope(
A, b, Tr, Tr_shift, round_value = P.rounding(method, True)
except ImportError as e:
A, b, Tr, Tr_shift, round_value = P.rounding(method, False)
-
+
return A, b, Tr, Tr_shift
+
@staticmethod
def sample_from_fva_output(
- min_fluxes,
- max_fluxes,
- biomass_function,
- max_biomass_objective,
- S,
+ model,
opt_percentage=100,
ess=1000,
psrf=False,
@@ -305,11 +308,7 @@ def sample_from_fva_output(
"""A static function to sample steady states when the output of FVA is given.
Keyword arguments:
- min_fluxes -- minimum values of the fluxes, i.e., a n-dimensional vector
- max_fluxes -- maximum values for the fluxes, i.e., a n-dimensional vector
- biomass_function -- the biomass objective function
- max_biomass_objective -- the maximum value of the biomass objective function
- S -- stoichiometric matrix
+ model -- a dingo.MetabolicNetwork() object
opt_percentage -- consider solutions that give you at least a certain
percentage of the optimal solution (default is to consider
optimal solutions only)
@@ -319,16 +318,18 @@ def sample_from_fva_output(
num_threads -- the number of threads to use for parallel mmcs
"""
+ min_fluxes, max_fluxes, opt_vector, opt_value = model._fva()
+
A, b, Aeq, beq = get_matrices_of_low_dim_polytope(
- S, min_fluxes, max_fluxes, opt_percentage, tol
+ model.S, min_fluxes, max_fluxes, opt_percentage, model._parameters["tol"]
)
- A = np.vstack((A, -biomass_function))
+ A = np.vstack((A, -model.biomass_function))
b = np.append(
b,
-(opt_percentage / 100)
- * self._parameters["tol"]
- * math.floor(max_biomass_objective / self._parameters["tol"]),
+ * model._parameters["tol"]
+ * math.floor(opt_value / model._parameters["tol"]),
)
A, b, N, N_shift = get_matrices_of_full_dim_polytope(A, b, Aeq, beq)
@@ -350,6 +351,17 @@ def sample_from_fva_output(
return steady_states
+ @staticmethod
+ def samples_as_df(model, samples):
+ """A static function to convert the samples numpy ndarray to a pandas DataFrame with model's reactions as indices
+
+ Keyword arguments:
+ model --
+ samples --
+ """
+ samples_df = pd.DataFrame(samples, index = model.reactions)
+ return samples_df
+
@property
def A(self):
return self._A
@@ -411,3 +423,4 @@ def set_tol(self, value):
def set_opt_percentage(self, value):
self._parameters["opt_percentage"] = value
+
diff --git a/dingo/__init__.py b/dingo/__init__.py
index 49f7d043..b2986eee 100644
--- a/dingo/__init__.py
+++ b/dingo/__init__.py
@@ -169,7 +169,7 @@ def dingo_main():
else:
raise Exception("An unknown format file given.")
- result_obj = model.fba()
+ result_obj = model._fba()
with open("dingo_fba_" + name + ".pckl", "wb") as dingo_fba_file:
pickle.dump(result_obj, dingo_fba_file)
diff --git a/dingo/loading_models.py b/dingo/loading_models.py
index 84efceb5..d1d714d9 100644
--- a/dingo/loading_models.py
+++ b/dingo/loading_models.py
@@ -9,6 +9,7 @@
import json
import numpy as np
import cobra
+import pandas as pd
def read_json_file(input_file):
"""A Python function to Read a Bigg json file and returns,
@@ -23,7 +24,7 @@ def read_json_file(input_file):
input_file -- a json file that contains the information about a mettabolic network, for example see http://bigg.ucsd.edu/models
"""
- try:
+ try:
cobra.io.load_matlab_model( input_file )
except:
cobra_config = cobra.Configuration()
@@ -45,7 +46,7 @@ def read_mat_file(input_file):
Keyword arguments:
input_file -- a mat file that contains a MATLAB structure with the information about a mettabolic network, for example see http://bigg.ucsd.edu/models
"""
- try:
+ try:
cobra.io.load_matlab_model( input_file )
except:
cobra_config = cobra.Configuration()
@@ -56,8 +57,8 @@ def read_mat_file(input_file):
return (parse_cobra_model( model ))
def read_sbml_file(input_file):
- """A Python function, based on the cobra.io.read_sbml_model() function of cabrapy
- and the extract_polytope() function of PolyRound
+ """A Python function, based on the cobra.io.read_sbml_model() function of cobrapy
+ and the extract_polytope() function of PolyRound
(https://gitlab.com/csb.ethz/PolyRound/-/blob/master/PolyRound/static_classes/parse_sbml_stoichiometry.py)
to read an SBML file (.xml) and return:
(a) lower/upper flux bounds
@@ -68,10 +69,10 @@ def read_sbml_file(input_file):
(f) the objective function to maximize the biomass pseudoreaction
Keyword arguments:
- input_file -- a xml file that contains an SBML model with the information about a mettabolic network, for example see:
+ input_file -- a xml file that contains an SBML model with the information about a mettabolic network, for example see:
https://github.com/VirtualMetabolicHuman/AGORA/blob/master/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/Abiotrophia_defectiva_ATCC_49176.xml
"""
- try:
+ try:
cobra.io.read_sbml_model( input_file )
except:
cobra_config = cobra.Configuration()
@@ -137,9 +138,11 @@ def parse_cobra_model(cobra_model):
exchanges.append(reac.id)
- return lb, ub, S, metabolites, reactions, biomass_index, biomass_function, medium, inter_medium, exchanges
-
-
-
-
+ # Map ids to complete names for convenience
+ reactions_map = pd.DataFrame( [x.name for x in cobra_model.reactions], [x.id for x in cobra_model.reactions])
+ reactions_map.columns = ["reaction_name"]
+ metabolites_map = pd.DataFrame( [x.name for x in cobra_model.metabolites], [x.id for x in cobra_model.metabolites])
+ metabolites_map.columns = ["metabolite_name"]
+ return lb, ub, S, metabolites, reactions, biomass_index, biomass_function, \
+ medium, inter_medium, exchanges, reactions_map, metabolites_map
diff --git a/dingo/utils.py b/dingo/utils.py
index 8074b155..546774e9 100644
--- a/dingo/utils.py
+++ b/dingo/utils.py
@@ -34,14 +34,14 @@ def compute_copula(flux1, flux2, n):
rng = range((j*math.floor(N/n)),((j+1)*math.floor(N/n)))
grouped_flux1[I1[rng]] = j
grouped_flux2[I2[rng]] = j
-
+
for i in range(n):
for j in range(n):
copula[i,j] = sum((grouped_flux1==i) *( grouped_flux2==j))
-
+
copula = copula / N
return copula
-
+
def apply_scaling(A, b, cs, rs):
"""A Python function to apply the scaling computed by the function `gmscale` to a convex polytope
diff --git a/tests/fast_implementation_test.py b/tests/fast_implementation_test.py
index a4b1995f..ec2992a8 100644
--- a/tests/fast_implementation_test.py
+++ b/tests/fast_implementation_test.py
@@ -12,6 +12,7 @@
from dingo.gurobi_based_implementations import fast_inner_ball
class TestFastMethods(unittest.TestCase):
+
def test_fast_max_bal_computation(self):
m = 2
@@ -26,6 +27,7 @@ def test_fast_max_bal_computation(self):
self.assertTrue(abs(max_ball[1] - 1) < 1e-08)
+
def test_fast_fva(self):
current_directory = os.getcwd()
@@ -34,12 +36,18 @@ def test_fast_fva(self):
model = MetabolicNetwork.from_json(input_file_json)
model.set_fast_mode()
- res = model.fva()
+ res = model._fva()
self.assertTrue(abs(res[3] - 0.8739215069684305) < 1e-08)
self.assertEqual(res[0].size, 95)
self.assertEqual(res[1].size, 95)
+ fva_df = model.fva()
+ biomass_function = model.reactions[model.biomass_index]
+ self.assertTrue(fva_df.loc[biomass_function]["maximum"] - fva_df.loc[biomass_function]["minimum"] <= 1e-03)
+ self.assertEqual(fva_df.shape, (95,2))
+
+
def test_ecoli_to_full_dimensional_polytope(self):
current_directory = os.getcwd()
@@ -55,9 +63,9 @@ def test_ecoli_to_full_dimensional_polytope(self):
self.assertEqual(sampler.A.shape[0], 26)
self.assertEqual(sampler.A.shape[1], 24)
-
self.assertEqual(steady_states.shape[0], 95)
+
def test_fast_fba(self):
current_directory = os.getcwd()
@@ -65,11 +73,13 @@ def test_fast_fba(self):
model = MetabolicNetwork.from_json(input_file_json)
model.set_fast_mode()
-
- res = model.fba()
-
+ res = model._fba()
self.assertTrue(abs(res[1] - 0.8739215069684305) < 1e-08)
+ fba_df = model.fba()
+ biomass_function = model.reactions[model.biomass_index]
+ self.assertTrue(abs(fba_df.loc[biomass_function]["fluxes"] - 0.8739215069684305) < 1e-08)
+
if __name__ == "__main__":
unittest.main()
diff --git a/tests/fba.py b/tests/fba.py
index 94dc5b6f..1f517f66 100644
--- a/tests/fba.py
+++ b/tests/fba.py
@@ -16,39 +16,51 @@ def test_fba_json(self):
input_file_json = os.getcwd() + "/ext_data/e_coli_core.json"
model = MetabolicNetwork.from_json(input_file_json)
model.set_slow_mode()
- res = model.fba()
-
+ res = model._fba()
self.assertTrue(abs(res[1] - 0.8739215067486387) < 1e-03)
+ biomass_function = model.reactions[model.biomass_index]
+ fba_df = model.fba()
+ self.assertTrue(abs(fba_df.loc[biomass_function]["fluxes"] - 0.8739215067486387) < 1e-03)
+
def test_fba_mat(self):
input_file_mat = os.getcwd() + "/ext_data/e_coli_core.mat"
model = MetabolicNetwork.from_mat(input_file_mat)
model.set_slow_mode()
-
- res = model.fba()
-
+ res = model._fba()
self.assertTrue(abs(res[1] - 0.8739215067486387) < 1e-03)
+ biomass_function = model.reactions[model.biomass_index]
+ fba_df = model.fba()
+ self.assertTrue(abs(fba_df.loc[biomass_function]["fluxes"] - 0.8739215067486387) < 1e-03)
+
def test_fba_sbml(self):
input_file_sbml = os.getcwd() + "/ext_data/e_coli_core.xml"
model = MetabolicNetwork.from_sbml(input_file_sbml)
model.set_slow_mode()
-
- res = model.fba()
-
+ res = model._fba()
self.assertTrue(abs(res[1] - 0.8739215067486387) < 1e-03)
+ biomass_function = model.reactions[model.biomass_index]
+ fba_df = model.fba()
+ self.assertTrue(abs(fba_df.loc[biomass_function]["fluxes"] - 0.8739215067486387) < 1e-03)
+
def test_modify_medium(self):
input_file_sbml = os.getcwd() + "/ext_data/e_coli_core.xml"
model = MetabolicNetwork.from_sbml(input_file_sbml)
- model.set_slow_mode()
+
+ # Check if script is running in GitHub action
+ if os.getenv('CI', 'false').lower() == 'true':
+ model.set_slow_mode()
initial_medium = model.medium
- initial_fba = model.fba()[-1]
+ initial_fba = model._fba()[-1]
+ initial_fba_df = model.fba()
+ # Original indices of the exchange reactions
e_coli_core_medium_compound_indices = {
"EX_co2_e" : 46,
"EX_glc__D_e" : 51,
@@ -59,25 +71,35 @@ def test_modify_medium(self):
"EX_pi_e" : 60
}
- glc_index = model.reactions.index("EX_glc__D_e")
- o2_index = model.reactions.index("EX_o2_e")
-
+ # Edit and assign new medium to model
new_media = initial_medium.copy()
- new_media["EX_glc__D_e"] = 1.5
- new_media["EX_o2_e"] = -0.5
-
+ new_media["EX_glc__D_e"] = 35
+ new_media["EX_o2_e"] = 0.5
model.medium = new_media
+ # Check if indices are affected
updated_media = model.medium
updated_medium_indices = {}
for reac in updated_media:
updated_medium_indices[reac] = model.reactions.index(reac)
+ glc_index = model.reactions.index("EX_glc__D_e")
+ o2_index = model.reactions.index("EX_o2_e")
+
self.assertTrue(updated_medium_indices == e_coli_core_medium_compound_indices)
+ self.assertTrue(model.lb[glc_index] == -35 and model.lb[o2_index] == -0.5)
+
+ # Check if optimal value is affected
+ self.assertTrue(
+ abs((model._fba()[-1] - initial_fba) - 0.1172) <= 1e-03
+ )
- self.assertTrue(model.lb[glc_index] == -1.5 and model.lb[o2_index] == 0.5)
- self.assertTrue(initial_fba - model.fba()[-1] > 0)
+ biomass_function = model.reactions[model.biomass_index]
+ fba_df = model.fba()
+ self.assertTrue(
+ abs((fba_df.loc[biomass_function]["fluxes"] - initial_fba_df.loc[biomass_function]["fluxes"]) - 0.1172) <= 1e-03
+ )
if __name__ == "__main__":
diff --git a/tests/sampling.py b/tests/sampling.py
index bf731715..39302a36 100644
--- a/tests/sampling.py
+++ b/tests/sampling.py
@@ -23,7 +23,7 @@ def test_sample_json(self):
steady_states = sampler.generate_steady_states(ess = 20000, psrf = True)
self.assertTrue( steady_states.shape[0] == 95 )
- self.assertTrue( abs( steady_states[12].mean() - 2.504 ) < 1e-02 )
+ self.assertTrue( abs( steady_states.loc["PPC"].mean() - 2.504 ) < 1e-02 )
def test_sample_mat(self):
@@ -32,10 +32,10 @@ def test_sample_mat(self):
model = MetabolicNetwork.from_mat(input_file_mat)
sampler = PolytopeSampler(model)
- steady_states = sampler.generate_steady_states(ess = 20000, psrf = True)
+ steady_states = sampler.generate_steady_states(ess = 20000, psrf = True)
self.assertTrue( steady_states.shape[0] == 95 )
- self.assertTrue( abs( steady_states[12].mean() - 2.504 ) < 1e-02 )
+ self.assertTrue( abs( steady_states.loc["PPC"].mean() - 2.504 ) < 1e-02 )
def test_sample_sbml(self):
@@ -47,7 +47,7 @@ def test_sample_sbml(self):
steady_states = sampler.generate_steady_states(ess = 20000, psrf = True)
self.assertTrue( steady_states.shape[0] == 95 )
- self.assertTrue( abs( steady_states[12].mean() - 2.504 ) < 1e-02 )
+ self.assertTrue( abs( steady_states.loc["PPC"].mean() - 2.504 ) < 1e-02 )
diff --git a/tutorials/dingo_tutorial.ipynb b/tutorials/dingo_tutorial.ipynb
index 1b1b7e5a..5c7799ee 100644
--- a/tutorials/dingo_tutorial.ipynb
+++ b/tutorials/dingo_tutorial.ipynb
@@ -689,7 +689,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 3,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
@@ -714,7 +714,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 2,
"metadata": {
"id": "FahzUpQFommK"
},
@@ -735,13 +735,13 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 4,
"metadata": {
"id": "NNSJM_Z3McSs"
},
"outputs": [],
"source": [
- "model = MetabolicNetwork.from_json('ext_data/e_coli_core.json')"
+ "model = MetabolicNetwork.from_json('../ext_data/e_coli_core.json')"
]
},
{
@@ -1043,7 +1043,7 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 5,
"metadata": {
"id": "kWFbcJumShgh"
},
@@ -1053,6 +1053,120 @@
"fba_output = model.fba()"
]
},
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fba_output_df = model.fba_to_df()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " fluxes | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | PFK | \n",
+ " 7.477382 | \n",
+ "
\n",
+ " \n",
+ " | PFL | \n",
+ " 0.000000 | \n",
+ "
\n",
+ " \n",
+ " | PGI | \n",
+ " 4.860861 | \n",
+ "
\n",
+ " \n",
+ " | PGK | \n",
+ " -16.023526 | \n",
+ "
\n",
+ " \n",
+ " | PGL | \n",
+ " 4.959985 | \n",
+ "
\n",
+ " \n",
+ " | ... | \n",
+ " ... | \n",
+ "
\n",
+ " \n",
+ " | NADH16 | \n",
+ " 38.534610 | \n",
+ "
\n",
+ " \n",
+ " | NADTRHD | \n",
+ " 0.000000 | \n",
+ "
\n",
+ " \n",
+ " | NH4t | \n",
+ " 4.765319 | \n",
+ "
\n",
+ " \n",
+ " | O2t | \n",
+ " 21.799493 | \n",
+ "
\n",
+ " \n",
+ " | PDH | \n",
+ " 9.282533 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
95 rows × 1 columns
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " minimum | \n",
+ " maximum | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | PFK | \n",
+ " 176.609999 | \n",
+ " 1.766100e+02 | \n",
+ "
\n",
+ " \n",
+ " | PFL | \n",
+ " 0.000000 | \n",
+ " 6.666667e-07 | \n",
+ "
\n",
+ " \n",
+ " | PGI | \n",
+ " 9.999999 | \n",
+ " 1.000000e+01 | \n",
+ "
\n",
+ " \n",
+ " | PGK | \n",
+ " -20.000000 | \n",
+ " -2.000000e+01 | \n",
+ "
\n",
+ " \n",
+ " | PGL | \n",
+ " 0.000000 | \n",
+ " 1.333333e-06 | \n",
+ "
\n",
+ " \n",
+ " | ... | \n",
+ " ... | \n",
+ " ... | \n",
+ "
\n",
+ " \n",
+ " | NADH16 | \n",
+ " 99.999999 | \n",
+ " 1.000000e+02 | \n",
+ "
\n",
+ " \n",
+ " | NADTRHD | \n",
+ " 19.999999 | \n",
+ " 2.000000e+01 | \n",
+ "
\n",
+ " \n",
+ " | NH4t | \n",
+ " 0.000000 | \n",
+ " 1.777778e-07 | \n",
+ "
\n",
+ " \n",
+ " | O2t | \n",
+ " 60.000000 | \n",
+ " 6.000000e+01 | \n",
+ "
\n",
+ " \n",
+ " | PDH | \n",
+ " 19.999999 | \n",
+ " 2.000000e+01 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
95 rows × 2 columns
\n",
+ "
\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " 0 | \n",
+ " 1 | \n",
+ " 2 | \n",
+ " 3 | \n",
+ " 4 | \n",
+ " 5 | \n",
+ " 6 | \n",
+ " 7 | \n",
+ " 8 | \n",
+ " 9 | \n",
+ " ... | \n",
+ " 4990 | \n",
+ " 4991 | \n",
+ " 4992 | \n",
+ " 4993 | \n",
+ " 4994 | \n",
+ " 4995 | \n",
+ " 4996 | \n",
+ " 4997 | \n",
+ " 4998 | \n",
+ " 4999 | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | PFK | \n",
+ " 7.477382 | \n",
+ " 7.477380 | \n",
+ " 7.477381 | \n",
+ " 7.477385 | \n",
+ " 7.477383 | \n",
+ " 7.477383 | \n",
+ " 7.477387 | \n",
+ " 7.477399e+00 | \n",
+ " 7.477389 | \n",
+ " 7.477392e+00 | \n",
+ " ... | \n",
+ " 7.477376 | \n",
+ " 7.477384 | \n",
+ " 7.477380e+00 | \n",
+ " 7.477389e+00 | \n",
+ " 7.477388 | \n",
+ " 7.477392 | \n",
+ " 7.477388 | \n",
+ " 7.477384 | \n",
+ " 7.477385e+00 | \n",
+ " 7.477379 | \n",
+ "
\n",
+ " \n",
+ " | PFL | \n",
+ " 0.000002 | \n",
+ " 0.000002 | \n",
+ " 0.000002 | \n",
+ " 0.000003 | \n",
+ " 0.000002 | \n",
+ " 0.000005 | \n",
+ " 0.000009 | \n",
+ " 9.758857e-07 | \n",
+ " 0.000001 | \n",
+ " 4.354554e-06 | \n",
+ " ... | \n",
+ " 0.000003 | \n",
+ " 0.000001 | \n",
+ " 1.003201e-07 | \n",
+ " 2.492911e-06 | \n",
+ " 0.000004 | \n",
+ " 0.000004 | \n",
+ " 0.000004 | \n",
+ " 0.000001 | \n",
+ " 4.600969e-07 | \n",
+ " 0.000002 | \n",
+ "
\n",
+ " \n",
+ " | PGI | \n",
+ " 4.860848 | \n",
+ " 4.860838 | \n",
+ " 4.860842 | \n",
+ " 4.860858 | \n",
+ " 4.860858 | \n",
+ " 4.860857 | \n",
+ " 4.860859 | \n",
+ " 4.860871e+00 | \n",
+ " 4.860868 | \n",
+ " 4.860880e+00 | \n",
+ " ... | \n",
+ " 4.860839 | \n",
+ " 4.860861 | \n",
+ " 4.860854e+00 | \n",
+ " 4.860868e+00 | \n",
+ " 4.860872 | \n",
+ " 4.860889 | \n",
+ " 4.860876 | \n",
+ " 4.860865 | \n",
+ " 4.860866e+00 | \n",
+ " 4.860845 | \n",
+ "
\n",
+ " \n",
+ " | PGK | \n",
+ " -16.023522 | \n",
+ " -16.023519 | \n",
+ " -16.023521 | \n",
+ " -16.023526 | \n",
+ " -16.023526 | \n",
+ " -16.023526 | \n",
+ " -16.023527 | \n",
+ " -1.602353e+01 | \n",
+ " -16.023529 | \n",
+ " -1.602353e+01 | \n",
+ " ... | \n",
+ " -16.023520 | \n",
+ " -16.023527 | \n",
+ " -1.602352e+01 | \n",
+ " -1.602353e+01 | \n",
+ " -16.023531 | \n",
+ " -16.023536 | \n",
+ " -16.023532 | \n",
+ " -16.023528 | \n",
+ " -1.602353e+01 | \n",
+ " -16.023521 | \n",
+ "
\n",
+ " \n",
+ " | PGL | \n",
+ " 4.959998 | \n",
+ " 4.960008 | \n",
+ " 4.960004 | \n",
+ " 4.959988 | \n",
+ " 4.959988 | \n",
+ " 4.959989 | \n",
+ " 4.959987 | \n",
+ " 4.959975e+00 | \n",
+ " 4.959978 | \n",
+ " 4.959966e+00 | \n",
+ " ... | \n",
+ " 4.960007 | \n",
+ " 4.959985 | \n",
+ " 4.959992e+00 | \n",
+ " 4.959978e+00 | \n",
+ " 4.959974 | \n",
+ " 4.959957 | \n",
+ " 4.959970 | \n",
+ " 4.959981 | \n",
+ " 4.959980e+00 | \n",
+ " 4.960000 | \n",
+ "
\n",
+ " \n",
+ " | ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ " ... | \n",
+ "
\n",
+ " \n",
+ " | NADH16 | \n",
+ " 38.534624 | \n",
+ " 38.534630 | \n",
+ " 38.534628 | \n",
+ " 38.534622 | \n",
+ " 38.534622 | \n",
+ " 38.534614 | \n",
+ " 38.534627 | \n",
+ " 3.853460e+01 | \n",
+ " 38.534611 | \n",
+ " 3.853461e+01 | \n",
+ " ... | \n",
+ " 38.534623 | \n",
+ " 38.534618 | \n",
+ " 3.853463e+01 | \n",
+ " 3.853462e+01 | \n",
+ " 38.534622 | \n",
+ " 38.534613 | \n",
+ " 38.534616 | \n",
+ " 38.534611 | \n",
+ " 3.853461e+01 | \n",
+ " 38.534620 | \n",
+ "
\n",
+ " \n",
+ " | NADTRHD | \n",
+ " 0.000024 | \n",
+ " 0.000031 | \n",
+ " 0.000021 | \n",
+ " 0.000007 | \n",
+ " 0.000039 | \n",
+ " 0.000013 | \n",
+ " 0.000025 | \n",
+ " 7.147399e-06 | \n",
+ " 0.000010 | \n",
+ " 6.149728e-07 | \n",
+ " ... | \n",
+ " 0.000023 | \n",
+ " 0.000008 | \n",
+ " 2.870932e-05 | \n",
+ " 9.178589e-07 | \n",
+ " 0.000015 | \n",
+ " 0.000022 | \n",
+ " 0.000028 | \n",
+ " 0.000011 | \n",
+ " 4.563802e-06 | \n",
+ " 0.000014 | \n",
+ "
\n",
+ " \n",
+ " | NH4t | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317e+00 | \n",
+ " 4.765318 | \n",
+ " 4.765317e+00 | \n",
+ " ... | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317e+00 | \n",
+ " 4.765317e+00 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317 | \n",
+ " 4.765317e+00 | \n",
+ " 4.765317 | \n",
+ "
\n",
+ " \n",
+ " | O2t | \n",
+ " 21.799498 | \n",
+ " 21.799500 | \n",
+ " 21.799500 | \n",
+ " 21.799499 | \n",
+ " 21.799499 | \n",
+ " 21.799494 | \n",
+ " 21.799504 | \n",
+ " 2.179949e+01 | \n",
+ " 21.799494 | \n",
+ " 2.179950e+01 | \n",
+ " ... | \n",
+ " 21.799496 | \n",
+ " 21.799497 | \n",
+ " 2.179950e+01 | \n",
+ " 2.179950e+01 | \n",
+ " 21.799502 | \n",
+ " 21.799500 | \n",
+ " 21.799500 | \n",
+ " 21.799494 | \n",
+ " 2.179950e+01 | \n",
+ " 21.799495 | \n",
+ "
\n",
+ " \n",
+ " | PDH | \n",
+ " 9.282543 | \n",
+ " 9.282553 | \n",
+ " 9.282558 | \n",
+ " 9.282555 | \n",
+ " 9.282532 | \n",
+ " 9.282546 | \n",
+ " 9.282532 | \n",
+ " 9.282541e+00 | \n",
+ " 9.282552 | \n",
+ " 9.282542e+00 | \n",
+ " ... | \n",
+ " 9.282555 | \n",
+ " 9.282543 | \n",
+ " 9.282538e+00 | \n",
+ " 9.282542e+00 | \n",
+ " 9.282544 | \n",
+ " 9.282552 | \n",
+ " 9.282563 | \n",
+ " 9.282535 | \n",
+ " 9.282538e+00 | \n",
+ " 9.282569 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
95 rows × 5000 columns
\n",
+ "