CBM_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cobra\n",
    "import cobra.core\n",
    "from cobra.core import Model, Reaction, Metabolite\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [Full cobrapy documentation](https://cobrapy.readthedocs.io/en/latest/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: Warming up\n",
    "\n",
    "## Objective: \n",
    "\n",
    "To get familiar with COBRA library by creating a toy model (figure 1) and manipulating it.\n",
    "\n",
    "![title](img/toy_model.png)\n",
    "\n",
    "Figure 1. Toy model with three metabolites (A, B y C), four reactions\n",
    "(v1-v4) and three exchange fluxes (b1-b3). a) Model chart; b) Stoichiometric matrix\n",
    "\n",
    "Documentation: [https://cobrapy.readthedocs.io/en/latest/building_model.html](https://cobrapy.readthedocs.io)/en/latest/building_model.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "toy_model = Model('Toymodel')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = Metabolite(\"A\")\n",
    "A.compartment = 'cytosol'\n",
    "B = Metabolite(\"B\")\n",
    "B.compartment = 'cytosol'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excercise 1.1:\n",
    "\n",
    "Create Metabolite C\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######################\n",
    "# Create Metabolite C\n",
    "######################\n",
    "\n",
    "## TODO\n",
    "## Write your code below\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Add the metabolites to the model\n",
    "toy_model.add_metabolites([A, B, C])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the reactions of a given metabolite\n",
    "print(toy_model.metabolites.A.reactions)\n",
    "# We get an empty set because we haven't created any reaction yet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating the reactions\n",
    "\n",
    "# Create reaction with id b1\n",
    "b1 = Reaction(\"b1\")\n",
    "\n",
    "# To add metabolites to the reactions we need to pass\n",
    "# a dictionary with metabolites as keys\n",
    "# and the stoichiometric coefficients as values\n",
    "b1.add_metabolites({A: 1})\n",
    "\n",
    "# The same is done for the other reactions\n",
    "b2 = Reaction(\"b2\")\n",
    "b2.add_metabolites({B: -1})\n",
    "\n",
    "b3 = Reaction(\"b3\")\n",
    "b3.add_metabolites({C: -1})\n",
    "\n",
    "v1 = Reaction(\"v1\")\n",
    "v1.add_metabolites({A:-1, B:1})\n",
    "\n",
    "v2 = Reaction(\"v2\")\n",
    "v2.add_metabolites({A:-1, C:1})\n",
    "\n",
    "v3 = Reaction(\"v3\")\n",
    "v3.add_metabolites({A:1, C:-1})\n",
    "\n",
    "# Create v4 (Excercise 1.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excercice 1.2: \n",
    "\n",
    "Create reactions v4 (and add its metabolites)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "########################\n",
    "# Create reactions v4\n",
    "########################\n",
    "\n",
    "## TODO\n",
    "## Write your code below\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adding the reactions to the toy model\n",
    "toy_model.add_reactions([b1,b2,b3,v1,v2,v3,v4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write the model in SBML format\n",
    "\n",
    "1. First import the corresponding function\n",
    "2. Write the model\n",
    "3. Optional: You can inspect the SBML using some plain text editor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cobra.io\n",
    "from cobra.io import write_sbml_model\n",
    "\n",
    "# Saving the model to a file in SBML format\n",
    "write_sbml_model(toy_model, \"out/toymodel.sbml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excercice 1.3: Optimisation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting the limits on the inputs\n",
    "toy_model.reactions.b1.upper_bound = 1\n",
    "\n",
    "# Setting a reaction to be optimized as the objective or target\n",
    "toy_model.reactions.b2.objective_coefficient = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To compute a FBA on the model we use the following function:\n",
    "solution = toy_model.optimize()\n",
    "\n",
    "# The result is a solution object which containsthe the following attributes:\n",
    "# objective_value: the objective value of the optimized function (biomass!)\n",
    "print(\"Objective value: %.2f\\n\" % solution.objective_value)\n",
    "# solution.status: shows the status of the solution, it should be optimal\n",
    "# if it is unfeasible, this means that there is no feasible solution\n",
    "print(\"Status: %s\\n\" % solution.status)\n",
    "# solution.fluxes: is a datagrame (pandas) storing the reactions (index) and\n",
    "# their flux values found in the optimal solution.\n",
    "print(\"Fluxes:\\n\")\n",
    "print(solution.fluxes)\n",
    "\n",
    "# Saving the solution into tab-separed-value (tsv) format (plain text)\n",
    "solution.fluxes.to_csv(\"out/toymodel_fba.tsv\", sep=\"\\t\")\n",
    "# Inspect the file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Genome-scale modelling\n",
    "\n",
    "In this part we are gonna use a genome-scale metabolic model of Escherichia coli named iJO1366\n",
    "The file has already been stored in the data folder and its path is data/iJO1366.xml\n",
    "\n",
    "Alternatively, you can also access it here:\n",
    "- [http://bigg.ucsd.edu/models/iJO1366](http://bigg.ucsd.edu/models/iJO1366)\n",
    "\n",
    "to download the model and to see other metadata (citation, description, etc)\n",
    "\n",
    "## Part 2.1: Studying the model\n",
    "\n",
    "### Read the SBML model\n",
    "\n",
    "First we need to import the function read_sbml_model from the cobra.io modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra.io import read_sbml_model\n",
    "\n",
    "# State the path to the file iJO1366.xml\n",
    "sbml_fname = './data/iJO1366.xml'\n",
    "\n",
    "# Read the model\n",
    "model = read_sbml_model(sbml_fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Exercise 2.1: Inspecting the model's numbers\n",
    "\n",
    "First print model description\n",
    "\n",
    "1. print(model)\n",
    "2. Print the total number of reactions: print len(model.reactions)\n",
    "3. Print the total number of metabolties: print len(model.metabolites)\n",
    "4. Print the total number of genes: print len(model.genes)\n",
    "5. Access a particular reaction:\n",
    "   * You can do it directly with: rxn = model.reactions.ENO\n",
    "   * Or you can do use the function get_by_id: rxn = model.reactions.get_by_id('ENO')\n",
    "6. Inspect the reaction by printing:\n",
    "   1. rxn.name\n",
    "   2. rxn.id\n",
    "   3. rxn.bounds\n",
    "   4. rxn.reaction\n",
    "   5. rxn.gene_reaction_rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO\n",
    "## Write your code below\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Exercise 2.2: Inspecting the genes\n",
    "\n",
    "First print the model description\n",
    "\n",
    "1. Access a particular reaction:\n",
    "   * You can do it directly with: gene = model.genes.b0720\n",
    "   * Or you can use the function get_by_id: gene = model.genes.get_by_id('b0720')\n",
    "6. Inspect the reaction by printing:\n",
    "   1. gene.name\n",
    "   2. gene.id\n",
    "   3. gene.reactions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO\n",
    "## Write your code below\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.3: Inspecting the system's boundaries\n",
    "\n",
    "* see the exchange fluxes\n",
    "* see the objective function (the reaction set to be optimized)\n",
    "\n",
    "Use print(model.summary())\n",
    "\n",
    "You can also find the objective function using the following filtering technique:\n",
    "* [r for r in model.reactions if r.objective_coefficient == 1]\n",
    "* the reaction id of the biomass is Ec_biomass_iJO1366_WT_53p95M\n",
    "and the exchange fluxes can be accessed using:\n",
    "* model.boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO\n",
    "## Write your code below\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2.2 Running a Flux Balance Analysis (FBA)\n",
    "\n",
    "Documentation: [https://cobrapy.readthedocs.io/en/latest/simulating.html](https://cobrapy.readthedocs.io/en/latest)/simulating.html\n",
    "\n",
    "By default, the model boundary condition (growth medium) is M9 aerobic (glucose minimal)\n",
    "\n",
    "1.  Check the medium by inspecting the lower_bound of the following reactions:\n",
    "  * EX\\_glc\\_e\\_.lower_bound\n",
    "  * EX\\_o2\\_e\\_.lower_bound\n",
    "2.  Optimize biomass using: \n",
    "  * solution = model.optimize()\n",
    "  \n",
    "3.  Inspect the solution as we did previously in Part 1.2 Optimization.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = model.optimize()\n",
    "\n",
    "print(\"Objective value: %.2f\\n\" % solution.objective_value)\n",
    "print(\"Status: %s\\n\" % solution.status)\n",
    "\n",
    "print(\"Fluxes:\\n\")\n",
    "print(solution.fluxes)\n",
    "\n",
    "# Saving the solution into tab-separed-value (tsv) format (plain text)\n",
    "solution.fluxes.to_csv(\"out/iJO1366_fba.tsv\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2.4: Identify reactions in the iJO1366 simulation\n",
    "\n",
    "Inspect the flux value of the following reactions\n",
    "* The glucose consumption: EX_glc_e_\n",
    "* The oxygen consumption: EX_o2_e_\n",
    "* The biomass reaction: Ec_biomass_iJO1366_WT_53p95M\n",
    "\n",
    "HINT: use the solution object -> solution.fluxes.reaction_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO\n",
    "## Write your code below\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Knock out studies.\n",
    "\n",
    "### Exercise 3.1: Single knock out study\n",
    "\n",
    "Documentation: [https://cobrapy.readthedocs.io/en/latest/deletions.html#Knocking-out-single-genes-and-reactions](https://cobrapy.readthedocs.io/en/latest/deletions.html#Knocking-out-single-genes-and-reactions)\n",
    "\n",
    "We will use gene b0720 as an example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cobra.manipulation import find_gene_knockout_reactions\n",
    "\n",
    "# Pick a gene of interest\n",
    "gene =  model.genes.b0720\n",
    "\n",
    "# Inspect the reactions associated to b0720\n",
    "print(\"id\\treaction_name\")\n",
    "for r in gene.reactions: \n",
    "    print(\"%s \\t%s\" % (r.id,r.name))\n",
    "\n",
    "print()\n",
    "# We can also check the genes associated to this reaction\n",
    "reaction = model.reactions.CS\n",
    "print(\"GPR:\",reaction.gene_reaction_rule)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, COBRA can find the proper reaction to be disabled when a gene is knocked out as follows:\n",
    "\n",
    "gene =  model.genes.b0720\n",
    " \n",
    "     with model:\n",
    "         gene.knock_out()\n",
    "         ko_solution = model.optimize()\n",
    "\n",
    "This codes knocks out the gene b0720, recalculates the FBA and stores the new solution in ko_solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If we perform the knockout in the \"with\" block we don't need to care\n",
    "# about restoring the knocked out gene afterwards; it is automatically restored out of the \"with\" block.\n",
    "with model:\n",
    "    gene.knock_out()\n",
    "    ko_solution = model.optimize()\n",
    "\n",
    "########################################################################################    \n",
    "# TODO\n",
    "# Check the growth value (Hint: ko_solution.fluxes.Ec_biomass_iJO1366_WT_53p95M or ko_solution.objective_values)\n",
    "# What happened?\n",
    "\n",
    "## write your code below\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Go to the Ecocyc database and check the in vivo experimental result for the knockout of b0720\n",
    "by accessing the following link:\n",
    "* [https://ecocyc.org/gene?orgid=ECOLI&id=EG10402](https://ecocyc.org/gene?orgid=ECOLI&id=EG10402)\n",
    "\n",
    "Is b0720 essential or not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3.2: Systems-wide knock out study of *E. coli*\n",
    "    \n",
    "COBRA has a special function to run the single gene knock outs of a list of genes. \n",
    "\n",
    "The function's name is single_gene_deletion\n",
    "\n",
    "First import the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the function single_gene_deletion\n",
    "from cobra.flux_analysis import single_gene_deletion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Then get the list of all the genes\n",
    "all_genes = [g.id for g in model.genes]\n",
    "\n",
    "# Running in silico (takes a while)\n",
    "knockout = single_gene_deletion(model, gene_list=all_genes)\n",
    "\n",
    "# This is a fix to get the gene's id as the index\n",
    "index_mapper = {i:list(i)[0] for i in knockout.index}\n",
    "knockout = knockout.rename(mapper=index_mapper, axis=0)\n",
    "\n",
    "# The output of the function single_gene_deletion is a dataframe\n",
    "print(knockout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We define a threshold to define whether the reduction on the biomass flux is considered lethal.\n",
    "threshold = 0.01\n",
    "\n",
    "# Use this threshold to find which set of genes' knock out reduce the predicted growth below the threshold.\n",
    "insilico_lethals = set(knockout.index[knockout.growth< threshold])\n",
    "# The set of non-essential genes are the genes with a growth value above the threshold.\n",
    "insilico_non_lethals = set(knockout.index[knockout.growth > threshold])\n",
    "\n",
    "print(\"in silico lethals:\", len(insilico_lethals))\n",
    "print(\"in silico non lethals:\", len(insilico_non_lethals))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we need to experimentally verify essential and non-essential gene sets.\n",
    "\n",
    "# Read the set of essential genes in vivo\n",
    "import json\n",
    "fname = './data/m9_invivo_lethals.json'\n",
    "with open(fname) as fh:\n",
    "    invivo_lethals = json.load(fh)\n",
    "    invivo_lethals = set(invivo_lethals)\n",
    "\n",
    "    \n",
    "# Convert the list of all model genes into a set.\n",
    "all_genes = set([g.id for g in model.genes])\n",
    "\n",
    "# We can use the difference to obtain the list of in vivo non-lethals\n",
    "invivo_non_lethals = all_genes - invivo_lethals\n",
    "\n",
    "# Print the size of both sets\n",
    "print(\"in vivo lethals:\", len(invivo_lethals))\n",
    "print(\"in vivo non lethals:\", len(invivo_non_lethals))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://en.wikipedia.org/wiki/Receiver_operating_characteristic\n",
    "\n",
    "# True Positives, genes predicted as essentials that are essentials in vivo (correctly predicted)\n",
    "TP =  insilico_lethals & invivo_lethals\n",
    "\n",
    "# True Negatives, genes predicted as NON-essentials that are NON-essential in vivo (correctly predicted)\n",
    "TN = insilico_non_lethals & invivo_non_lethals\n",
    "\n",
    "# False Positives, wrongly predicted as NON-essential genes\n",
    "FN = insilico_non_lethals & invivo_lethals\n",
    "\n",
    "# False Positives, wrongly predicted as essential genes\n",
    "FP = insilico_lethals & invivo_non_lethals\n",
    "\n",
    "# True in vivo esssential genes\n",
    "P = TP | FN\n",
    "# True in vivo NON-esssential genes\n",
    "N = TN | FP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluated Exercises\n",
    "\n",
    "\n",
    "### Excercise 1:\n",
    "Complete the following table\n",
    "\n",
    "| In vivo \\ In silico        | in silico lethal | in silico non-lethal |\n",
    "| -------------------------- |:----------------:| --------------------:|\n",
    "| <b>in vivo lethal</b>      |               ?  |                   ?  |\n",
    "| <b>in vivo non-lehtal</b>  |               ?  |                   ?  |\n",
    "\n",
    "\n",
    "### Excercise 2:\n",
    "Acces the following link:\n",
    "\n",
    "https://en.wikipedia.org/wiki/Sensitivity_and_specificity\n",
    "\n",
    "Get the formulas and calculate the following metrics:\n",
    "* sensitivity\n",
    "* specificity\n",
    "* precision\n",
    "* accuracy\n",
    "\n",
    "### Excercise 3:\n",
    "In one paragraph, comment the predictive capacity of the model and briefly discuss the possible sources of errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sensitivity, recall, hit rate, or true positive rate (TPR)\n",
    "# We computed the sensitivity as follows\n",
    "sensitivity = len(TP) / len(P) \n",
    "\n",
    "# TODO\n",
    "# complete the following code\n",
    "\n",
    "# Specificity, selectivity or true negative rate (TNR)\n",
    "specificity = ## COMPLETE HERE \n",
    "\n",
    "# Precision or positive predictive value (PPV)\n",
    "precision =  ## COMPLETE HERE\n",
    "\n",
    "# Accuracy (ACC)\n",
    "accuracy =  ## COMPLETE HERE\n",
    "\n",
    "# Print the different values and discuss their meanings"
   ]
  }
 ],
 "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.6.7"
  },
  "latex_metadata": {
   "affiliation": "Barcelona Supercomputing Center, Barcelona Spain",
   "author": "Miguel Ponce de Leon",
   "title": "Constraint-based modeling tutorial"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}