2. Quick start¶
The INCAWrapper python wrapper to the MATLAB application INCA. INCA can apply multiple great 13C-MFA algorithms. Unfortunately, it is meant to be used through a GUI and therefore it is laborious to handle many samples and measurements and is hard in incorporate into high throughput pipelines. The INCAWrapper operates on both sides of INCA; it helps input model and data into INCA, it can run INCA algorithms and finally provides a workflow to export the results from INCA after execution of algorithms. This is a quick walk through of how to use the INCAWrapper to give interested users an overview of the capabilities.
[13]:
import pandas as pd
import pathlib
import incawrapper
from incawrapper import run_inca
import ast
The input data to the INCAWrapper is given as pandas dataframes. These has to comply to specified dataschemas (See :ref:input_data.ipynb for a thorough walk through of the requirements). The minimum inputs to run MFA in INCA is
Reaction data
Tracer data
For stationary 13C-MFA least one of
Flux measurements
MS measurements
For Isotopically Non-stationary 13C-MFA additionally measurements of the metabolite pool are can be used.
Let’s load some data to see the format of the input. All input data is validated for column name and data types, thus it will raise an error if either of those are incorrect.
[14]:
data_folder = pathlib.Path("./examples/Literature data/simple model")
tracers_data = pd.read_csv(data_folder / "tracers.csv",
converters={'atom_mdv':ast.literal_eval, 'atom_ids':ast.literal_eval} # a trick to read lists from csv
)
reactions_data = pd.read_csv(data_folder / "reactions.csv")
flux_data = pd.read_csv(data_folder / "flux_measurements.csv")
ms_data = pd.read_csv(data_folder / "ms_measurements.csv",
converters={'labelled_atom_ids': ast.literal_eval} # a trick to read lists from csv
)
We use a simple toy model with 5 reactions [1,2]. In this data set we have two experiments exp1 and exp2. They were conducted with different tracers. For exp1 we have both flux and MS measurements, while for exp2 the data set only contains MS measurements.
[15]:
reactions_data.head()
[15]:
| model | rxn_id | rxn_eqn | |
|---|---|---|---|
| 0 | simple_model | R1 | A (abc) -> B (abc) |
| 1 | simple_model | R2 | B (abc) <-> D (abc) |
| 2 | simple_model | R3 | B (abc) -> C (bc) + E (a) |
| 3 | simple_model | R4 | B (abc) + C (de) -> D (bcd) + E (a) + E (e) |
| 4 | simple_model | R5 | D (abc) -> F (abc) |
[16]:
tracers_data.head()
[16]:
| experiment_id | met_id | tracer_id | atom_ids | ratio | atom_mdv | enrichment | |
|---|---|---|---|---|---|---|---|
| 0 | exp1 | A | [2-13C]A | [2] | 1.0 | [0, 1] | 1 |
| 1 | exp2 | A | [1,2-13C]A | [1, 2] | 0.5 | [0.05, 0.95] | 1 |
[17]:
flux_data.head()
[17]:
| experiment_id | rxn_id | flux | flux_std_error | |
|---|---|---|---|---|
| 0 | exp1 | R1 | 10.0 | 0.00001 |
[18]:
ms_data.head()
[18]:
| experiment_id | met_id | ms_id | measurement_replicate | labelled_atom_ids | unlabelled_atoms | mass_isotope | intensity | intensity_std_error | time | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | exp1 | F | F1 | 1 | [1, 2, 3] | NaN | 0 | 0.0001 | 0.000002 | 0 |
| 1 | exp1 | F | F1 | 1 | [1, 2, 3] | NaN | 1 | 0.8008 | 0.016016 | 0 |
| 2 | exp1 | F | F1 | 1 | [1, 2, 3] | NaN | 2 | 0.1983 | 0.003966 | 0 |
| 3 | exp1 | F | F1 | 1 | [1, 2, 3] | NaN | 3 | 0.0009 | 0.000018 | 0 |
| 4 | exp2 | F | F1 | 1 | [1, 2, 3] | NaN | 0 | 0.0002 | 0.000002 | 0 |
With the data loaded we are ready to create the INCA script that specifies the model, data and algorithms that we wish to apply. We will specify that the estimation algorithm should do 5 restarts and that INCA should not do natural abundance correction. Finally, specify that we wish to run three of INCA’s algorithms: estimate, simulate and continuate and that the model including the results should be saved in a file called simple_model_quickstart.mat in the current working directory.
[19]:
output_file = pathlib.Path("./examples/Literature data/simple model/simple_model_quikstart.mat")
script = incawrapper.create_inca_script_from_data(reactions_data, tracers_data, flux_data, ms_data, experiment_ids=["exp1"])
script.add_to_block("options", incawrapper.define_options(fit_starts=5,sim_na=False))
script.add_to_block("runner", incawrapper.define_runner(output_file, run_estimate=True, run_simulation=True, run_continuation=True))
Now the script is ready to run in matlab.
[20]:
import dotenv
inca_directory = pathlib.Path(dotenv.get_key(dotenv.find_dotenv(), "INCA_base_directory"))
incawrapper.run_inca(script, INCA_base_directory=inca_directory)
INCA script saved to /var/folders/z6/mxpxh4k56tv0h0ff41vmx7gdwtlpvp/T/tmpkga0p7fp/inca_script.m.
Starting MATLAB engine...
ms_exp1 = 1x1 msdata object
fields: atoms id [idvs] more on state
F1
m = 1x1 model object
fields: [expts] [mets] notes [options] [rates] [states]
5 reactions (6 fluxes)
6 states (3 balanced, 1 source, 2 sink and 0 unbalanced)
6 metabolites
1 experiments
Directional
Iteration Residual Step-size derivative Lambda
0 9.99596e+11
1 8.14509e+11 0.0973 -9.02e+11 1.51438
2 1471.59 1 -311 1.51438
3 1686.79 1 2.43e+03 0.504795
4 1647.49 1 2.35e+03 1.00959
5 1451.25 1 1.97e+03 4.03836
6 219.371 1 517 7.70091
7 136.921 1 -3.99 2.56697
8 136.494 1 -0.0262 0.855657
Optimization terminated: Unconstrained optimum found.
Parameters converged to within tolerance.
Directional
Iteration Residual Step-size derivative Lambda
0 9.99432e+11
1 9.87336e+11 0.00607 -9.93e+11 20.194
2 683304 1 -1.48e+06 20.194
3 38384 1 -8.04e+04 6.73132
4 2747.89 1 -4.53e+03 2.24377
5 1582.96 1 2.13e+03 0.747925
6 3023.3 1 2.74e+03 0.701643
7 3007.12 1 2.73e+03 1.40329
8 2912.05 1 2.66e+03 5.61314
9 2176.43 1 1.99e+03 44.9051
10 522.602 1 -45 359.241
11 159.881 1 -49.3 265.304
12 139.894 1 -3.49 88.4348
13 136.623 1 -0.485 29.4783
14 136.493 1 -0.0147 9.82609
Optimization terminated: Unconstrained optimum found.
Parameters converged to within tolerance.
Directional
Iteration Residual Step-size derivative Lambda
0 9.98988e+11
1 9.95521e+11 0.0017 -1.02e+12 893.792
2 4.36045e+07 1 -3.46e+08 893.792
3 1.184e+07 1 -5.95e+06 297.931
4 2.18909e+06 1 -2.04e+06 269.491
5 121069 1 -2.74e+05 201.915
6 7335.28 1 -1.46e+04 67.305
7 1853.89 1 -530 22.435
8 2196.53 1 2.68e+03 7.47833
9 1958.04 1 2.05e+03 14.9567
10 1908.05 1 1.96e+03 59.8266
11 1510.62 1 1.31e+03 478.613
12 703.253 1 291 517.658
13 651.681 1 -1.66 209.579
14 638.339 1 -8.07 69.8595
15 640.145 1 344 23.2865
16 603.87 1 -22 46.573
17 5106.12 0.704 -2.4e+03 15.5243
18 894.548 1 1.78e+03 31.0487
19 575.861 1 -16 124.195
20 721.113 1 986 41.3982
21 517.689 1 -28.1 82.7965
22 5029.28 0.824 -688 27.5988
23 905.256 1 1.79e+03 55.1976
24 465.752 1 -26.7 220.791
25 605.331 1 821 73.5968
26 389.228 1 -23.5 147.194
27 2317.9 1 5.68e+03 49.0646
28 421.563 1 413 98.1291
29 324.497 1 -28.2 392.517
30 279.9 1 123 130.839
31 158.307 1 -19.3 130.84
32 137.439 1 -2.4 43.6134
33 136.491 1 -0.0161 14.5378
Optimization terminated: Unconstrained optimum found.
Parameters converged to within tolerance.
Directional
Iteration Residual Step-size derivative Lambda
0 9.99884e+11
1 9.99022e+11 0.000451 -9.78e+11 3225.6
2 5.96824e+08 1 -2.24e+10 3225.6
3 2.88806e+07 1 -6.28e+06 1075.2
4 1.18299e+07 1 -4.05e+06 358.4
5 762453 1 -1.69e+06 356.277
6 37981.6 1 -8.99e+04 118.759
7 3140.34 1 -4.27e+03 39.5863
8 1806.13 1 -89.6 13.1954
9 1946.85 1 1.85e+03 4.39848
10 1935.65 1 1.82e+03 8.79696
11 1904.62 1 1.78e+03 35.1879
12 1635.79 1 1.36e+03 281.503
13 758.985 1 364 378.714
14 695.907 1 -0.197 181.23
15 695.248 1 -0.338 60.4099
16 693.002 1 -1.32 20.1366
17 675.123 1 -16.8 6.71222
18 5193.72 0.277 -2.19e+04 2.23741
19 5191.48 0.561 -1.08e+04 4.47481
20 645.678 1 -22.3 17.8992
21 5155.13 0.433 -7.08e+03 5.96641
22 5151.7 0.8 -3.81e+03 11.9328
23 613.431 1 -20.9 47.7313
24 5116.55 0.772 -2.5e+03 15.9104
25 644.567 1 491 31.8209
26 592.081 1 -12.2 127.283
27 573.034 1 249 42.4278
28 3501.03 1 8.75e+03 45.2624
29 438.403 1 215 90.5249
30 272.772 1 134 90.3291
31 136.779 1 -2.94 89.3289
32 136.495 1 -0.0259 29.7763
Optimization terminated: Unconstrained optimum found.
Parameters converged to within tolerance.
Directional
Iteration Residual Step-size derivative Lambda
0 9.9856e+11
1 9.2428e+11 0.0379 -9.61e+11 31.4564
2 1.21036e+06 1 -4.17e+04 31.4564
3 83632.7 1 -1.51e+05 10.4855
4 4187.54 1 -9.23e+03 3.71603
5 1108.73 1 673 1.23868
6 2761.68 1 2.94e+03 0.412892
7 2750.64 1 2.93e+03 0.825784
8 2685.62 1 2.87e+03 3.30313
9 2167.77 1 2.39e+03 26.4251
10 711.742 1 536 211.401
11 147.484 1 73.4 211.682
12 136.832 1 0.35 70.5608
13 136.5 1 -0.0441 23.5203
Optimization terminated: Unconstrained optimum found.
Parameters converged to within tolerance.
Directional
Iteration Residual Step-size derivative Lambda
0 136.491
1 136.491 1 -2.81e-05 0.000358365
Optimization terminated: Unconstrained optimum found.
Parameters converged to within tolerance.
Estimation completed in 10.6800 seconds.
Preprocessing time: 7.1800 s
Computation time: 3.3500 s
Postprocessing time: 0.1500 s
========== Varying R1 upward from 10 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.768292 8.76523e-06 8.76523e-06 0.000000 0.000000 0 0 1 0.000358423
2 1.536584 1.23959e-05 3.63068e-06 -0.000000 0.000000 0 0 1 0.000119474
3 2.304875 1.51818e-05 2.78592e-06 0.000000 0.000000 0 0 1 3.98247e-05
4 3.073167 1.75305e-05 2.34864e-06 -0.000000 0.000000 0 0 1 1.32749e-05
5 3.841459 1.95996e-05 2.06919e-06 -0.000000 0.000000 0 0 1 4.42497e-06
6 4.609751 2.14703e-05 1.87069e-06 0.000000 0.000000 0 0 1 1.47499e-06
========== Varying R1 downward from 10 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.768292 8.76523e-06 8.76523e-06 0.000000 0.000000 0 0 1 0.000358423
2 1.536584 1.23959e-05 3.63068e-06 -0.000000 0.000000 0 0 1 0.000119474
3 2.304875 1.51818e-05 2.78592e-06 0.000000 0.000000 0 0 1 3.98247e-05
4 3.073167 1.75305e-05 2.34864e-06 0.000000 0.000000 0 0 1 1.32749e-05
5 3.841459 1.95996e-05 2.06919e-06 -0.000000 0.000000 0 0 1 4.42497e-06
6 4.609751 2.14703e-05 1.87069e-06 -0.000000 0.000000 0 0 1 1.47499e-06
========== Varying R2 net upward from 6.084 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.753210 0.0596636 0.0596636 0.004243 0.019324 1 0 1 0.000358423
2 1.518555 0.0845453 0.0248817 -0.000791 0.002156 1 0 1 0.000119474
3 2.285424 0.103567 0.0190212 -0.000535 0.000887 1 0 1 3.98247e-05
4 3.052810 0.119555 0.0159886 -0.000338 0.000568 1 0 1 1.32749e-05
5 3.820457 0.133608 0.0140525 -0.000218 0.000427 1 0 1 4.42497e-06
6 4.588261 0.146286 0.0126785 -0.000141 0.000347 1 0 1 1.47499e-06
========== Varying R2 net downward from 6.084 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.740402 0.0595483 0.0595483 -0.021366 0.006522 1 0 1 0.000358423
2 1.503313 0.0850441 0.0254958 -0.005295 0.000086 1 0 1 0.000119474
3 2.268190 0.104634 0.0195899 -0.003403 0.000012 1 0 1 3.98247e-05
4 3.033916 0.121179 0.0165447 -0.002561 0.000005 1 0 1 1.32749e-05
5 3.800120 0.135782 0.014603 -0.002084 0.000003 1 0 1 4.42497e-06
6 4.566633 0.149007 0.0132258 -0.001776 0.000003 1 0 1 1.47499e-06
========== Varying R2 exch upward from 6.621 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.667767 0.290301 0.290301 -0.079498 0.021026 1 0 1 0.000358423
2 1.410871 0.42865 0.138349 -0.023158 0.002030 1 0 1 0.000119474
3 2.162448 0.537269 0.10862 -0.015854 0.000861 1 0 1 3.98247e-05
4 2.917636 0.630696 0.0934262 -0.012564 0.000539 1 0 1 1.32749e-05
5 3.674891 0.714516 0.08382 -0.010646 0.000392 1 0 1 4.42497e-06
6 4.433503 0.791576 0.0770601 -0.009371 0.000308 1 0 1 1.47499e-06
========== Varying R2 exch downward from 6.621 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.770310 0.289411 0.289411 0.008753 0.006733 1 0 1 0.000358423
2 1.542661 0.403807 0.114396 0.004146 0.000088 1 0 1 0.000119474
3 2.314956 0.48932 0.0855128 0.004042 0.000038 1 0 1 3.98247e-05
4 3.087059 0.559935 0.0706148 0.003839 0.000028 1 0 1 1.32749e-05
5 3.858960 0.621063 0.0611278 0.003633 0.000024 1 0 1 4.42497e-06
========== Varying R4 upward from 1.958 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.740402 0.0297741 0.0297741 -0.021366 0.006522 1 0 1 0.000358423
2 1.503313 0.0425221 0.0127479 -0.005295 0.000086 1 0 1 0.000119474
3 2.268190 0.052317 0.00979493 -0.003403 0.000012 1 0 1 3.98247e-05
4 3.033916 0.0605893 0.00827235 -0.002561 0.000005 1 0 1 1.32749e-05
5 3.800120 0.0678908 0.00730148 -0.002084 0.000003 1 0 1 4.42497e-06
6 4.566633 0.0745037 0.0066129 -0.001776 0.000003 1 0 1 1.47499e-06
========== Varying R4 downward from 1.958 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.753210 0.0298318 0.0298318 0.004243 0.019324 1 0 1 0.000358423
2 1.518555 0.0422727 0.0124409 -0.000791 0.002156 1 0 1 0.000119474
3 2.285424 0.0517833 0.00951062 -0.000535 0.000887 1 0 1 3.98247e-05
4 3.052810 0.0597776 0.00799428 -0.000338 0.000568 1 0 1 1.32749e-05
5 3.820457 0.0668038 0.00702625 -0.000218 0.000427 1 0 1 4.42497e-06
6 4.588261 0.0731431 0.00633925 -0.000141 0.000347 1 0 1 1.47499e-06
========== Varying R5 upward from 8.042 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.753210 0.0298318 0.0298318 0.004243 0.019324 1 0 1 0.000358423
2 1.518555 0.0422727 0.0124409 -0.000791 0.002156 1 0 1 0.000119474
3 2.285424 0.0517833 0.00951062 -0.000535 0.000887 1 0 1 3.98247e-05
4 3.052810 0.0597776 0.00799428 -0.000338 0.000568 1 0 1 1.32749e-05
5 3.820457 0.0668038 0.00702625 -0.000218 0.000427 1 0 1 4.42497e-06
6 4.588261 0.0731431 0.00633925 -0.000141 0.000347 1 0 1 1.47499e-06
========== Varying R5 downward from 8.042 ==========
Delta Delta Predictor Corrector Corrector Failed
Iteration residual parameter Step-size error adjustment iterations steps h/hopt Lambda
1 0.740402 0.0297741 0.0297741 -0.021366 0.006522 1 0 1 0.000358423
2 1.503313 0.0425221 0.0127479 -0.005295 0.000086 1 0 1 0.000119474
3 2.268190 0.052317 0.00979493 -0.003403 0.000012 1 0 1 3.98247e-05
4 3.033916 0.0605893 0.00827235 -0.002561 0.000005 1 0 1 1.32749e-05
5 3.800120 0.0678908 0.00730148 -0.002084 0.000003 1 0 1 4.42497e-06
6 4.566633 0.0745037 0.0066129 -0.001776 0.000003 1 0 1 1.47499e-06
Continuation completed in 1.5000 seconds.
Preprocessing time: 0.3300 s
Computation time: 1.1700 s
Warning: Network is ill-conditioned.
Simulation completed in 0.0700 seconds.
Preprocessing time: 0.0600 s
Computation time: 0.0100 s
Postprocessing time: 0.0000 s
--- 23.097320795059204 seconds -
We can now read the results from INCA using the INCAResults object.
[21]:
res = incawrapper.INCAResults(output_file)
Typically, the first thing we want to inspect is the fitted fluxes. They can be displayed as a Pandas dataframe.
[22]:
res.fitdata.fitted_parameters
[22]:
| type | id | eqn | val | std | lb | ub | unit | free | alf | chi2s | cont | cor | cov | vals | base | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Net flux | R1 | A -> B | 10.000000 | 0.000010 | 9.99998 | 10.00002 | [] | 0 | 0.05 | [141.10066093755805, 140.33236917370687, 139.5... | 0 | [1.0, 8.94704371251548e-05, 2.002130779029279e... | [9.999999786569624e-11, 6.084217157864558e-11,... | [9.999978529670276, 9.999980400360155, 9.99998... | {'id': []} |
| 1 | Net flux | R2 net | B <-> D | 6.084230 | 0.068003 | 5.947702 | 6.218202 | [] | 1 | 0.05 | [141.0575435245828, 140.2910305858517, 139.524... | 0 | [8.947043202321765e-05, 1.0, 0.928424443085881... | [6.084216810919862e-11, 0.004624345816403641, ... | [5.935222845778996, 5.948448640360109, 5.96305... | {'id': []} |
| 2 | Exch flux | R2 exch | B <-> D | 6.620853 | 0.330687 | 6.001071 | 7.352861 | [] | 1 | 0.05 | [140.34987009287295, 139.5779692296173, 138.80... | 0 | [2.0021303593639435e-05, 0.9284244520811099, 0... | [6.620785170108334e-11, 0.020877998964188516, ... | [5.999789999353027, 6.060917783460922, 6.13153... | {'id': []} |
| 3 | Net flux | R3 | B -> C + E | 1.957885 | 0.034001 | 1.890899 | 2.026149 | [] | 1 | 0.05 | [141.07917098749863, 140.31136692375648, 139.5... | 0 | [5.758289131530709e-05, -0.9999999891876596, -... | [1.957891487824881e-11, -0.0023121728777807347... | [1.8847417919539917, 1.8910810441689598, 1.898... | {'id': []} |
| 4 | Net flux | R4 | B + C -> D + E + E | 1.957885 | 0.034001 | 1.890899 | 2.026149 | [] | 0 | 0.05 | [141.07917098749863, 140.31136692375648, 139.5... | 0 | [5.758289131530709e-05, -0.9999999891876596, -... | [1.957891487824881e-11, -0.0023121728777807347... | [1.8847417919539917, 1.8910810441689598, 1.898... | {'id': []} |
| 5 | Net flux | R5 | D -> F | 8.042115 | 0.034001 | 7.973851 | 8.109101 | [] | 0 | 0.05 | [141.05754359730764, 140.2910306695598, 139.52... | 0 | [0.00023652374955743415, 0.9999999891876609, 0... | [8.042108298744743e-11, 0.0023121729386229063,... | [7.967611420633636, 7.974224317908821, 7.98152... | {'id': []} |
| 6 | Norm | exp1 F1 exp1_F1_0_0_1 norm | [] | 0.985739 | 0.016453 | [] | [] | [] | 1 | 0.05 | [] | 0 | [-7.775915103023788e-11, 0.5822932413163197, 0... | [-1.2793585635328952e-17, 0.000651489883093479... | [] | {'id': []} |
[23]:
res.fitdata.get_goodness_of_fit()
Fit accepted: False
Confidence level: 0.05
Chi-square value (SSR): 136.49091035304212
Expected chi-square range: [9.82069117e-04 5.02388619e+00]
[24]:
import incawrapper.visualization as incaviz
incaviz.plot_norm_prob(res)
[24]:
(<Figure size 640x480 with 1 Axes>,
<Axes: title={'center': 'Probability Plot'}, xlabel='Theoretical quantiles', ylabel='Ordered Values'>)
2.1. References¶
[1] M. R. Antoniewicz, J. K. Kelleher, and G. Stephanopoulos, “Determination of confidence intervals of metabolic fluxes estimated from stable isotope measurements,” Metabolic Engineering, vol. 8, no. 4, pp. 324–337, Jul. 2006, doi: 10.1016/j.ymben.2006.01.004.
[2] M. R. Antoniewicz, J. K. Kelleher, and G. Stephanopoulos, “Elementary metabolite units (EMU): A novel framework for modeling isotopic distributions,” Metabolic Engineering, vol. 9, no. 1, pp. 68–86, Jan. 2007, doi: 10.1016/j.ymben.2006.09.001.