{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# INCAWrapper validation case medium size model\n",
"In this notebook, we use the INCAWrapper to fit a simulated data. This notebook serves as an integration test the INCAWrapper that can be run when changes are made to the codebase to ensure that the INCAWrapper performs consistently. Furthermore, this case act as a validation case the show that the INCAWrapper actually works and that it is capable of estimating correct flux distribution from simulated data.\n",
"\n",
"This notebook is not meant as a tutorial and therefore code description is a bit more sparse. For a proper tutorial see the other examples at https://incawrapper.readthedocs.io/en/latest/examples/index.html.\n",
"\n",
"The model we use is from Alagesan, S., Minton, N.P. & Malys, N. 13C-assisted metabolic flux analysis to investigate heterotrophic and mixotrophic metabolism in Cupriavidus necator H16. Metabolomics 14, 9 (2018). https://doi.org/10.1007/s11306-017-1302-z, which is also used for one of the tutorials.\n",
"\n",
"The simulation mimics a two parallel experiment where C. necator is grown with different labelled fructose and labelled glycerol. We simulated MS measurements of the amino acids and a few exchanges fluxes are measured. To increase the information about the systems we simulate measurements of CO2 exchange flux. To do this we added one additional reaction to the original model, i.e. CO2 -> CO2.ext. For more details about the simulation see the [simulation script](https://github.com/biosustain/incawrapper/blob/main/docs/examples/Literature%20data/Cupriavidus%20necator%2520%20Alagesan%202017/c_necator_simulation.py)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"import dotenv\n",
"import ast\n",
"import pandera as pa\n",
"import incawrapper\n",
"from incawrapper import utils\n",
"from incawrapper import visualization\n",
"import pathlib\n",
"import matplotlib.pyplot as plt\n",
"import pytest"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# import environment variables\n",
"INCA_base_directory = dotenv.get_key(dotenv.find_dotenv(), \"INCA_base_directory\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# set up path to data\n",
"working_dir = pathlib.Path(dotenv.find_dotenv()).parent\n",
"data_directory = working_dir / 'docs' / 'examples' / 'Literature data' / 'Cupriavidus necator Alagesan 2017' / 'simulated_data'\n",
"results_file = data_directory / 'fit_simulated_fructose.mat'"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Reading the reactions, tracers, and simulated measurements\n",
"rxn = pd.read_csv(data_directory / 'reactions_processed.csv')\n",
"tracers = pd.read_csv(\n",
" data_directory / 'tracer_info.csv', \n",
" converters={\"atom_ids\": ast.literal_eval, \"atom_mdv\": ast.literal_eval}\n",
")\n",
"flux_measurements = pd.read_csv(data_directory / 'flux_measurements_no_noise.csv')\n",
"mdv_measurements = pd.read_csv(\n",
" data_directory / 'mdv_no_noise.csv',\n",
" converters={\"labelled_atom_ids\": ast.literal_eval}\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Looking at the traces dataframe, we can see the design of the two simulated parallel experiments"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" experiment_id \n",
" met_id \n",
" tracer_id \n",
" atom_ids \n",
" atom_mdv \n",
" enrichment \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" simulation1 \n",
" FRU.ext \n",
" D-[1-13C]fructose \n",
" [1] \n",
" [0.0, 1.0] \n",
" 1 \n",
" \n",
" \n",
" 1 \n",
" simulation1 \n",
" GLY.ext \n",
" [1,2-13C]glycerol \n",
" [1, 2] \n",
" [0.0, 1.0] \n",
" 1 \n",
" \n",
" \n",
" 2 \n",
" simulation2 \n",
" FRU.ext \n",
" [1,6-13C]fructose \n",
" [1, 6] \n",
" [0.0, 1.0] \n",
" 1 \n",
" \n",
" \n",
" 3 \n",
" simulation2 \n",
" GLY.ext \n",
" [1,2-13C]glycerol \n",
" [1, 2] \n",
" [0.0, 1.0] \n",
" 1 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"yaxis_range = [-250, 400]\n",
"plot_df = true_and_fitted_fluxes.copy()\n",
"# replace inf with large number\n",
"plot_df['lb'].replace(np.inf, yaxis_range[0], inplace=True)\n",
"plot_df['ub'].replace(np.inf, yaxis_range[1], inplace=True)\n",
"fig, ax = plt.subplots(figsize=(10, 5))\n",
"errbars = plot_df[['lb', 'ub']].subtract(plot_df['val_fitted'], axis=0).abs().T\n",
"ax.scatter(x=plot_df['id'], y=plot_df['val_fitted'], color='black', label='fitted value')\n",
"ax.errorbar(x=plot_df['id'], y=plot_df['val_fitted'], yerr=errbars, color='black', fmt='none', label='95% CI of fitted value')\n",
"ax.scatter(x=plot_df['id'], y=plot_df['val_true'], color='red', label='literature value', s=8)\n",
"# rotate x-axis labels\n",
"ax.tick_params(axis='x', labelrotation=90)\n",
"ax.legend()\n",
"ax.set_ylim(yaxis_range)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that most of the true flux value lies close to the estimated flux value. However, some of the flux estimates have large uncertainties. Therefor, we investigate how many fluxes are \"well determined\"."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def well_determined_fluxes(val, lb, ub, cutoff):\n",
" '''Finds well defined fluxes, these are fluxes where the relative uncertainty span is less than the cutoff.'''\n",
" if val == 0:\n",
" return (ub - lb) < cutoff\n",
" relative_uncertainty_span = (ub - lb) / np.abs(val)\n",
" return relative_uncertainty_span < cutoff\n",
"\n",
"assert well_determined_fluxes(1, 0.9, 1.1, 0.21) == True\n",
"assert well_determined_fluxes(1, 0.9, 1.1, 0.1) == False\n",
"assert well_determined_fluxes(0, -0.49, 0.49, 1) == True"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.37333333333333335"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"true_and_fitted_fluxes['well_determined'] = true_and_fitted_fluxes.apply(\n",
" lambda x: well_determined_fluxes(x['val_true'], x['lb'], x['ub'], 1), axis=1\n",
")\n",
"true_and_fitted_fluxes['well_determined'].sum() / true_and_fitted_fluxes['well_determined'].size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that only 37% of the fluxes are well determined using a cutoff=1. This shows that the simulated experiments do not contain sufficient information about the flux distribution to gain an accurate estimate. This could likely be improve by introducing more simulated parallel experiments, using other labelled substrates. However, some fluxes, e.g. the exchange flux of reversible reactions, can be difficult to estimate at all.\n",
"\n",
"We will calculate the fraction of well determined net fluxes"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.45"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"net_fluxes = true_and_fitted_fluxes['type'] == 'Net flux'\n",
"true_and_fitted_fluxes.loc[net_fluxes, 'well_determined'].sum() / net_fluxes.sum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is still not super impressive, however, visual inspection of the plot shows that the flux estimate is in good agreement with the simulated flux distribution. We strongly believe that this performance is due lack of information in the simulated data, rather that errors in the INCAWrapper.\n",
"\n",
"Finally, we have a few tests which validates that the analysis was consistent with earlier runs for testing purposes."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"## Setup validation tests that test that the results are similar to the previous run\n",
"# Test that the fraction of correct fluxes similar to ealier run\n",
"assert fraction_correct_fluxes == 1.0\n",
"\n",
"# Test that the chi2 is similar to ealier run, this test may be fragile because the flux \n",
"# estimation is stochastic. Therefore, we only test that the chi2 very small (i.e. < 5).\n",
"assert res.fitdata.chi2 < 5\n",
"\n",
"# Test that the fraction of well determined fluxes is similar to ealier run\n",
"assert true_and_fitted_fluxes['well_determined'].sum() / true_and_fitted_fluxes['well_determined'].size == pytest.approx(0.373, 0.005)\n",
"assert true_and_fitted_fluxes.loc[net_fluxes, 'well_determined'].sum() / net_fluxes.sum() == pytest.approx(0.45, 0.005)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "incawrapper-dev",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}