{ "cells": [ { "cell_type": "markdown", "id": "f11cff02", "metadata": {}, "source": [ "# Libraries and useful functions" ] }, { "cell_type": "markdown", "id": "8cbb2a33", "metadata": {}, "source": [ "**Important**\n", "\n", "We are going to use two different models for small E. Coli: One FBA model (.mat file) and the thermo version we created in the previous excercise (.json)" ] }, { "cell_type": "code", "execution_count": 16, "id": "0ddc0742", "metadata": {}, "outputs": [], "source": [ "from cobra.io import load_matlab_model\n", "from gurobipy import GRB \n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import os\n", "from pytfa.io.json import load_json_model\n", "import pandas as pd\n", "import math" ] }, { "cell_type": "code", "execution_count": 2, "id": "b7d8186c", "metadata": {}, "outputs": [], "source": [ "models_path = './models/'" ] }, { "cell_type": "code", "execution_count": 3, "id": "2e426b1a", "metadata": {}, "outputs": [], "source": [ "def setSolver(model, solver_name):\n", " \"\"\"\n", " Set the solver for a given COBRA model.\n", "\n", " Parameters:\n", " model (cobra.Model): The COBRA model to set the solver for.\n", " solver_name (str): The name of the solver to use ('cplex' or 'gurobi').\n", "\n", " Returns:\n", " None: The function modifies the model in place.\n", " \"\"\"\n", " model.solver = solver_name\n", "\n", "def printSolverParams(model, solver):\n", " if solver == \"gurobi\":\n", " print(model.solver.problem.Params)\n", " elif solver == \"cplex\":\n", " print(\"feasibility:\" + str(model.solver.configuration.tolerances.feasibility))\n", " print(\"integrality:\"+ str(model.solver.configuration.tolerances.integrality))\n", " print(\"optimality\" + str(model.solver.configuration.tolerances.optimality))\n", " elif solver == \"glpk\":\n", " print(\"feasibility:\" + str(model.solver.configuration.tolerances.feasibility))\n", " print(\"integrality:\"+ str(model.solver.configuration.tolerances.integrality))\n", " else:\n", " raise ValueError(\"Cannot modify tolerance params for solver \" + solver)\n", " pass\n", "\n", "def setSolverParams(model, solver, FeasibilityTol =1e-9, IntFeasTol=1e-9, OptimalityTol=1e-9):\n", " \"\"\"\n", " Set solver parameters for a given COBRA model.\n", " Parameters:\n", " model (cobra.Model): The COBRA model to set the solver parameters for.\n", " solver (str): The name of the solver to use ('cplex' or 'gurobi').\n", " FeasibilityTol (float): The feasibility tolerance to set.\n", " IntFeasTol (float): The integrality tolerance to set.\n", " OptimalityTol (float): The optimality tolerance to set.\n", "\n", " Returns:\n", " None: The function modifies the model in place.\n", " \"\"\"\n", " if solver == \"gurobi\":\n", " model.solver.problem.Params.FeasibilityTol = FeasibilityTol\n", " model.solver.problem.Params.IntFeasTol = IntFeasTol\n", " model.solver.problem.Params.OptimalityTol = OptimalityTol\n", " model.solver.problem.Params.Presolve = 1\n", " model.solver.problem.Params.Timelimit = 3600\n", " model.solver.problem.setParam(GRB.Param.NumericFocus, 1)\n", " model.solver.problem.setParam(GRB.Param.ScaleFlag, 2)\n", " elif solver == \"cplex\":\n", " model.solver.configuration.tolerances.feasibility = FeasibilityTol\n", " model.solver.configuration.tolerances.integrality = IntFeasTol\n", " model.solver.configuration.tolerances.optimality = OptimalityTol\n", " model.solver.problem.parameters.read_scale = -1\n", " model.solver.problem.parameters.emphasis.numerical = 1\n", " model.solver.problem.parameters.timelimit = 3600\n", " model.solver.configuration.presolve = True\n", " #model_engro2_thermo.solver.problem.parameters.mip.display.set(3)\n", " #model_engro2_thermo.solver.problem.set_results_stream()\n", " elif solver == \"glpk\":\n", " model.solver.configuration.tolerances.feasibility = FeasibilityTol\n", " model.solver.configuration.tolerances.integrality = IntFeasTol\n", " else:\n", " raise ValueError(\"Cannot modify tolerance params for solver \" + solver)\n", " \n", " printSolverParams(model, solver)\n", " pass" ] }, { "cell_type": "markdown", "id": "f1e1bbb6", "metadata": {}, "source": [ "# Part 1" ] }, { "cell_type": "code", "execution_count": 11, "id": "f0409a0f", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "This model seems to have metCharge instead of metCharges field. Will use metCharge for what metCharges represents.\n", "No defined compartments in model model_red. Compartments will be deduced heuristically using regular expressions.\n", "Using regular expression found the following compartments:c, e, p\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "feasibility:1e-09\n", "integrality:1e-09\n", "optimality1e-09\n" ] } ], "source": [ "model = load_matlab_model(os.path.join(models_path, 'small_ecoli.mat'))\n", "setSolver(model, 'cplex')\n", "setSolverParams(model, \"cplex\", FeasibilityTol =1e-9, IntFeasTol=1e-9, OptimalityTol=1e-9)" ] }, { "cell_type": "code", "execution_count": 12, "id": "4ca5d8b5", "metadata": {}, "outputs": [ { "data": { "text/html": [ "Optimal solution with objective value 0.811
\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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
fluxesreduced_costs
DM_4CRSOL0.0001810.000000e+00
DM_5DRIB0.0001870.000000e+00
DM_AMOB0.0000020.000000e+00
DM_MTHTHF0.0010870.000000e+00
Ec_biomass_iJO1366_WT_53p95M0.8109623.569410e-15
.........
LMPD_250_trp-L_c0.044793-2.775558e-17
LMPD_251_tyr-L_c0.1086630.000000e+00
LMPD_252_udcpdp_c0.000045-6.661338e-16
LMPD_253_utp_c0.1136170.000000e+00
LMPD_254_val-L_c0.3334550.000000e+00
\n", "

599 rows × 2 columns

\n", "
" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.optimize()" ] }, { "cell_type": "code", "execution_count": 7, "id": "4dc55081", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "DM_4CRSOL\n", "DM_5DRIB\n", "DM_AMOB\n", "DM_MTHTHF\n", "DM_5mtr_e\n", "DM_ac_e\n", "DM_acald_e\n", "DM_acser_e\n", "DM_akg_e\n", "DM_ala-D_e\n", "DM_ala-L_e\n", "DM_alaala_e\n", "DM_arg-L_e\n", "DM_asn-L_e\n", "DM_asp-L_e\n", "DM_ca2_e\n", "DM_cbl1_e\n", "DM_cit_e\n", "DM_cl_e\n", "DM_co2_e\n", "DM_cobalt2_e\n", "DM_colipa_e\n", "DM_cu2_e\n", "DM_cys-L_e\n", "DM_dha_e\n", "DM_enter_e\n", "DM_etoh_e\n", "DM_fe2_e\n", "DM_fe3_e\n", "DM_for_e\n", "DM_glc_e\n", "DM_glcn_e\n", "DM_glu-L_e\n", "DM_gly_e\n", "DM_glyc_e\n", "DM_glyc-R_e\n", "DM_glycogenn1_c\n", "DM_glyclt_e\n", "DM_gthrd_e\n", "DM_h_e\n", "DM_h2o_e\n", "DM_his-L_e\n", "DM_ile-L_e\n", "DM_k_e\n", "DM_lac-D_e\n", "DM_leu-L_e\n", "DM_lys-L_e\n", "DM_mal-L_e\n", "DM_meoh_e\n", "DM_mg2_e\n", "DM_mn2_e\n", "DM_mobd_e\n", "DM_nh4_e\n", "DM_ni2_e\n", "DM_o2_e\n", "DM_phe-L_e\n", "DM_pheme_e\n", "DM_pi_e\n", "DM_pro-L_e\n", "DM_pyr_e\n", "DM_ser-L_e\n", "DM_so4_e\n", "DM_spmd_e\n", "DM_succ_e\n", "DM_thr-L_e\n", "DM_trp-L_e\n", "DM_tyr-L_e\n", "DM_val-L_e\n", "DM_zn2_e\n" ] } ], "source": [ "for rxn in model.reactions:\n", " if (rxn.id.startswith(\"DM_\")):\n", " print(rxn.id)" ] }, { "cell_type": "code", "execution_count": 4, "id": "6bd40687", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2025-09-26 13:50:14,298 - thermomodel_tutorial_ecoli - INFO - # Model initialized with units kcal/mol and temperature 298.15 K\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "feasibility:1e-09\n", "integrality:1e-09\n", "optimality1e-09\n" ] } ], "source": [ "tmodel = load_json_model(os.path.join(models_path, 'small_ecoli_thermo.json'))\n", "setSolver(tmodel, 'cplex')\n", "setSolverParams(tmodel, \"cplex\", FeasibilityTol =1e-9, IntFeasTol=1e-9, OptimalityTol=1e-9)" ] }, { "cell_type": "code", "execution_count": 5, "id": "ae805c3c", "metadata": {}, "outputs": [ { "data": { "text/html": [ "Optimal solution with objective value 0.811
\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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
fluxesreduced_costs
DM_4CRSOL0.000181None
DM_5DRIB0.000187None
DM_AMOB0.000002None
DM_MTHTHF0.001087None
Ec_biomass_iJO1366_WT_53p95M0.810997None
.........
LMPD_250_trp-L_c0.044795None
LMPD_251_tyr-L_c0.108668None
LMPD_252_udcpdp_c0.000045None
LMPD_253_utp_c0.113622None
LMPD_254_val-L_c0.333469None
\n", "

599 rows × 2 columns

\n", "
" ], "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tmodel.optimize()" ] }, { "cell_type": "code", "execution_count": 9, "id": "309495c2", "metadata": {}, "outputs": [], "source": [ "id_biomass = \"Ec_biomass_iJO1366_WT_53p95M\"\n", "\n", "dict_carbon_sources = {\n", " 'DM_glc_e':'D-glucose',\n", " 'DM_lac-D_e':'D-lactate',\n", " 'DM_ac_e':'Acetate', \n", " 'DM_etoh_e':'Ethanol'\n", "}" ] }, { "cell_type": "code", "execution_count": 13, "id": "a49718f5", "metadata": {}, "outputs": [], "source": [ "# We check if default bounds are the same in both models\n", "for rxn in model.reactions:\n", " if(rxn.lower_bound != tmodel.reactions.get_by_id(rxn.id).lower_bound or rxn.upper_bound != tmodel.reactions.get_by_id(rxn.id).upper_bound):\n", " print(\"Different bounds for reaction \" + rxn.id + \": FBA model [\" + str(rxn.lower_bound) + \",\" + str(rxn.upper_bound) + \"] - Thermo model [\" + str(tmodel.reactions.get_by_id(rxn.id).lower_bound) + \",\" + str(tmodel.reactions.get_by_id(rxn.id).upper_bound) + \"]\")" ] }, { "cell_type": "code", "execution_count": 14, "id": "26f4e6e0", "metadata": {}, "outputs": [], "source": [ "#let's save original bounds\n", "original_bounds = {}\n", "for rxn in model.reactions:\n", " original_bounds[rxn.id] = (rxn.lower_bound, rxn.upper_bound)" ] }, { "cell_type": "code", "execution_count": 15, "id": "8fe13eae", "metadata": {}, "outputs": [], "source": [ "#let's fix objective function\n", "model.objective = id_biomass\n", "tmodel.objective = id_biomass\n", "#maximize biomass\n", "model.objective_direction = 'max'\n", "tmodel.objective_direction = 'max'" ] }, { "cell_type": "code", "execution_count": 17, "id": "0ec68338", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
FBA growth rate (aerobic)FBA growth rate (anaerobic)Thermo growth rate (aerobic)Thermo growth rate (anaerobic)
DM_glc_eNaNNaNNaNNaN
DM_lac-D_eNaNNaNNaNNaN
DM_ac_eNaNNaNNaNNaN
DM_etoh_eNaNNaNNaNNaN
\n", "
" ], "text/plain": [ " FBA growth rate (aerobic) FBA growth rate (anaerobic) \\\n", "DM_glc_e NaN NaN \n", "DM_lac-D_e NaN NaN \n", "DM_ac_e NaN NaN \n", "DM_etoh_e NaN NaN \n", "\n", " Thermo growth rate (aerobic) Thermo growth rate (anaerobic) \n", "DM_glc_e NaN NaN \n", "DM_lac-D_e NaN NaN \n", "DM_ac_e NaN NaN \n", "DM_etoh_e NaN NaN " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#let's create our dataframe to store results\n", "df_res = pd.DataFrame(columns=['FBA growth rate (aerobic)', 'FBA growth rate (anaerobic)', 'Thermo growth rate (aerobic)', 'Thermo growth rate (anaerobic)'], index=dict_carbon_sources.keys())\n", "df_res" ] }, { "cell_type": "code", "execution_count": 19, "id": "7b021640", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Working with model: FBA\n", "\n", "Condition: aerobic\n", "\n", "Condition: anaerobic\n", "\n", "Working with model: Thermo\n", "\n", "Condition: aerobic\n", "\n", "Condition: anaerobic\n" ] } ], "source": [ "# First we work with the FBA model\n", "for model_selected in [\"FBA\", \"Thermo\"]:\n", "\n", " if(model_selected == \"Thermo\"):\n", " selected_model = tmodel\n", " else:\n", " selected_model = model\n", " print(\"\\nWorking with model: \" + model_selected)\n", "\n", " for condition in [\"aerobic\", \"anaerobic\"]:\n", "\n", " print(\"\\nCondition: \" + condition)\n", "\n", " if condition == \"aerobic\":\n", " selected_model.reactions.get_by_id(\"DM_o2_e\").lower_bound = -20 # We allow unlimited oxygen uptake\n", " else:\n", " selected_model.reactions.get_by_id(\"DM_o2_e\").lower_bound = 0 # We do not allow oxygen uptake\n", "\n", "\n", " for key, val in dict_carbon_sources.items():\n", "\n", " #allow uptake of one carbon source\n", " selected_model.reactions.get_by_id(key).lower_bound = -10\n", "\n", " #fix all others to 0\n", " for key2, val2 in dict_carbon_sources.items():\n", " if(key2 != key):\n", " selected_model.reactions.get_by_id(key2).lower_bound = 0\n", " \n", " #optimize\n", " sol = selected_model.optimize()\n", " if(sol.status != 'optimal'):\n", " print(\"No optimal solution found for carbon source \" + val + \" (\" + key + \")\")\n", " sol_f = math.nan\n", " else:\n", " sol_f = sol.objective_value\n", "\n", " df_res.loc[key, model_selected + \" growth rate (\" + condition + \")\"] = sol_f\n", "\n", " selected_model = None\n", "\n", "#restore original bounds\n", "for rxn in model.reactions:\n", " rxn.lower_bound = original_bounds[rxn.id][0]\n", " rxn.upper_bound = original_bounds[rxn.id][1]\n", "for rxn in tmodel.reactions:\n", " rxn.lower_bound = original_bounds[rxn.id][0]\n", " rxn.upper_bound = original_bounds[rxn.id][1]\n", " \n", " " ] }, { "cell_type": "code", "execution_count": 20, "id": "d1c7c3c9", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
FBA growth rate (aerobic)FBA growth rate (anaerobic)Thermo growth rate (aerobic)Thermo growth rate (anaerobic)
DM_glc_e0.9938260.2731290.9938690.273155
DM_lac-D_e0.4392140.0247510.4392320.024753
DM_ac_e0.2594010.00.2594110.0
DM_etoh_e0.4323350.00.4323520.0
\n", "
" ], "text/plain": [ " FBA growth rate (aerobic) FBA growth rate (anaerobic) \\\n", "DM_glc_e 0.993826 0.273129 \n", "DM_lac-D_e 0.439214 0.024751 \n", "DM_ac_e 0.259401 0.0 \n", "DM_etoh_e 0.432335 0.0 \n", "\n", " Thermo growth rate (aerobic) Thermo growth rate (anaerobic) \n", "DM_glc_e 0.993869 0.273155 \n", "DM_lac-D_e 0.439232 0.024753 \n", "DM_ac_e 0.259411 0.0 \n", "DM_etoh_e 0.432352 0.0 " ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_res" ] }, { "cell_type": "markdown", "id": "1d5d9898", "metadata": {}, "source": [ "# Part 2" ] }, { "cell_type": "code", "execution_count": 23, "id": "f734edb1", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\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", " \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", " \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", " \n", " \n", " \n", " \n", "
MetaboliteConcentration in mmolStandard deviation
model ID
akg_c2-Oxoglutarate0.030100NaN
3pg_c3-Phosphoglycerate0.5060000.190000
ade_cAdenine0.1580000.064600
adn_cAdenosine0.0837000.032900
adp_cADP0.5620000.255000
............
gln-L_cGln0.6414500.387912
gtp_cGTP0.6272290.306649
pep_cPhosphoenolpyruvate0.1200000.047200
ptrc_c1.4-Butanediamine0.1135800.033971
pydx5p_cPyridoxal 5-phosphate0.244000NaN
\n", "

72 rows × 3 columns

\n", "
" ], "text/plain": [ " Metabolite Concentration in mmol Standard deviation\n", "model ID \n", "akg_c 2-Oxoglutarate 0.030100 NaN\n", "3pg_c 3-Phosphoglycerate 0.506000 0.190000\n", "ade_c Adenine 0.158000 0.064600\n", "adn_c Adenosine 0.083700 0.032900\n", "adp_c ADP 0.562000 0.255000\n", "... ... ... ...\n", "gln-L_c Gln 0.641450 0.387912\n", "gtp_c GTP 0.627229 0.306649\n", "pep_c Phosphoenolpyruvate 0.120000 0.047200\n", "ptrc_c 1.4-Butanediamine 0.113580 0.033971\n", "pydx5p_c Pyridoxal 5-phosphate 0.244000 NaN\n", "\n", "[72 rows x 3 columns]" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "metabolomics = pd.read_csv('metabolomics_data.csv', index_col=0, sep=';')\n", "metabolomics" ] }, { "cell_type": "code", "execution_count": 50, "id": "b245ba80", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "There are duplicated entries in the metabolomics file\n", " Metabolite Concentration in mmol Standard deviation\n", "model ID \n", "datp_c dADP 0.0352 0.0205\n", "datp_c dATP 0.0403 0.0284\n" ] } ], "source": [ "#check for duplicate entries\n", "if(metabolomics.index.duplicated().any()):\n", " print(\"There are duplicated entries in the metabolomics file\")\n", " print(metabolomics[metabolomics.index.duplicated(keep=False)])" ] }, { "cell_type": "code", "execution_count": 51, "id": "2c41ea23", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\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", " \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", " \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", " \n", " \n", " \n", " \n", "
MetaboliteConcentration in mmolStandard deviation
model ID
akg_c2-Oxoglutarate0.030100NaN
3pg_c3-Phosphoglycerate0.5060000.190000
ade_cAdenine0.1580000.064600
adn_cAdenosine0.0837000.032900
adp_cADP0.5620000.255000
............
gln-L_cGln0.6414500.387912
gtp_cGTP0.6272290.306649
pep_cPhosphoenolpyruvate0.1200000.047200
ptrc_c1.4-Butanediamine0.1135800.033971
pydx5p_cPyridoxal 5-phosphate0.244000NaN
\n", "

71 rows × 3 columns

\n", "
" ], "text/plain": [ " Metabolite Concentration in mmol Standard deviation\n", "model ID \n", "akg_c 2-Oxoglutarate 0.030100 NaN\n", "3pg_c 3-Phosphoglycerate 0.506000 0.190000\n", "ade_c Adenine 0.158000 0.064600\n", "adn_c Adenosine 0.083700 0.032900\n", "adp_c ADP 0.562000 0.255000\n", "... ... ... ...\n", "gln-L_c Gln 0.641450 0.387912\n", "gtp_c GTP 0.627229 0.306649\n", "pep_c Phosphoenolpyruvate 0.120000 0.047200\n", "ptrc_c 1.4-Butanediamine 0.113580 0.033971\n", "pydx5p_c Pyridoxal 5-phosphate 0.244000 NaN\n", "\n", "[71 rows x 3 columns]" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#let's keep the first entry for each duplicated metabolite\n", "metabolomics = metabolomics[~metabolomics.index.duplicated(keep='first')] \n", "metabolomics" ] }, { "cell_type": "code", "execution_count": 37, "id": "2687859b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd.isna(metabolomics.loc[\"akg_c\"][\"Standard deviation\"])" ] }, { "cell_type": "code", "execution_count": 27, "id": "a50842dd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of metabolites found in model: 60 out of 72\n" ] } ], "source": [ "mets_file = metabolomics.index.tolist()\n", "mets_model_metabolomics = [met for met in mets_file if met in tmodel.metabolites]\n", "print(\"Number of metabolites found in model: \" + str(len(mets_model_metabolomics)) + \" out of \" + str(len(mets_file)))" ] }, { "cell_type": "code", "execution_count": 28, "id": "d6391b9c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['akg_c',\n", " '3pg_c',\n", " 'adp_c',\n", " 'ala-L_c',\n", " 'arg-L_c',\n", " 'asn-L_c',\n", " 'asp-L_c',\n", " 'atp_c',\n", " 'spmd_c',\n", " 'ctp_c',\n", " 'datp_c',\n", " 'datp_c',\n", " 'dctp_c',\n", " 'dhap_c',\n", " 'fdp_c',\n", " 'f6p_c',\n", " 'fum_c',\n", " 'g1p_c',\n", " 'g6p_c',\n", " 'glu-L_c',\n", " 'gthrd_c',\n", " 'gly_c',\n", " 'his-L_c',\n", " 'ile-L_c',\n", " 'leu-L_c',\n", " 'lys-L_c',\n", " 'mal-L_c',\n", " 'met-L_c',\n", " 'nad_c',\n", " 'nadp_c',\n", " 'nadph_c',\n", " 'phe-L_c',\n", " 'pro-L_c',\n", " 'pyr_c',\n", " 'r5p_c',\n", " 'ru5p-D_c',\n", " 'amet_c',\n", " 's7p_c',\n", " 'ser-L_c',\n", " 'succ_c',\n", " 'thmpp_c',\n", " 'thr-L_c',\n", " 'trp-L_c',\n", " 'tyr-L_c',\n", " 'udcpdp_c',\n", " 'utp_c',\n", " 'val-L_c',\n", " 'glc-D_e',\n", " 'etoh_e',\n", " 'ac_e',\n", " 'lac-D_e',\n", " 'succ_e',\n", " 'pyr_e',\n", " 'dttp_c',\n", " 'fad_c',\n", " 'gln-L_c',\n", " 'gtp_c',\n", " 'pep_c',\n", " 'ptrc_c',\n", " 'pydx5p_c']" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mets_model_metabolomics" ] }, { "cell_type": "code", "execution_count": 29, "id": "c05ce4b4", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
FBA growth rate (aerobic)FBA growth rate (anaerobic)Thermo growth rate (aerobic)Thermo growth rate (anaerobic)
DM_glc_eNaNNaNNaNNaN
DM_lac-D_eNaNNaNNaNNaN
DM_ac_eNaNNaNNaNNaN
DM_etoh_eNaNNaNNaNNaN
\n", "
" ], "text/plain": [ " FBA growth rate (aerobic) FBA growth rate (anaerobic) \\\n", "DM_glc_e NaN NaN \n", "DM_lac-D_e NaN NaN \n", "DM_ac_e NaN NaN \n", "DM_etoh_e NaN NaN \n", "\n", " Thermo growth rate (aerobic) Thermo growth rate (anaerobic) \n", "DM_glc_e NaN NaN \n", "DM_lac-D_e NaN NaN \n", "DM_ac_e NaN NaN \n", "DM_etoh_e NaN NaN " ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#let's create our dataframe to store results\n", "df_res2 = pd.DataFrame(columns=['FBA growth rate (aerobic)', 'FBA growth rate (anaerobic)', 'Thermo growth rate (aerobic)', 'Thermo growth rate (anaerobic)'], index=dict_carbon_sources.keys())\n", "df_res2" ] }, { "cell_type": "code", "execution_count": 32, "id": "7f292dc0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LC_10fthf_c\n", "LC_13dpg_c\n", "LC_2ddg6p_c\n", "LC_2dmmq8_c\n", "LC_2dmmql8_c\n", "LC_2fe2s_c\n", "LC_2h3oppan_c\n", "LC_2pg_c\n", "LC_3pg_c\n", "LC_4crsol_c\n", "LC_4fe4s_c\n", "LC_5drib_c\n", "LC_5mthf_c\n", "LC_6pgc_c\n", "LC_6pgl_c\n", "LC_ac_c\n", "LC_acald_c\n", "LC_accoa_c\n", "LC_acon-C_c\n", "LC_actp_c\n", "LC_adp_c\n", "LC_adpglc_c\n", "LC_akg_c\n", "LC_ala-L_c\n", "LC_amet_c\n", "LC_amob_c\n", "LC_amp_c\n", "LC_arg-L_c\n", "LC_asn-L_c\n", "LC_asp-L_c\n", "LC_atp_c\n", "LC_bglycogen_c\n", "LC_bmocogdp_c\n", "LC_btn_c\n", "LC_chor_c\n", "LC_cit_c\n", "LC_co2_c\n", "LC_coa_c\n", "LC_colipa_c\n", "LC_ctp_c\n", "LC_cys-L_c\n", "LC_datp_c\n", "LC_db4p_c\n", "LC_dctp_c\n", "LC_dgtp_c\n", "LC_dha_c\n", "LC_dhap_c\n", "LC_dhor-S_c\n", "LC_dttp_c\n", "LC_e4p_c\n", "LC_etoh_c\n", "LC_f6p_c\n", "LC_fad_c\n", "LC_fdp_c\n", "LC_fe2_c\n", "LC_fe3_c\n", "LC_for_c\n", "LC_fum_c\n", "LC_g1p_c\n", "LC_g3p_c\n", "LC_g6p_c\n", "LC_glc-D_c\n", "LC_gln-L_c\n", "LC_glu-L_c\n", "LC_glx_c\n", "LC_gly_c\n", "LC_glyc_c\n", "LC_glyc-R_c\n", "LC_glyc3p_c\n", "LC_glyclt_c\n", "LC_glycogen_c\n", "LC_glycogenn1_c\n", "LC_gthrd_c\n", "LC_gtp_c\n", "LC_h_c\n", "LC_h2o_c\n", "LC_h2o2_c\n", "LC_his-L_c\n", "LC_iasp_c\n", "LC_icit_c\n", "LC_ile-L_c\n", "LC_lac-D_c\n", "LC_leu-L_c\n", "LC_lipopb_c\n", "LC_lys-L_c\n", "LC_mal-L_c\n", "LC_malcoa_c\n", "LC_met-L_c\n", "LC_mlthf_c\n", "LC_mql8_c\n", "LC_mqn8_c\n", "LC_mththf_c\n", "LC_nad_c\n", "LC_nadh_c\n", "LC_nadp_c\n", "LC_nadph_c\n", "LC_nh4_c\n", "LC_o2_c\n", "LC_oaa_c\n", "LC_orot_c\n", "LC_pe160_c\n", "LC_pe161_c\n", "LC_pe181_c\n", "LC_pep_c\n", "LC_pg160_c\n", "LC_pg161_c\n", "LC_pg181_c\n", "LC_phe-L_c\n", "LC_pi_c\n", "LC_ppi_c\n", "LC_pro-L_c\n", "LC_ptrc_c\n", "LC_pydx5p_c\n", "LC_pyr_c\n", "LC_q8_c\n", "LC_q8h2_c\n", "LC_r5p_c\n", "LC_ribflv_c\n", "LC_ru5p-D_c\n", "LC_s17bp_c\n", "LC_s7p_c\n", "LC_ser-L_c\n", "LC_so4_c\n", "LC_spmd_c\n", "LC_succ_c\n", "LC_succoa_c\n", "LC_thf_c\n", "LC_thmpp_c\n", "LC_thr-L_c\n", "LC_trdox_c\n", "LC_trdrd_c\n", "LC_trp-L_c\n", "LC_tyr-L_c\n", "LC_udcpdp_c\n", "LC_utp_c\n", "LC_val-L_c\n", "LC_xu5p-D_c\n", "LC_5mtr_e\n", "LC_ac_e\n", "LC_acald_e\n", "LC_acser_e\n", "LC_akg_e\n", "LC_ala-D_e\n", "LC_ala-L_e\n", "LC_alaala_e\n", "LC_arg-L_e\n", "LC_asn-L_e\n", "LC_asp-L_e\n", "LC_cit_e\n", "LC_co2_e\n", "LC_colipa_e\n", "LC_cys-L_e\n", "LC_dha_e\n", "LC_etoh_e\n", "LC_fe2_e\n", "LC_fe3_e\n", "LC_for_e\n", "LC_glc-D_e\n", "LC_glcn_e\n", "LC_glu-L_e\n", "LC_gly_e\n", "LC_glyc_e\n", "LC_glyc-R_e\n", "LC_glyclt_e\n", "LC_gthrd_e\n", "LC_h_e\n", "LC_h2o_e\n", "LC_his-L_e\n", "LC_ile-L_e\n", "LC_lac-D_e\n", "LC_leu-L_e\n", "LC_lys-L_e\n", "LC_mal-L_e\n", "LC_meoh_e\n", "LC_nh4_e\n", "LC_o2_e\n", "LC_phe-L_e\n", "LC_pi_e\n", "LC_pro-L_e\n", "LC_pyr_e\n", "LC_ser-L_e\n", "LC_so4_e\n", "LC_spmd_e\n", "LC_succ_e\n", "LC_thr-L_e\n", "LC_trp-L_e\n", "LC_tyr-L_e\n", "LC_val-L_e\n", "LC_5mtr_p\n", "LC_ac_p\n", "LC_acald_p\n", "LC_acser_p\n", "LC_akg_p\n", "LC_ala-D_p\n", "LC_ala-L_p\n", "LC_alaala_p\n", "LC_arg-L_p\n", "LC_asn-L_p\n", "LC_asp-L_p\n", "LC_cit_p\n", "LC_clpn160_p\n", "LC_clpn161_p\n", "LC_clpn181_p\n", "LC_co2_p\n", "LC_colipa_p\n", "LC_cys-L_p\n", "LC_dha_p\n", "LC_etoh_p\n", "LC_fe2_p\n", "LC_fe3_p\n", "LC_for_p\n", "LC_glc-D_p\n", "LC_glcn_p\n", "LC_glu-L_p\n", "LC_gly_p\n", "LC_glyc_p\n", "LC_glyc-R_p\n", "LC_glyclt_p\n", "LC_gthrd_p\n", "LC_h_p\n", "LC_h2o_p\n", "LC_his-L_p\n", "LC_ile-L_p\n", "LC_lac-D_p\n", "LC_leu-L_p\n", "LC_lys-L_p\n", "LC_mal-L_p\n", "LC_meoh_p\n", "LC_murein3p3p_p\n", "LC_murein3px4p_p\n", "LC_murein4p4p_p\n", "LC_murein4px4p_p\n", "LC_murein4px4px4p_p\n", "LC_nh4_p\n", "LC_o2_p\n", "LC_pe160_p\n", "LC_pe161_p\n", "LC_pe181_p\n", "LC_pg160_p\n", "LC_pg161_p\n", "LC_pg181_p\n", "LC_phe-L_p\n", "LC_pi_p\n", "LC_pro-L_p\n", "LC_pyr_p\n", "LC_ser-L_p\n", "LC_so4_p\n", "LC_spmd_p\n", "LC_succ_p\n", "LC_thr-L_p\n", "LC_trp-L_p\n", "LC_tyr-L_p\n", "LC_val-L_p\n", "LC_ade_p\n" ] } ], "source": [ "#let's integrate metabolomcis data into tmodel, here you can see the Log Concentration avriables. Each one of them has a lower and upper bound we can modify based on out omics data.\n", "for var in tmodel.variables:\n", " if(\"LC_\" in var.name):\n", " print(var.name)" ] }, { "cell_type": "code", "execution_count": 52, "id": "a4b932f8", "metadata": {}, "outputs": [], "source": [ "for met_id in mets_model_metabolomics:\n", "\n", " val = metabolomics.loc[met_id, 'Concentration in mmol'] \n", "\n", " # we can derive a LB and UB by using the deviation column \n", "\n", " err = metabolomics.loc[met_id, 'Standard deviation']\n", "\n", " lb = metabolomics.loc[met_id, 'Concentration in mmol'] - err\n", " ub = metabolomics.loc[met_id, 'Concentration in mmol'] + err\n", "\n", " variable = tmodel.variables.get(\"LC_\" + met_id)\n", "\n", " if not pd.isna(err):\n", " if (val - err) > 0:\n", " variable.ub = np.log(1e-3 * (val + err))\n", " variable.lb = np.log(1e-3 * (val - err))\n", " elif (val - err) <= 0 and variable.lb <= np.log(1e-3 * (val + err)): # convert mmol to mol\n", " variable.ub = np.log(1e-3 * (val + err))\n", " else:\n", " variable.ub = np.log(1e-3 * (val * 1.5))\n", " variable.lb = np.log(1e-3 * (val * 0.5))\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 53, "id": "8d7fb4b4", "metadata": {}, "outputs": [ { "data": { "text/html": [ "Optimal solution with objective value 0.811
\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", " \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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
fluxesreduced_costs
DM_4CRSOL0.000181None
DM_5DRIB0.000187None
DM_AMOB0.000002None
DM_MTHTHF0.001087None
Ec_biomass_iJO1366_WT_53p95M0.810997None
.........
LMPD_250_trp-L_c0.044795None
LMPD_251_tyr-L_c0.108668None
LMPD_252_udcpdp_c0.000045None
LMPD_253_utp_c0.113622None
LMPD_254_val-L_c0.333469None
\n", "

599 rows × 2 columns

\n", "
" ], "text/plain": [ "" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#check feasibility\n", "tmodel.optimize()" ] }, { "cell_type": "code", "execution_count": 54, "id": "0276e71c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Working with model: FBA\n", "\n", "Condition: aerobic\n", "\n", "Condition: anaerobic\n", "\n", "Working with model: Thermo\n", "\n", "Condition: aerobic\n", "\n", "Condition: anaerobic\n" ] } ], "source": [ "# First we work with the FBA model\n", "for model_selected in [\"FBA\", \"Thermo\"]:\n", "\n", " if(model_selected == \"Thermo\"):\n", " selected_model = tmodel\n", " else:\n", " selected_model = model\n", " print(\"\\nWorking with model: \" + model_selected)\n", "\n", " for condition in [\"aerobic\", \"anaerobic\"]:\n", "\n", " print(\"\\nCondition: \" + condition)\n", "\n", " if condition == \"aerobic\":\n", " selected_model.reactions.get_by_id(\"DM_o2_e\").lower_bound = -20 # We allow unlimited oxygen uptake\n", " else:\n", " selected_model.reactions.get_by_id(\"DM_o2_e\").lower_bound = 0 # We do not allow oxygen uptake\n", "\n", "\n", " for key, val in dict_carbon_sources.items():\n", "\n", " #allow uptake of one carbon source\n", " selected_model.reactions.get_by_id(key).lower_bound = -10\n", "\n", " #fix all others to 0\n", " for key2, val2 in dict_carbon_sources.items():\n", " if(key2 != key):\n", " selected_model.reactions.get_by_id(key2).lower_bound = 0\n", " \n", " #optimize\n", " sol = selected_model.optimize()\n", " if(sol.status != 'optimal'):\n", " print(\"No optimal solution found for carbon source \" + val + \" (\" + key + \")\")\n", " sol_f = math.nan\n", " else:\n", " sol_f = sol.objective_value\n", "\n", " df_res2.loc[key, model_selected + \" growth rate (\" + condition + \")\"] = sol_f\n", "\n", " selected_model = None\n", "\n", "#restore original bounds\n", "for rxn in model.reactions:\n", " rxn.lower_bound = original_bounds[rxn.id][0]\n", " rxn.upper_bound = original_bounds[rxn.id][1]\n", "for rxn in tmodel.reactions:\n", " rxn.lower_bound = original_bounds[rxn.id][0]\n", " rxn.upper_bound = original_bounds[rxn.id][1]\n", " " ] }, { "cell_type": "code", "execution_count": 55, "id": "d6284cad", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
FBA growth rate (aerobic)FBA growth rate (anaerobic)Thermo growth rate (aerobic)Thermo growth rate (anaerobic)
DM_glc_e0.9938260.2731290.9938690.18888
DM_lac-D_e0.4392140.0247510.00.0
DM_ac_e0.2594010.00.00.0
DM_etoh_e0.4323350.00.00.0
\n", "
" ], "text/plain": [ " FBA growth rate (aerobic) FBA growth rate (anaerobic) \\\n", "DM_glc_e 0.993826 0.273129 \n", "DM_lac-D_e 0.439214 0.024751 \n", "DM_ac_e 0.259401 0.0 \n", "DM_etoh_e 0.432335 0.0 \n", "\n", " Thermo growth rate (aerobic) Thermo growth rate (anaerobic) \n", "DM_glc_e 0.993869 0.18888 \n", "DM_lac-D_e 0.0 0.0 \n", "DM_ac_e 0.0 0.0 \n", "DM_etoh_e 0.0 0.0 " ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_res2" ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "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.18" } }, "nbformat": 4, "nbformat_minor": 5 }