Course GEM and FBA _3 Manipulating models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Manipulating models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Genome-scale metabolic models can be modified to simulate, for example, the effects on fluxes when a gene is knocked out or mutated.This is extremely practical and time-saving for effective strain development and increased understanding of an organism.\n",
    "\n",
    "First, the ecoli_core_model is loaded into the jupyter notebook again using COBRApy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: cobra in c:\\programdata\\anaconda3\\lib\\site-packages (0.18.1)\n",
      "Requirement already satisfied: depinfo in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (1.5.3)\n",
      "Requirement already satisfied: pandas>=0.17.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (1.0.1)\n",
      "Requirement already satisfied: future in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (0.18.2)\n",
      "Requirement already satisfied: six in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (1.14.0)\n",
      "Requirement already satisfied: ruamel.yaml>=0.16 in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (0.16.10)\n",
      "Requirement already satisfied: python-libsbml-experimental==5.18.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (5.18.1)\n",
      "Requirement already satisfied: swiglpk in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (4.65.1)\n",
      "Requirement already satisfied: optlang>=1.4.2 in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (1.4.4)\n",
      "Requirement already satisfied: numpy>=1.13 in c:\\programdata\\anaconda3\\lib\\site-packages (from cobra) (1.18.1)\n",
      "Requirement already satisfied: pipdeptree>=0.13.2 in c:\\programdata\\anaconda3\\lib\\site-packages (from depinfo->cobra) (1.0.0)\n",
      "Requirement already satisfied: pytz>=2017.2 in c:\\programdata\\anaconda3\\lib\\site-packages (from pandas>=0.17.0->cobra) (2019.3)\n",
      "Requirement already satisfied: python-dateutil>=2.6.1 in c:\\programdata\\anaconda3\\lib\\site-packages (from pandas>=0.17.0->cobra) (2.8.1)\n",
      "Requirement already satisfied: ruamel.yaml.clib>=0.1.2; platform_python_implementation == \"CPython\" and python_version < \"3.9\" in c:\\programdata\\anaconda3\\lib\\site-packages (from ruamel.yaml>=0.16->cobra) (0.2.0)\n",
      "Requirement already satisfied: sympy>=1.0.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from optlang>=1.4.2->cobra) (1.5.1)\n",
      "Requirement already satisfied: pip>=6.0.0 in c:\\programdata\\anaconda3\\lib\\site-packages (from pipdeptree>=0.13.2->depinfo->cobra) (20.0.2)\n",
      "Requirement already satisfied: mpmath>=0.19 in c:\\programdata\\anaconda3\\lib\\site-packages (from sympy>=1.0.0->optlang>=1.4.2->cobra) (1.1.0)\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "pip install cobra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra.io import read_sbml_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "model=read_sbml_model('e_coli_core.xml')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making temporary changes to the model\n",
    "\n",
    "Usually one relies on making copies if objects need to be changed but the original state needs to be retained. Unfortunately, making copies of models is time consuming.\n",
    "\n",
    "To repeat: In this case the stoichiometric processes of the e.coli_core_model are considered, which are responsible for biomass production and this production should be increased.\n",
    "\n",
    "The following command displays all reactions and flows that are optimised to achieve this goal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong><em>Optimal</em> solution with objective value 0.874</strong><br><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fluxes</th>\n",
       "      <th>reduced_costs</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>PFK</th>\n",
       "      <td>7.477382</td>\n",
       "      <td>-1.387779e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PFL</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.081668e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PGI</th>\n",
       "      <td>4.860861</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PGK</th>\n",
       "      <td>-16.023526</td>\n",
       "      <td>6.938894e-18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PGL</th>\n",
       "      <td>4.959985</td>\n",
       "      <td>1.517883e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NADH16</th>\n",
       "      <td>38.534610</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NADTRHD</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>-2.546243e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NH4t</th>\n",
       "      <td>4.765319</td>\n",
       "      <td>-1.387779e-17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>O2t</th>\n",
       "      <td>21.799493</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PDH</th>\n",
       "      <td>9.282533</td>\n",
       "      <td>-1.387779e-17</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>95 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "<Solution 0.874 at 0x1271a39a188>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 192 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "copy_of_model=model.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes, even milliseconds add up pretty quickly if you need to run many simulation (e.g. if you need to knock out every single gene individually in the model to check if it is essential or not)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 12.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for gene in model.genes:\n",
    "    mutant=model.copy()\n",
    "    mutant.genes.get_by_id(gene.id).knock_out()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For that reason cobrapy provides a mechanism that is less time consuming. Almost all methods that make changes to the mdoel such as knocking-out genes, reactions, adding or removing metabolites, reactions etc can be automatically reverted upon exit from a python context. How this works is probably best understood by looking at an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 24 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with model:\n",
    "    for gene in model.genes:\n",
    "        gene.knock_out()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the 'with model' statements starts the context and changes done to the model one indentation level to the right, are automatically recorded. When that block finishes, the context manager is requested to roll-back all changes leaving the model looking exactly as it did before all the changes.\n",
    "\n",
    "Changing flux bounds can as indicated also be done reversibly. For example let's set the lower and upper bound of phosphoglycerate kinase (PGK) to 0 (effectively knocking out the reaction)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGK's bounds inside the with statement\n",
      "0 (0, 0)\n",
      "Mutant growth rate:  -5.487661025838073e-16\n",
      "PGK's bounds outside the with statement\n",
      "(-1000.0, 1000.0)\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    model.reactions.PGK.bounds = 0, 0\n",
    "    print(\"PGK's bounds inside the with statement\")\n",
    "    print(model.reactions.PGK.lower_bound, model.reactions.PGK.bounds)\n",
    "    print('Mutant growth rate: ', model.optimize().objective_value)\n",
    "print(\"PGK's bounds outside the with statement\")\n",
    "print(model.reactions.PGK.bounds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slim versus full optimize\n",
    "\n",
    "Mathematical solvers are now so fast that for many small to mid-size models computing the solution can be even faster than it takes us to collect the values from the solver and convert that to objects that are usable for in python. When we use model.optimize we gather values for all reactions and metabolites and that can take some time. If we are only interested in the flux value of a single reaction or the objective, it is faster to instead use model.slim_optimize which only does the optimization and returns the objective value leaving it up to you to fetch other values that you may need. For example, let's optimize and get the flux value of the ATPM reaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 8 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "8.39"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "solution=model.optimize()\n",
    "solution.fluxes['ATPM']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 0 ns\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "8.39"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "model.slim_optimize()\n",
    "model.reactions.ATPM.flux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing the medium\n",
    "\n",
    "The media composition can also be displayed and, if necessary, optimised in relation to the specified objective of the flux balance analysis. One can access the medium condition using model.medium. The indicated bound is the effective upper uptake bound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'EX_co2_e': 1000.0,\n",
       " 'EX_glc__D_e': 10.0,\n",
       " 'EX_h_e': 1000.0,\n",
       " 'EX_h2o_e': 1000.0,\n",
       " 'EX_nh4_e': 1000.0,\n",
       " 'EX_o2_e': 1000.0,\n",
       " 'EX_pi_e': 1000.0}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.medium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8739215069684324\n"
     ]
    }
   ],
   "source": [
    "solution=model.optimize()\n",
    "print(solution.fluxes['BIOMASS_Ecoli_core_w_GAM'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Changing the carbon source in the medium can be achieved by adjusting the flux bounds of the respective exchange reactions appropriately. For example, the following code block removes glucose from the medium and adds succinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3975630154277627\n"
     ]
    }
   ],
   "source": [
    "medium=model.medium\n",
    "with model:\n",
    "    medium['EX_glc__D_e']=0\n",
    "    medium['EX_succ_e']=10\n",
    "    model.medium=medium\n",
    "    solution=model.optimize()\n",
    "    print(solution.fluxes['BIOMASS_Ecoli_core_w_GAM'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Changing the carbon source to succinate led to a significant drop in growth rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise (10 min)\n",
    "Change the carbon source in the medium to a different carbon source. What is the difference in the growth rate observed? How about growing E. coli under anaerobic conditions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "medium=model.medium\n",
    "with model:\n",
    "    medium['EX_o2_e']=0\n",
    "    model.medium=medium\n",
    "    solution=model.optimize()\n",
    "    print(solution.fluxes['BIOMASS_Ecoli_core_w_GAM'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding reactions and pathways"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To add reactions and pathways the tools \"Reaction\" and \"Metabolite\" from COBRApy are used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra import Reaction, Metabolite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, let's create a new reactions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_reaction=Reaction('alchemy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This reaction is going to convert water into gold (unfortunately lead is not part of E. coli metabolism; creating wine would be blasphemy). So we need to create a new metabolite, since gold is not yet part of E. coli's native metabolism."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "gold=Metabolite(id='gold_c', compartment='c')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we're going to specify the reaction's stoichiometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_reaction.add_metabolites({model.metabolites.h2o_c: -1, gold: 1})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Printing the reaction reveals that the reaction indeed converts water into gold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h2o_c --> gold_c\n"
     ]
    }
   ],
   "source": [
    "print(new_reaction.build_reaction_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's add the new reaction to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.add_reactions([new_reaction])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quickly check that the reaction was indeed added to the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <table>\n",
       "            <tr>\n",
       "                <td><strong>Reaction identifier</strong></td><td>alchemy</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Name</strong></td><td></td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Memory address</strong></td>\n",
       "                <td>0x01271e186588</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Stoichiometry</strong></td>\n",
       "                <td>\n",
       "                    <p style='text-align:right'>h2o_c --> gold_c</p>\n",
       "                    <p style='text-align:right'>H2O H2O --> </p>\n",
       "                </td>\n",
       "            </tr><tr>\n",
       "                <td><strong>GPR</strong></td><td></td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Lower bound</strong></td><td>0.0</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Upper bound</strong></td><td>1000.0</td>\n",
       "            </tr>\n",
       "        </table>\n",
       "        "
      ],
      "text/plain": [
       "<Reaction alchemy at 0x1271e186588>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.reactions.alchemy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's produce some gold then!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.objective=model.reactions.alchemy\n",
    "model.optimize().objective_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happened? Forgot to add an exchange reaction so that gold can leave the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <table>\n",
       "            <tr>\n",
       "                <td><strong>Reaction identifier</strong></td><td>DM_gold_c</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Name</strong></td><td> demand</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Memory address</strong></td>\n",
       "                <td>0x01271e2f25c8</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Stoichiometry</strong></td>\n",
       "                <td>\n",
       "                    <p style='text-align:right'>gold_c --> </p>\n",
       "                    <p style='text-align:right'> --> </p>\n",
       "                </td>\n",
       "            </tr><tr>\n",
       "                <td><strong>GPR</strong></td><td></td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Lower bound</strong></td><td>0</td>\n",
       "            </tr><tr>\n",
       "                <td><strong>Upper bound</strong></td><td>1000.0</td>\n",
       "            </tr>\n",
       "        </table>\n",
       "        "
      ],
      "text/plain": [
       "<Reaction DM_gold_c at 0x1271e2f25c8>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.add_boundary(model.metabolites.gold_c, type='demand')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1000.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.objective=model.reactions.alchemy\n",
    "model.optimize().objective_value"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}