QMiBi_GSMM-cobrapy_L2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/uliebal/RWTH-QMB1/blob/master/2001_QuantMiBi1_CobraPy_SimVal.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "dPnMzwQ_IWwZ"
   },
   "source": [
    "# Seminar Quantitativ Microbiology 2: \n",
    "# Simulation of Genome Scale Metabolic Models with CobraPy\n",
    "\n",
    "## Introduction\n",
    "The seminar provides a guide of how to work with genome scale metabolic models (GSMM) of micro-organisms. This seminar extends the introduction to CobraPy from the seminar in Quantitative Microbiology 1. The goal of this seminar is to identify minimal medium composition, extract information about and selecting the appropriate biomass composition formula and testing the reproduction of experimental data by the GSMM. We examine a recent model of *P. pastoris* ([Tomas-Gamisans et al., 2017](https://dx.doi.org/10.1111/1751-7915.12871)). This organism is of biotechnological relevance because it can glycosylate recombinant proteins for use as drugs and it can metabolize methanol as a potential alternative to petrol based chemistry ([Liebal et al., 2018](https://dx.doi.org/10.1016/j.mec.2018.e00075)).\n",
    "\n",
    "We will analyse the GSMM of *P. pastoris* to estimate for the exchange reactions the range of permissible flux values. This indicates which phenotypes we can expect during standard cultivation and can highlight easily overproduced metabolites. Subsequently, we will identify the minimal medium composition. An important property of GSMM is their ability to predict growth rates. We will extract experimentally measured growth rates for various substrates and compare them with predictions of the model. The data reproduction is exemplified with the substrates of methanol and glycerol and the self-learning task is to supplement growth rates based on glucose uptake rate from literature ([Lehnen et al., 2017](https://dx.doi.org/10.1016/j.meteno.2017.07.001)).\n",
    "\n",
    "The seminar uses Jupyter notebooks, a new way to use and visualize code in the web. Such a notebook is composed of a sequence of cells. Cellls can be either text/comments, like this introduction, or it contains python code to be run. In this example the code is evaluated by the cloud service [Binder](https://mybinder.org/). The output for each code-cell is shown directly below it. For a overview on Jupyter notebooks read [this review](https://www.nature.com/articles/d41586-018-07196-1). Another usefull resource to develop Jupyter notebooks is via [Google Colaboration](https://colab.research.google.com).\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "## Tutorial Steps\n",
    "  * Set up of Python environment\n",
    "    * Basic libraries(sys, pandas, numpy, matplotlib, zipfile, cobrapy)\n",
    "  * Analysis of Genome Scale Metabolic Model\n",
    "    * Retrieval of GSMM for *P. pastoris*\n",
    "    * Flux variability of exchange reactions\n",
    "    * Minimal medium composition\n",
    "  * Experimental growth rate reproduction\n",
    "    * Familiarizing with biomass composition reactions\n",
    "    * Defining functions for correct biomass equation switch\n",
    "    * Data retrieval\n",
    "    * Simulation loop\n",
    "    * Graphical output\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XRP3j145RXOi"
   },
   "source": [
    "## Set-up compute environment\n",
    "\n",
    "Before we can analyse GSMM, we have adjust the python environment that it integrates the cobrapy toolbox and downloading the GSMM."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ewTIEAb4Rb2G"
   },
   "source": [
    "### Basic Python libraries \n",
    "Some libraries that facilitate data manipulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6wxIUGOEReDX"
   },
   "outputs": [],
   "source": [
    "import sys # loading commands to control/navigate within the system architecture\n",
    "# Loading pandas, a library for data manipulation\n",
    "import pandas as pd\n",
    "\n",
    "# Loading numpy, a library fo manipulation of numbers\n",
    "import numpy as np\n",
    "\n",
    "# loading matplotlib, a library for visualization\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# in order to extraact the GSMM we need to process zip-files\n",
    "import zipfile\n",
    "\n",
    "# loading cobrapy, a library dedicated to the analysis of genome scale metabolic models\n",
    "from cobra.io import read_sbml_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4JJ24oZJx0Tr"
   },
   "source": [
    "## Analysis of Genome Scale Metabolic Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "lTzR-LCwR18N"
   },
   "source": [
    "### Retrieval of GSMM for *P. pastoris*\n",
    "\n",
    "*P. pastoris* is an important biotechnological organism and several GSMMs have been generated. We will use the model generated by [Tomas-Gamisans et al., 2017](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5743807/#!po=67.8571). Please click on the link and find the additional data file for the model iMT1026v3 in the chapter \"Support Information\", right click on \"Click here for additional data file\" and select \"Copy link address\". In the next cell, replace the '?..?' after \"!wget\" with the link address by pressing Ctrl+V.  Then press Ctrl+Enter again to to download the zipped model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "J29I6wOhR261"
   },
   "outputs": [],
   "source": [
    "!wget ?..?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "9XwIzjaSOyb1"
   },
   "source": [
    "**When using Binder:**\n",
    "Go back to \"Home\"-Tab in your browser. In the list of files you should see the zip-file \"MBT2-11-224-s003.zip\". If this is not the case, click on the button with the two right turning arrows (Refresh notebook files). On the basis of this you can check once again whether the required zip file of the desired model has been saved. Now go back to the simulation-Tab. Replace the ?..? in the window below with the name of the zip file and then press Ctrl+Enter to unzip this file.\n",
    "\n",
    "**When using Colab:**\n",
    "Click on \"Files\" on the left side of this window. You should now see the file MBT2-11-224-s003.zip. If this is not the case, click on \"Refresh\". Right-click on this zip file and copy it by clicking \"Copy path\". To extract the zip file, replace the '?..?' in the window below with the file path by pressing Ctrl+V. Then press Ctrl+Enter again to run the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "er5OlV2xNCeX"
   },
   "outputs": [],
   "source": [
    "with zipfile.ZipFile('?...?', 'r') as zip_ref:\n",
    "  zip_ref.extractall()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "n6sbsuSJS2bb"
   },
   "source": [
    "**When using Binder:**\n",
    "Go back to the \"Home\"-Tab in the browser. In the list of files you should see a new folder with the name \"folder\". If this is not the case, click on the button with the two right turning arrows (Refresh notebook files). This folder contains the xml file of the model. Click on the folder to see the name of the xml-file. Now go back to the simulation-Tab. Replace '?...?' in the cell below by entering the xml-file address as \"folder/filename of xml\". Then press Ctrl+Enter to load the model.\n",
    "\n",
    "**When using Colab:**\n",
    "Click again on \"Files\" on the left side of the window. You should now see a new folder with the name \"folder\". If this is not the case, please click on \"Refresh\". This folder contains the xml file of the model. Click on the folder. Then right click on the xml-file and choose \"Copy path\" again. In the cell below, replace the '?..?' with the copied file path with Ctrl+V which results in \"folder/filename of xml\". Then press Ctrl+Enter to load the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "baip9ogPSXxn"
   },
   "outputs": [],
   "source": [
    "from IPython.utils import io\n",
    "# suppress output, because hundreds of warnings are generated\n",
    "with io.capture_output() as captured:\n",
    "  # generating cobra variable from SBML/xml file\n",
    "  model = read_sbml_model('folder/?...?');\n",
    "print('model loaded!')   \n",
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Jn805mCObF2k"
   },
   "source": [
    "### Flux variability of exchange reactions\n",
    "\n",
    "Flux balance analysis provides a single optimal solution. Mostly, there exist a number of alternative flux distributions around the optimum, which can be physiologically relevant. To identify the variability of exchange fluxes around the optimum solution 'flux variability analysis' can be performed ([Mahadevan & Schilling, 2003](http://dx.doi.org/10.1016/j.ymben.2003.09.002)). Use the following command to identify minimum and maximum flux ranges of the model for the exchange reactions.\n",
    "\n",
    "model.summary(fva=.95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2HiWV97jWhKI"
   },
   "outputs": [],
   "source": [
    "model.summary(fva=.95) # additional argument specifies allowed deviation from the optimum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cvSOmwrZbNIG"
   },
   "source": [
    "### Media test\n",
    "The availability of nutrients has a major impact on metabolic fluxes and cobrapy provides some helpers to manage the exchanges between the external environment and the metabolic model. More detailed descriptions: [here](https://cobrapy.readthedocs.io/en/latest/media.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mEHqAD7nbL8n"
   },
   "outputs": [],
   "source": [
    "from cobra.medium import minimal_medium\n",
    "\n",
    "max_growth = model.slim_optimize()\n",
    "minimal_medium(model, max_growth, minimize_components=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xpFf1wR0QX4L"
   },
   "source": [
    "## Experimental growth rate reproduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JsDpC6gVQigc"
   },
   "source": [
    "### Familiarizing with biomass composition reactions\n",
    "\n",
    "Microorganisms adapt to their substrate. Different substrates provide different energy content and require different cellular resources to become metabolized. In GSMM these differences may be represented by different equations/reactions for the substrates. In iMR1026 for *P. pastoris* there are various biomass equations for glucose, glycerol, glucose-glycerol mixtures, and methanol. When simulating a model, we have to make sure we use the right biomass equation fitting with the substrate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DBNmaT-fRpQz"
   },
   "outputs": [],
   "source": [
    "# List of all reactions with 'BIOMASS' in their name\n",
    "model.reactions.query('?..?')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "82QMmzXYR4br"
   },
   "outputs": [],
   "source": [
    "# Looking in detail to biomass with methanol; use the reaction name given to you from the list of all biomass reactions;\n",
    "model.reactions.?..?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N_eSuR-EyQ-M"
   },
   "source": [
    "### Defining functions for correct biomass equation switch\n",
    "\n",
    "For each substrate, the boundary exchange fluxes are activated and the reactions of competing substrates are disabled.\n",
    "\n",
    "Remember to write a specific function to adapt the model for glucose utilization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XqM2cco0Qpn0"
   },
   "outputs": [],
   "source": [
    "def AdaptMethanol(model, meoh_up):\n",
    "  model.objective = 'BIOMASS_meoh'\n",
    "  # setting uptake reactions right\n",
    "  model.reactions.Ex_glc_D.lower_bound = 0\n",
    "  model.reactions.Ex_glyc.lower_bound = 0\n",
    "  model.reactions.ATPM.lower_bound = .4\n",
    "  model.reactions.Ex_meoh.lower_bound = -np.abs(meoh_up)\n",
    "  # setting additional biomass composition\n",
    "  model.reactions.LIPIDS_meoh.upper_bound = 1000\n",
    "  model.reactions.PROTEINS_meoh.upper_bound = 1000\n",
    "  model.reactions.STEROLS_meoh.upper_bound = 1000\n",
    "  model.reactions.BIOMASS_meoh.upper_bound = 1000\n",
    "  # deactivating Glyc-based biomass composition\n",
    "  model.reactions.LIPIDS_glyc.upper_bound = 0\n",
    "  model.reactions.PROTEINS_glyc.upper_bound = 0\n",
    "  model.reactions.STEROLS_glyc.upper_bound = 0\n",
    "  model.reactions.BIOMASS_glyc.upper_bound = 0  \n",
    "  # deactivating Glc-based biomass composition\n",
    "  model.reactions.LIPIDS.upper_bound = 0\n",
    "  model.reactions.PROTEINS.upper_bound = 0\n",
    "  model.reactions.STEROLS.upper_bound = 0\n",
    "  model.reactions.BIOMASS.upper_bound = 0  \n",
    "  return model\n",
    "\n",
    "def AdaptGlycerol(model, glyc_up):\n",
    "  model.objective = 'BIOMASS_glyc'\n",
    "  # setting uptake reactions right\n",
    "  model.reactions.Ex_meoh.lower_bound = 0\n",
    "  model.reactions.Ex_glc_D.lower_bound = 0\n",
    "  model.reactions.ATPM.lower_bound = 2.5\n",
    "  model.reactions.Ex_glyc.lower_bound = -np.abs(glyc_up)\n",
    "  # setting additional biomass composition\n",
    "  model.reactions.LIPIDS_glyc.upper_bound = 1000\n",
    "  model.reactions.PROTEINS_glyc.upper_bound = 1000\n",
    "  model.reactions.STEROLS_glyc.upper_bound = 1000\n",
    "  model.reactions.BIOMASS_glyc.upper_bound = 1000  \n",
    "  # deactivating MeOH-based biomass composition\n",
    "  model.reactions.LIPIDS_meoh.upper_bound = 0\n",
    "  model.reactions.PROTEINS_meoh.upper_bound = 0\n",
    "  model.reactions.STEROLS_meoh.upper_bound = 0\n",
    "  model.reactions.BIOMASS_meoh.upper_bound = 0\n",
    "  # deactivating Glc-based biomass composition\n",
    "  model.reactions.LIPIDS.upper_bound = 0\n",
    "  model.reactions.PROTEINS.upper_bound = 0\n",
    "  model.reactions.STEROLS.upper_bound = 0\n",
    "  model.reactions.BIOMASS.upper_bound = 0\n",
    "  return model\n",
    "\n",
    "def AdaptGlucose(model, glc_up):\n",
    "  # setting uptake reactions right\n",
    "  model.reactions.Ex_meoh.lower_bound = 0\n",
    "  model.reactions.Ex_glyc.lower_bound = 0\n",
    "  model.reactions.Ex_glc_D.lower_bound = -np.abs(glc_up)\n",
    "  # setting additional biomass composition\n",
    "  model.reactions.LIPIDS.upper_bound = 1000\n",
    "  model.reactions.PROTEINS.upper_bound = 1000\n",
    "  model.reactions.STEROLS.upper_bound = 1000\n",
    "  model.reactions.BIOMASS.upper_bound = 1000  \n",
    "  # deactivating Glyc-based biomass composition\n",
    "  model.reactions.LIPIDS_glyc.upper_bound = 0\n",
    "  model.reactions.PROTEINS_glyc.upper_bound = 0\n",
    "  model.reactions.STEROLS_glyc.upper_bound = 0\n",
    "  model.reactions.BIOMASS_glyc.upper_bound = 0  \n",
    " # deactivating meoh-based biomass composition\n",
    "  model.reactions.LIPIDS_meoh.upper_bound = 0\n",
    "  model.reactions.PROTEINS_meoh.upper_bound = 0\n",
    "  model.reactions.STEROLS_meoh.upper_bound = 0\n",
    "  model.reactions.BIOMASS_meoh.upper_bound = 0\n",
    "  model.objective = 'BIOMASS'\n",
    "  return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "K045st1gSllg"
   },
   "source": [
    "### Data retrieval\n",
    "\n",
    "For evaluation of the growth rate prediction of the *P. pastoris* model we use experimental data from the closely related organism *Ogataea polymorpha*. The measurements in the table are extracted from [van Dijken et al. 1976](https://dx.doi.org/10.1007/bf00446560) for methanol and from [de Koning et al., 1987](https://dx.doi.org/10.1007/BF00456710) and [Moon et al., 2003](https://dx.doi.org/10.1385/ABAB:111:2:65).\n",
    "\n",
    "Data Address: [here](https://rwth-aachen.sciebo.de/s/o72jwWQWh3ame1e/download)\n",
    "\n",
    "**When using Binder:**\n",
    "Click on \"here\" and save the file. Go back to \"Home\"-Tab in the browser afterwards and navigate back to the root directory of your cloud system (i.e. you see the Jupyter notebook files). Click on the \"Upload\" button on the right side of the window and select the file you just saved. Now you should see the file \"Opol-expt-grwth_MeOH-Glyc.csv\", click again on \"Upload\" in the line of the corresponding file. Now go back to the simulation-Tab. Replace '?..?' in the cell below by entering the name of the csv file you downloaded. Then press Ctrl+Enter to load the experimental data.\n",
    "\n",
    "**When using Colab:**\n",
    "Click on \"here\" and save the file. Then click on \"Files\" on the left side of the window. Upload the file you have just saved by clicking on \"Upload\" and selecting this file. Now you should see the file \"Opol-expt-grwth_MeOH-Glyc.csv\" in the \"Files\" section. If this is not the case, click on \"Refresh\". Then right-click on this file and select \"Copy path\". Replace '?..?' in the cell below with the file by pressing Ctrl+V. Then evaluate the cell by Ctrl+Enter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mfOlv0uiCNvI"
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv('?...?', delimiter=',|;', engine='python')\n",
    "data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aMLKO9MNe1fq"
   },
   "source": [
    "Adding measurements for glucose from [Lehnen et al., 2017](https://doi.org/10.1016/j.meteno.2017.07.001) to our basic data table. The required information is the growth rate on glucose for *H. polymorpha*. Extract the necessary information from Table 3 in the article. The reaction name for the 'Exchange'-reaction is 'Ex_glc_D', replace this in the corresponding position of '?..?' in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4mTMJvd0ZHFv"
   },
   "outputs": [],
   "source": [
    "data = data.append({'Substrate':'Glucose', 'Exchange':'?..?', 'uptake rate (mmol/gCDW/h)':?..?, 'growth rate (/h)':?..?, 'source':'Lehnen et al.'}, ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Wxe3Csg0TtYV"
   },
   "source": [
    "### Simulation loop\n",
    "For-Loop over all experimental data points.\n",
    "\n",
    "Remember to add a decision when you include glucose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "TPyYPHrXUMs1"
   },
   "outputs": [],
   "source": [
    "growth_simulated = [];\n",
    "# test_model = model.copy()\n",
    "# iteration over all rows in 'data'\n",
    "for index, row in data.iterrows():\n",
    "  with model as test_model:\n",
    "    print(index) # printing the row number to get feedback that everything is working\n",
    "    # selecting the right substrate in the model based on 'Substrate' in 'data'\n",
    "    if row['Substrate'] == 'Methanol':\n",
    "      model = AdaptMethanol(test_model, row['uptake rate (mmol/gCDW/h)'])\n",
    "    elif row['Substrate'] == 'Glycerol':\n",
    "      model = AdaptGlycerol(test_model, row['uptake rate (mmol/gCDW/h)'])\n",
    "    elif row['Substrate'] == 'Glucose':\n",
    "      test_model = AdaptGlucose(test_model, row['uptake rate (mmol/gCDW/h)'])\n",
    "    else:\n",
    "      print('substrate not considered')      \n",
    "  #     model.optimize()\n",
    "    growth_simulated.append(model.slim_optimize())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aN5SDn5yWJER"
   },
   "source": [
    "### Graphical output\n",
    "\n",
    "Remember to add glucose to the visualization. Add the right axis labels, and a file name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1DG-7h3qWLRi"
   },
   "outputs": [],
   "source": [
    "plt.scatter(data['growth rate (/h)'][0:7], growth_simulated[0:7], s=50, c='k', marker='o');\n",
    "plt.scatter(data['growth rate (/h)'][7], growth_simulated[7], s=50, c='k', marker='s');\n",
    "plt.scatter(data['growth rate (/h)'][8], growth_simulated[8], s=50, c='k', marker='d');\n",
    "plt.xlabel('?..?');\n",
    "plt.ylabel('?..?');\n",
    "myline = np.linspace(0,np.max(growth_simulated),10);\n",
    "plt.plot(myline,myline,'k--');\n",
    "plt.title('Growth rate comparison');\n",
    "plt.legend(['Optimum', 'Methanol (van Dijken)', 'Glycerol (deKoning, Moon)', 'Glucose (Lehnen et al.)'], loc=2);\n",
    "plt.style.use('seaborn-paper')\n",
    "\n",
    "# Saving figure\n",
    "plt.savefig('?..?.png')"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "2001_QuantMiBi1_CobraPy_SimVal.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "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": 1
}