models.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "models.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Gibbons-Lab/isb_course_2021/blob/main/models.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BXWCzV6yAz-D"
      },
      "source": [
        "# 🦠 Metabolic modeling of bacterial isolates 🧫\n",
        "\n",
        "This notebook will accompany the second day of the of the [2021 ISB Virtual Microbiome series](https://isbscience.org/microbiome2021/). The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2021/models). \n",
        "\n",
        "Save your own local copy of this notebook by using `File > Save a copy in Drive`. At some point you may be promted to certify the notebook is safe. We promise that it is 🤞\n",
        "\n",
        "**Disclaimer:**\n",
        "\n",
        "The Google Colab notebook environment will interpret any command as Python code by default. If we want to run bash commands we will have to prefix them by `!`. So any command you see with a leading `!` is a bash command and if you wanted to run it in your terminal you would omit the leading `!`. For example, if in the Colab notebook you ran `!wget` you would just run `wget` in your terminal. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6I42MY9QWpfi"
      },
      "source": [
        "## Setup\n",
        "\n",
        "The setup is a bit simpler this time. We will install a small Python package that was built for the very purpose of this course and installing it will also pull in some additional commands and packages we need, for instance the metabolic modeling package [COBRAPY](https://opencobra.github.io/cobrapy) which we will use extensively."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JS-fUTqAAwlk"
      },
      "source": [
        "!pip install -q git+https://github.com/gibbons-lab/2021_microbiome_course_utils\n",
        "!bootstrap"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G3GbbjGOWgnI"
      },
      "source": [
        "That's it. We are ready to start.\n",
        "\n",
        "# 🧬 From Genome to Model "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kL85Tq0Lnq0F"
      },
      "source": [
        "## The Assembly\n",
        "For our first step we will get our assembled genome. As noted before, these assemblies will not give us a singe contigious piece of DNA, but rather a set of fragmented *contigs*.\n",
        "\n",
        "We will use the course Python package to get our own personal bacterial assembly. **Before running the code below, substitute `[YOUR NAME]` with your full name or nickname.** You can call yourself whatever you'd like, just make sure it's fairly unique. This name will be used to randomly assign you a specific genome assembly from the database. \n",
        "\n",
        "What organism do you get?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EVf8zp5NDha2"
      },
      "source": [
        "from course import gimme_genome\n",
        "\n",
        "assembly = gimme_genome(\"[YOUR NAME]\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-QtV1SRhYYP1"
      },
      "source": [
        "This will download your assembled genome and tell you which organism you are dealing with. Let's take a peak at your assembly. What do you expect it to look like?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xFaxWGXWGBlF"
      },
      "source": [
        "!head *.fna"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yZVa-6SiY4Zd"
      },
      "source": [
        "This output is divided into the ID of the contig and some additional information. `flag` denotes the contig's shape in the contig graph `flag=1` denotes the contig is standalone, `flag=2` denotes a looped path, and `flag=0` is for all other contigs. `multi` is an approximation of the coverage of the contig (how many reads cover each base on average). And `len`, unsurprisingly, is the length of the contig. Your organism will have a whole bunch of these contigs.\n",
        "\n",
        "## The Genes\n",
        "\n",
        "To map our genome to enzymes we will have to identify the genes encoding for these enzymes. For this, the first step is to identify all of the open reading frames (ORFs), or genes, in the genome. Luckily, this is not super hard for bacterial genomes and we can use [prodigal](https://doi.org/10.1186/1471-2105-11-119).\n",
        "\n",
        "What is this weird `${ASSEMBLY}`? Well this is a variable that will be filled in automatically with the ID of your assembly. This way we don't have to copy and paste the name of the assembly everytime."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zZvpMJfBYyFm"
      },
      "source": [
        "!prodigal -p single -i ${ASSEMBLY}.fna \\\n",
        "    -o ${ASSEMBLY}.gff -a ${ASSEMBLY}.faa -q"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M1ypzoh6l0dM"
      },
      "source": [
        "This will predict coding regions in the genome and generate a genome feature file (GFF), along with the amino acid sequences for the predicted proteins. This is what we will use to look for enzymes. Can you guess why we are using the protein sequences and not the nucleotide sequences?\n",
        "\n",
        "Let's have a look at the proteins."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ujgbNx33Z9Ht"
      },
      "source": [
        "!head ${ASSEMBLY}.faa"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xZ_eNnSBmirU"
      },
      "source": [
        "The FASTA headers now indicate where in the contig the gene starts and ends and gives some more info about the coding region.\n",
        "\n",
        "## The Reconstruction\n",
        "\n",
        "Now we have our proteins and can start building the metabolic model. Let's recall the the basics of metabolic reconstruction and what we'll need.\n",
        "\n",
        "![reconstruction methods](https://github.com/Gibbons-Lab/isb_course_2021/raw/main/docs/models/assets/reconstruction.png)\n",
        "\n",
        "Roughly, we need to run the following steps:\n",
        "\n",
        "1. map the identified proteins to our universal model\n",
        "2. find the model most consistent with the genomic evidence that can also produce biomass \n",
        "3. curate the model with various growth media\n",
        "\n",
        "CARVEME can run all of these steps, but we will need to supply the required input data:\n",
        "\n",
        "1. a list of protein sequences found within the genome\n",
        "2. the specific growth medium we want to gapfill the model on\n",
        "\n",
        "But what media should we use for curation? Well, the isolates all come from the human gut, so having a description of the available nutrients there might be helpful. We've chosen to use a medium based on the standard European diet (i.e. a 'Western diet'). Additionally, these isolates were grown in some medium in the lab in order to pick individual colonies for sequencing. Therefore, we will use the anaerobic M9 medium, which is a fairly rich medium commonly used for picking gut isolates. We have pre-prepared these media for you, so let's get them."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vfIBOHRPaIy2"
      },
      "source": [
        "from course import gimme_media\n",
        "\n",
        "media = gimme_media()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JnyyoFKtrhSC"
      },
      "source": [
        "Let's take a look what we got here."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "a50rIKxEHC0q"
      },
      "source": [
        "media"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jGMPySIbsBhk"
      },
      "source": [
        "So we actually have two environmental conditions here. One describes the metabolites expected to be present in the human gut after an average European diet has passed through the small intestine. The other is a well-defined anaerobic M9 medium. We also have some flux bounds that limit the maximum uptake per hour for each of these media.\n",
        "\n",
        "So we got our proteins and the growth media. This is all we need to build the reconstruction...right? \n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "**Houston, we have a problem 🚀**\n",
        "\n",
        "Well it turns out we cannot run CARVEME on Google Colab, because it uses the commercial [CPLEX](https://www.ibm.com/analytics/cplex-optimizer) solver from IBM. IBM does provide free licenses to students and faculty through the [Academic Inititative](https://www.ibm.com/academic/faqs/faqs), but we can't add it to the Colab environment due to licensing restrictions 😭 \n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "But no worries, we have an awesome robot friend (🤖) that can impersonate CARVEME very well. For instance, the following command looks indistinguishable from the actual command you'd normally run in CARVEME:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "W6akG9jOr6w6"
      },
      "source": [
        "!carve -h"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z5diMu84t8BN"
      },
      "source": [
        "So let's use it. We will have to supply the protein FASTA file, the growth media database, and the names of all the growth conditions we want to enable our bacterium to grow on."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dn_-mRjFtsLd"
      },
      "source": [
        "!carve \\\n",
        "    --mediadb carveme_media.tsv --gapfill \"western_diet_gut,M9[-O2]\" \\\n",
        "    --fbc -o ${ASSEMBLY}.xml ${ASSEMBLY}.faa "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mWlYNlMa1fc4"
      },
      "source": [
        "Our CARVEME robot did a great job! Looks just like the real deal. And you know, if it walks like a 🦆 and talks like a 🦆...\n",
        "\n",
        "This generated a file ending in `.xml` which is an SBML file. [SBML](https://sbml.org) is the community standard for sharing metabolic models. And while we are talking about community standards let's also look into model checking a bit."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q9iQxKKNna6z"
      },
      "source": [
        "## The Model Check\n",
        "\n",
        "We talked a bit about model curation and validation and learned that the best way to test a model is using real-world data. But, believe us, there are **plenty** of ways the model can be wrong even before you compare it to data. \n",
        "\n",
        "This has been known in the metabolic modeling community for quite some time, and has been discussed extensively [here](https://www.embopress.org/doi/full/10.15252/msb.20156157) and [here](https://doi.org/10.1038/s41587-020-0446-y), for example. \n",
        "\n",
        "Luckily, based on this prior feedback from the community, there is now a nice tool called [MEMOTE](https://memote.io), which allows you to upload an SBML model and run it through a whole battery of checks and tests. If we were to do this during the course, we would spend most of our time waiting around, so for now, we have supplied you with the pre-computed MEMOTE report for your reconstruction. To find it, look through the files on the left (📁) and download the file ending in `.html`. You can open it with any web browser and each section has a detailed explanation of what was tested and why. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4AjEfZYGptgr"
      },
      "source": [
        "# 💻 From Model to Fluxes\n",
        "\n",
        "Okay, now we have our metabolic reconstruction, so let's finally get to some actual metabolic modeling! We will make use of the [COBRAPY](https://opencobra.github.io/cobrapy/) package, which allows us to run a bunch of different analyses and simulations with our metabolic model.\n",
        "\n",
        "COPRAPY is already installed in your notebook. If you need help with COBRAPY in the future, you can check out [its documentation](https://cobrapy.readthedocs.io/en/latest/)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pDDbGUtprZUM"
      },
      "source": [
        "## Metabolic models and COBRAPY\n",
        "\n",
        "Before we actually start simulating the model let's have a look at it. We have learned that FBA builds on top of a lot of inferred information. So let's try to find this in the model we built. Here a short overview what ingredients we need to cook up an FBA result:\n",
        "\n",
        "![flux balance analysis](https://github.com/Gibbons-Lab/isb_course_2021/raw/main/docs/models/assets/fba.png)\n",
        "\n",
        "The first thing we have to to is to actually read the model into memory, and maybe inspect it a bit to see what is going on. **Before running the code cell, substitute `[YOUR ASSEMBLY]` with the ID of your assembly (including the `.asm1` part).** You can find the file name ending in .xml on the left (📁)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zm0N7BvXUthu"
      },
      "source": [
        "from cobra.io import read_sbml_model\n",
        "\n",
        "model = read_sbml_model(\"[YOUR ASSEMBLY].xml\")\n",
        "model"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MXHsl3zdsJv6"
      },
      "source": [
        "Okay, this gives us some very basic info about the model, like the number of metabolites and reactions, and the model objective, which in our case is growth. The model has accessions for metabolites and reactions, which can be used to look up a whole bunch of information and annotations. \n",
        "\n",
        "The objective is described, a bit opaquely, as growth. We talked a bit before about the biomass composition and we can have a look at this by taking a peak at the biomass reaction."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yzSjdmaAhsOa"
      },
      "source": [
        "model.reactions.Growth.reaction"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UOfC6dl9iL-E"
      },
      "source": [
        "This looks like what we'd expect. From the previous MEMOTE report, we know that the total consumed mass here should add up to something around 1gDW (± 0.05gDW). You can also see that the reactions are balanced. For instance, the consumed ATP is converted into ADP. \n",
        "\n",
        "The short names can be a bit cryptic. There are fairly human-readable here, but if we had used gapseq or ModelSEED, these features might have looked like `cpd00027_c`. Let's get some more info on, for instance, the 42nd metabolite in the model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "debRN0JgsC81"
      },
      "source": [
        "model.metabolites[42]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eP6Qa4lpjZn7"
      },
      "source": [
        "Okay, this is much more informative. You get the full name of the metabolite, the molecular formula, and the reactions in the model that include this metabolite. The same kind of output can be generated for individual reactions in the model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FuQ2SOPfsmoG"
      },
      "source": [
        "model.reactions[42]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WJeyrFUxs-4l"
      },
      "source": [
        "This gives us a lot of information about the reaction, including the substrates, products, reversibility, and bounds. The GPR is the *Gene-Protein-Reaction* rule and it states under what conditions the enzyme is functional.\n",
        "\n",
        "Do those `k*_*_*` look familiar? Exactly, those come from our assembly and gene predicton output. For instance, `k79_237_3` indicates the 3rd gene from the contig named k79_237. So we can track our contigs and genes all the way through the metabolic model.\n",
        "\n",
        "If you want more information about a reaction you can also check out the annotations (see below)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GkFvkFNJsyKE"
      },
      "source": [
        "model.reactions[42].annotation"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xavKRjIRtJtj"
      },
      "source": [
        "This will give you feature IDs for other databases, if available, which can be super-helpful when trying to map results back to experimental data. This is why the presence of sufficient annotations is a major criterion for MEMOTE scoring. \n",
        "\n",
        "## Flux Balance Analysis\n",
        "\n",
        "Okay, can we finally get to an actual simulation, please? 🥺 Yes! Running an FBA here is really simple. Just call `optimize()` on your model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ndbuS4ohtbXB"
      },
      "source": [
        "fba_solution = model.optimize()\n",
        "fba_solution"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ei3iO5eeuiPA"
      },
      "source": [
        "There is a lot of stuff going on under the hood here, so it's kind of cool that it runs so fast. \n",
        "\n",
        "Here you see the predicted growth rate and the fluxes. We won't spend too much time on the reduced costs here, but they are metrics that can be helpful sometimes. The two things to know about reduced costs are: \n",
        "\n",
        "1. They are only different from zero if the flux is zero.\n",
        "2. They quantify by how much the growth rate would be reduced if we force a flux through this reaction.\n",
        "\n",
        "So they are kind of the cost in growth rate the organism would have to pay to get flux through that previously shut-off reaction.\n",
        "\n",
        "This output may be a bit complex. For any solution we can also get a summary output:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UA7rK6E0uVXg"
      },
      "source": [
        "model.summary(fba_solution)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AG1NNm9PxAs3"
      },
      "source": [
        "This is still fairly long, but definitely looks nicer. We now see what your bacterium eats (uptake) and what proportion of the total carbon influx comes from that particular uptake (C-Flux). We get the same output for secretion. \n",
        "\n",
        "What is the major carbon source in for your organism in this simulation? \n",
        "\n",
        "Does the organism secrete carbon-rich metabolites? If yes, which ones?\n",
        "\n",
        "<br>\n",
        "\n",
        "---\n",
        "\n",
        "**Break time ☕  🥪** \n",
        "\n",
        "---\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vmfv6xayBVxD"
      },
      "source": [
        "# 🌍 Comparative Reverse Ecology\n",
        "\n",
        "Why are we interested in those uptake fluxes? Well, they are the recipe for an organism's growth, and thus, related to its evolutionary history (i.e. to the environmental conditions that it has adapted to live in). They are not solely what is present in the organism's surroundings, but the actual consumption preferences of that organism within a particular environment.\n",
        "\n",
        "When we started constructing our models, we had a description of the (metabolic) environment (the available nutrients) and our organism's genome. But now, suddenly, we learn something about how our organism interacts with the environment. And that means we learned something about its **ecology** 🤯\n",
        "\n",
        "Inferring information about an organism's ecology from genomic data is called [Reverse Ecology](https://en.wikipedia.org/wiki/Reverse_ecology). Check these [slides](https://gibbons-lab.github.io/isb_course_2021/models/#/26) for some more discussion on this point.\n",
        "\n",
        "Before we get there, however, we need to curate our model a bit..."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pvdGPq8uz-_t"
      },
      "source": [
        "## Getting realistic: PFBA and media\n",
        "\n",
        "As we noted, the particular flux solutions are often not unique and there may be many other solutions. To alleviate this let's use parsimonious FBA (pFBA), which minimizes the total import flux and the amount of enzymes required for growth.\n",
        "\n",
        "![parsimonious FBA](https://github.com/Gibbons-Lab/isb_course_2021/raw/main/docs/models/assets/pfba.png)\n",
        "\n",
        "pFBA is part of COBRAPY and will give us a new solution."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dc-zLuH92Pju"
      },
      "source": [
        "from cobra.flux_analysis import pfba\n",
        "\n",
        "pfba_solution = pfba(model)\n",
        "model.summary(pfba_solution)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rmTF8uIZ1k9z"
      },
      "source": [
        "We got our first real results -- awesome! But are they realistic? Remember that we stated that the biomass flux is essentially equivalent to the growth rate. We can convert the growth rate $\\mu$ to the doubling time with the formula\n",
        "\n",
        "$$\n",
        "\\tau = \\ln(2) / \\mu\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YsVweGSl0JN3"
      },
      "source": [
        "from math import log\n",
        "\n",
        "tau = log(2) / fba_solution.objective_value\n",
        "print(f\"The predicted dubling time is {tau:.2g} hours ({tau*60*60:.2g} seconds).\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UFfgwEuc1b9s"
      },
      "source": [
        "That probably looks a bit fast, right? So if you complained that this growth rate was too high in the Slack chat (or to yourself when we ran the first FBA), then you should take a few seconds to find the mirror closest to you and tell yourself that you are, indeed, a smart cookie, because you were absolutely right! No bacterium doubles that fast.\n",
        "\n",
        "But why is it predicted to grow so fast? \n",
        "\n",
        "## The Environment and Nutrient Preferences\n",
        "\n",
        "Well, we did not impose particular environmental conditions yet. So what did the model use? Let's look at the first uptake reaction."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uXcXT-831KOX"
      },
      "source": [
        "model.exchanges[0]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QrVULjR623kW"
      },
      "source": [
        "You see that the reaction is formulated into the export reaction (i.e. positive values denote the rate at which the compounds disappears). So, a negative flux denotes import and the lower bound dictates the largest possible import, which is 1000 mmol/(gDW h), in this case. So one mole per hour per gram of bacteria. So about $6\\cdot 10^{23}$ molecules per hour, hmmmmmm 🤨\n",
        "\n",
        "But we gapfilled the model with realistic media! So let's use those this time. Since the gut environment is probably a bit more interesting let's use this one. We will read in the environmental conditions using the `pandas` package and then format it a bit."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gBBIMkBs2ygl"
      },
      "source": [
        "import pandas as pd\n",
        "\n",
        "media_db = pd.read_csv(\"carveme_media.tsv\", sep=\"\\t\")\n",
        "media_db.index = media_db.reaction\n",
        "western_diet = media_db[media_db.medium == \"western_diet_gut\"]\n",
        "western_diet"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UsKO8vZN6he_"
      },
      "source": [
        "Here, the `flux` columns denotes the upper bound for the uptake flux. You can see that those are much smaller than 1000 mmol/(gDW h).\n",
        "\n",
        "Let's apply this medium to our model. It's easy! We just set the model medium to the flux bounds. However, some of the metabolites or transport reactions may not be present in the model. So we will first filter it by the exchange reactions that can be found in the model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GY4pywaI6d-s"
      },
      "source": [
        "exchange_ids = [r.id for r in model.exchanges]\n",
        "western_diet_found = western_diet[western_diet.reaction.isin(exchange_ids)]\n",
        "\n",
        "model.medium = western_diet_found.flux"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bhkbV9Ml7SvF"
      },
      "source": [
        "Now let's run the PFBA again."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NXjiBadI61fd"
      },
      "source": [
        "diet_sol = pfba(model)\n",
        "model.summary(diet_sol)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6PbeEhFc7yly"
      },
      "source": [
        "This should have changed the doubling time."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7QP0cBOp74e7"
      },
      "source": [
        "tau = log(2) / diet_sol.fluxes[\"Growth\"]\n",
        "print(f\"The new doubling time is {tau:.2g} h.\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iz-stBB88jAP"
      },
      "source": [
        "That changed quite a bit, didn't it? What are the major carbon sources this time? \n",
        "\n",
        "We can also subet the fluxes so we get only the uptake fluxes.\n",
        "\n",
        "**Remember:** In our model, uptakes are exchange reactions with negative flux values. \n",
        "\n",
        "In practice the solver will not be able to distinguish very small values from zero (due to [floating-point arithmetic](https://docs.python.org/3/tutorial/floatingpoint.html)). Therefore, we will consider uptakes to be 'active' if they are smaller than -1e-6. \n",
        "\n",
        "Calling `head()` will only show the first few fluxes."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cDs7B4Zw8Fu3"
      },
      "source": [
        "uptake = diet_sol.fluxes[diet_sol.fluxes.index.isin(exchange_ids) & (diet_sol.fluxes < -1e-6)]\n",
        "uptake.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PbLgcZrI-NS5"
      },
      "source": [
        "This looks interesting and quantifies the expected uptake of metabolites in the gut to achieve the optimal growth rate. These import values are actual proportions of consumed nutrients, which may be quite different for each individual organism in the gut. \n",
        "\n",
        "Let's make this a clean table and also add in the human-readable metabolite names."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ruWOa8gX_Q8T"
      },
      "source": [
        "uptake_nice = uptake.to_frame()\n",
        "uptake_nice[\"reaction\"] = uptake.index\n",
        "uptake_nice[\"metabolite\"] = [model.reactions.get_by_id(rid).reactants[0].name for rid in uptake.index]\n",
        "uptake_nice.sort_values(by=\"fluxes\").head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B5aTWeser4pE"
      },
      "source": [
        "## Comparative Reverse Ecology\n",
        "\n",
        "Now, let's compare all the uptake and secretion rates that everyone in the class has estimated. For that, let's just collect everybody's exchange results across the entire course. Easy peasy."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "h0q_O8DIC9Wv"
      },
      "source": [
        "from course import gimme_exchanges\n",
        "\n",
        "exchanges = gimme_exchanges()\n",
        "exchanges"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XWgWTC3WUu4L"
      },
      "source": [
        "Note that this includes uptake and secretion fluxes. As discussed before, transport reactions are formulated in the export direction, so negative fluxes denote uptake and positive fluxes denote secretion.\n",
        "\n",
        "Let's visualize this. We can start by drawing a heatmap of all the uptake fluxes. For this we will filter the table by uptake and use the absolute values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i6036zmWR349"
      },
      "source": [
        "uptake = exchanges[exchanges.flux < 0].copy()\n",
        "uptake.flux = uptake.flux.abs()\n",
        "uptake"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OB_lRhfzV-eg"
      },
      "source": [
        "Now let's convert this output into a matrix that we can visualize as a heatmap. One issue is that flux values between different metabolites can vary across many orders of magnitude. For instance water (H₂O) can often just diffuse passively into the cell with very large fluxes, but importing large molecules will always be much slower. There are many way to normalize these fluxes. For instance, one could scale it by molecular weight of a compound. However, a very easy way is to just log-transform the fluxes. We will stick with that for now."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YkrCt-QVR8_n"
      },
      "source": [
        "from seaborn import clustermap\n",
        "from numpy import log10\n",
        "\n",
        "uptake_matrix = uptake.pivot_table(\n",
        "    index=\"metabolite\", \n",
        "    columns=\"assembly\", \n",
        "    values=\"flux\", \n",
        "    fill_value=1e-6\n",
        ")\n",
        "clustermap(\n",
        "    log10(uptake_matrix), \n",
        "    cmap=\"magma\", \n",
        "    xticklabels=False, \n",
        "    yticklabels=True, \n",
        "    figsize=(20, 20)\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LhDLdHGzX3Hh"
      },
      "source": [
        "There are several things we can learn from this heatmap. Rows denote individual metabolites and columns denote individual isolates. The column dendrogram groups organisms with similar uptake fluxes together (i.e. organisms with similar metabolic requirements) and the row dendrogram groups metabolites that behave similarly across organisms (i.e. metabolites that tend to be taken up together).\n",
        "\n",
        "\n",
        "You can see that even though all isolates were subjected to the same environmental conditions they use very different strategies to produce biomass. Rows that are very \"patchy\" denote metabolites that are only necessary for a small number of isolates, while homogeneously colored rows denote metabolites that are required by all isolates to form biomass.\n",
        "\n",
        "One way to interpret the set of uptake fluxes for a single organism is as the realized growth niche of that organism in the gut of a human consuming an average European diet. One thing we could expect is that this niche is fairly conserved within a species. In ecology this is known as [*phylogenetic niche conservatism*](https://en.wikipedia.org/wiki/Phylogenetic_niche_conservatism). To investigate this phenomenon, let's have a look at *Phocaeicola vulgatus*, which is the recently renamed `Bacteroides vulgatus`, one of the most abundant taxa in the human gut. We will rotate the heatmap into a horizontal layout this time, so that the rows now denote isolates and columns denote metabolites. \n",
        "\n",
        "For this, we first generate a table that contains the annotatons for all isolates with the same order as our uptake matrix. We can then use this to filter the matrix by specific taxa."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "h-qgt4yQYAgu"
      },
      "source": [
        "isolates = uptake[[\"assembly\", \"genus\", \"species\", \"classification\"]].drop_duplicates()\n",
        "isolates.index = isolates.assembly\n",
        "isolates = isolates.loc[uptake_matrix.columns, :]\n",
        "\n",
        "phocaeicola = uptake_matrix.loc[:, isolates.species == \"Phocaeicola vulgatus\"]\n",
        "clustermap(\n",
        "    log10(phocaeicola.T), \n",
        "    cmap=\"magma\", \n",
        "    xticklabels=True, \n",
        "    yticklabels=True, \n",
        "    figsize=(22, 10)\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aVjYM2nNzZfR"
      },
      "source": [
        "At first glance, we see that the majority of uptakes seem to be identical for *Phocaeicola vulgatus* isolates. However, there are still some differences, indicating that an organisms' exact niche preferences can be strain-specific. \n",
        "\n",
        "\n",
        "\n",
        "But this is only one species in our data set. It would be much more interesting to look at all of them, so can we do better?\n",
        "\n",
        "## Niche maps\n",
        "\n",
        "Probably, but we might have to drop a bit of information in the visualization. If we move away from individual metabolites we can reduce our question to the following:\n",
        "\n",
        "> How close are organismal growth niches to one another?\n",
        "\n",
        "Just like a map can help you find the nearest grocery store to your house, we can build a map of niche distances for our growth utpake fluxes to identify organisms that are closest to one another in niche space. Thinking back to the course from yesterday, we already talked a bit about how to visualize multi-dimensional beta-diversity distances. \n",
        "\n",
        "We will employ a similar approach today, using a method called UMAP (\"map\" is even in the name). This is a non-linear approach that will try to estimate the best possible arrangement of uptake profiles in two dimensions (2D) that is consistent with actual multi-dimensional niche distances. Let's reduce each uptake profile to two coordinates in 2D space (like GPS coordinates). UMAP has a few parameters that can alter the visualization which are [explained here](https://umap-learn.readthedocs.io/en/latest/parameters.html). We will only adjust the minimum distance between points in the final projection because we already know that many uptake profiles are very similar and we don't want them to be piled on top of each other."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rfXqD1O3aBuo"
      },
      "source": [
        "from umap import UMAP\n",
        "\n",
        "results = UMAP(min_dist=1).fit_transform(log10(uptake_matrix.T))\n",
        "reduced = pd.DataFrame({\"UMAP 1\": results[:, 0], \"UMAP 2\": results[:, 1]}, index=uptake_matrix.columns)\n",
        "reduced"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "w0ZVAp9g1-5Z"
      },
      "source": [
        "Now we will merge this with the taxonomy information from the assemblies/isolates."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LruTiRqaaKGC"
      },
      "source": [
        "reduced = pd.concat([isolates, reduced], axis=1)\n",
        "reduced"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bR-Tx9yF2Gk3"
      },
      "source": [
        "This output can now be visualized on a 2D map and we can color the points using some additional information. For example, let's use the genus annotations of the isolates."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OUTWni3Mas3Z"
      },
      "source": [
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.figure(figsize=(15, 15))  # <-- this makes the figure larger\n",
        "\n",
        "reduced[\"genus\"] = reduced.species.str.split(\" \").str[0]\n",
        "sns.scatterplot(\n",
        "    x=\"UMAP 1\", \n",
        "    y=\"UMAP 2\", \n",
        "    hue=\"genus\", \n",
        "    data=reduced\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VOip4EEq3IYz"
      },
      "source": [
        "So we see that there is definitely a tendency for niches to be similar within a genus. But not always. Let's have a look at the genus *Bacteroides*, for instance."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jEka6Jnkl8uB"
      },
      "source": [
        "plt.figure(figsize=(10, 10))\n",
        "sns.scatterplot(\n",
        "    x=\"UMAP 1\", \n",
        "    y=\"UMAP 2\", \\\n",
        "    hue=\"species\", \n",
        "    data=reduced[reduced.genus == \"Bacteroides\"]\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9GOB7RyI3udI"
      },
      "source": [
        "As you can see there is a pretty discrete set of growth niches, but they do not separate strictly by species. Looking at this plot, what would be your prediction for the coexistence of these species in a single individual's gut?\n",
        "\n",
        "**Hint:** What happens if you and four of your friends all want the last chocolate chip cookie?\n",
        "\n",
        "**However**, there is a major drawback here that we need to consider. What could it be?\n",
        "\n",
        "# Exercise"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oNqAdUdB5Iy2"
      },
      "source": [
        "## Niche variability\n",
        "\n",
        "Let's finish by putting some hard numbers on the degree of niche conservatism across taxa. On the first day of the course we already learned how to run statistical tests on differences in beta-diversity distances between groups of samples: PERMANOVA. We can't rely on Qiime 2 here, so we will use [scikit-bio](http://scikit-bio.org/). \n",
        "\n",
        "We will start by actually calculating the distances between uptake profiles for each isolate (i.e. the equivalent to beta-diversity in metabolic niche space). here we use the `pdist` function from [scipy](https://www.scipy.org/) which uses Euclidean distances by default. So it's probably best to log-transform again."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YG1FReX95oSC"
      },
      "source": [
        "from scipy.spatial.distance import pdist\n",
        "from skbio.stats.distance import DistanceMatrix\n",
        "\n",
        "isolates.index = isolates.assembly\n",
        "D = DistanceMatrix(\n",
        "    pdist(log10(uptake_matrix).T), \n",
        "    ids=uptake_matrix.columns\n",
        ")\n",
        "D.shape"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-1NDkKWE6awX"
      },
      "source": [
        "Now we can run our PERMANOVA. Let's test whether the genus explains a significant fraction of variability in the growth niches. The `scikit-bio` PERMANOVA does not calculate the explained variance (R²) by default, so we will add it in here with a helper function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RW4XKU3f533s"
      },
      "source": [
        "from skbio.stats.distance import permanova\n",
        "from course import add_r_squared\n",
        "\n",
        "perm_results = permanova(D, grouping=isolates, column=\"genus\")\n",
        "add_r_squared(perm_results)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wTrVqTgi84Qk"
      },
      "source": [
        "How do you interpret the result? Also substitute the tested rank to species by changing the `column=\"...\"` part in the code above. \n",
        "\n",
        "- What do you observe for the R²? \n",
        "- How do you intepret the results when comparing the genus and species level?\n",
        "- What is your interpretation of the unexplained variance 1-R²? "
      ]
    }
  ]
}