workshop.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "workshop.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Gibbons-Lab/mmmih/blob/main/workshop.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C0vqP4LJ9y6K"
      },
      "source": [
        "# 🧫🦠 Modeling microbiota-wide metabolism with MICOM\n",
        "\n",
        "This notebook will accompany the 3rd Microbiome Movement - Maternal and Infant Health workshop. The presentation slides can be [found here](https://gibbons-lab.github.io/mmmih). \n",
        "\n",
        "You can save your own local copy of this notebook by using `File > Save a copy in Drive`. You may be prompted to cetify the notebook is safe. We promise that it is 🤞\n",
        "\n",
        "**Disclaimer:**\n",
        "The linear and quadratic programming problems MICOM has to solve are very large and very complicated. There are some very good commercial solvers that are very expensive (even though they are often free for academic use). To make this tutorial as accessible as possible we will use the Open Source solver [OSQP](https://osqp.org/). OSQP is amazing with quadratic programming problems (kudos!) but not as accurate for linear problems. Solvers usually only guarantee a solution within a certain numerical tolerance of the real solution. In order to make everything work with OSQP this tolerance has to be relaxed to about 10<sup>-3</sup>. This means that any result with an absolute value smaller than that might very well be zero so we should look at larger values only. Installing cost-free academic versions of commercial solvers like [IBM CPLEX](https://www.ibm.com/analytics/cplex-optimizer) or [Gurobi](https://www.gurobi.com/) would allow you to lower the tolerance to 10<sup>-6</sup>.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qgBBl4GtuTuX"
      },
      "source": [
        "# 📝 Setup\n",
        "\n",
        "MICOM installation is is usually pretty straight-forward and can be as easy as typing `pip install micom` into your Terminal.\n",
        "\n",
        "But first let's start by downloading the materials again and switching to the folder."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ckON4xr3_bW5"
      },
      "source": [
        "!git clone https://github.com/gibbons-lab/mmmih materials\n",
        "%cd materials"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "is6fmBUeorwv"
      },
      "source": [
        "## Basic Installation\n",
        "\n",
        "Installing MICOM is straight-forward in Python."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M_TeC5yrst3h"
      },
      "source": [
        "!pip install -q micom\n",
        "\n",
        "print(\"Done! 🎉 \")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oPy1f-WLI0lZ"
      },
      "source": [
        "Okay, all done. So let's get started 😁.\n",
        "\n",
        "# 💻 MICOM\n",
        "\n",
        "We will use the Python interface to MICOM since it plays nicely with Colaboratory. However, you could run the same steps within the QIIME 2 MICOM plugin ([q2-micom](https://library.qiime2.org/plugins/q2-micom/26/)). \n",
        "\n",
        "Here is an overview of all the steps and functions:\n",
        "![micom overview](https://github.com/micom-dev/q2-micom/raw/706f583a060b91c12c0cec7acea2354fdd0dd320/docs/assets/overview.png).\n",
        "\n",
        "MICOM starts from a combined abundance/taxonomy table, which MICOM abbreviates to a taxonomy table. From that we usually generate a set of community models, one for each sample. \n",
        "\n",
        "Because we don't have that much time we will start from a set of already generated models. Those are personalized models for 8 individuals: 4 healthy ones and 4 with recurrent *C. difficile* infections (illustration by Stephanie Swegle). \n",
        "\n",
        "![recurrent C. difficile cycle](https://github.com/Gibbons-Lab/mmmih/raw/main/docs/assets/cycle.png)\n",
        "\n",
        "Those were processed with QIIME 2 and then matched to the AGORA database. The entire build process is described in more details in the [MICOM documentation](https://micom-dev.github.io/micom/high_level.html#Building-and-models-and-simulating-growth).\n",
        "\n",
        "To get a quick overview over the models we will start by looking at the manifest which is a summary of the build step.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UV9SObSQkSZh"
      },
      "source": [
        "import pandas as pd\n",
        "\n",
        "manifest = pd.read_csv(\"models/manifest.csv\")\n",
        "manifest"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AEk7yfd1lbYp"
      },
      "source": [
        "We get a summary of the sample IDs, the number of taxa that were found in each sample and how many could be matched to AGORA. We only included taxa that constitute at least 2.5% of the total abundance in each sample here to make everything a bit quicker.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y4KAJkhIdspQ"
      },
      "source": [
        "\n",
        "\n",
        "So we now have our community models and can leverage MICOM fully by simulating community growth.\n",
        "\n",
        "## Simulating growth\n",
        "\n",
        "With our community models built, we can start to simulate growth with the cooperative tradeoff algorithm. Because we have no diet information for our samples, we will apply the same 'average Western Diet' to each individual. We will start by downloading this diet from the [MICOM data repository](https://zenodo.org/record/3755182/)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Xi4hZLiaPLNu"
      },
      "source": [
        "!wget -O western_diet_gut.qza https://zenodo.org/record/3755182/files/western_diet_gut.qza?download=1"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TVWMEu4nfKdL"
      },
      "source": [
        "This is a QIIME 2 artifact, which we can load into MICOM."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3XBGQlFnfO4W"
      },
      "source": [
        "from micom.qiime_formats import load_qiime_medium\n",
        "\n",
        "medium = load_qiime_medium(\"western_diet_gut.qza\")\n",
        "medium"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4s8R4WYUez4g"
      },
      "source": [
        "Many dietary components get absorbed in the small intestine. This medium was created by taking dietary components and depleting all nutrients absorbed in the small intestine by a factor of 10. \n",
        "\n",
        "Okay let's go right ahead and simulate growth. This will take a little while and is a great time to ask some questions 🤔"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "l-vJJwtdRZxP"
      },
      "source": [
        "from micom.workflows import grow\n",
        "import pickle\n",
        "\n",
        "growth_results = grow(manifest, \"models\", medium, tradeoff=0.5, threads=2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PrM7t9u9fGiC"
      },
      "source": [
        "If that takes too long or was aborted, we can read it in from the treasure chest."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DTHObChXfIRT"
      },
      "source": [
        "# Will only run if the previous step failed\n",
        "\n",
        "import pickle\n",
        "\n",
        "try:\n",
        "  growth_results\n",
        "except NameError:\n",
        "  growth_results = pickle.load(open(\"treasure_chest/growth.pickle\", \"rb\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "69gq9QfAzqxq"
      },
      "source": [
        "What kind of results did we get? Well, `grow` returns a tuple of 3 data sets:\n",
        "\n",
        "1. The predicted growth rate for all taxa in all samples\n",
        "2. The import and export fluxes for each taxon and the external environment\n",
        "3. Annotations for the fluxes mapping to other databases\n",
        "\n",
        "The growth rates are pretty straightforward."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xHQKoxgRzqWM"
      },
      "source": [
        "growth_results.growth_rates.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F5BK7DDv0UfA"
      },
      "source": [
        "More interesting are the exchange fluxes."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lQW2BBS10jdN"
      },
      "source": [
        "growth_results.exchanges"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Pu5XtkUl1YG1"
      },
      "source": [
        "So we see how much of each metabolite is either consumed or produced by each taxon in each sample. `tolerance` denotes the accuracy of the solver and tells you the smallest absolute flux that is likely different form zero (i.e. substantial flux). *All of the fluxes are normalized to 1g dry weight of bacteria*. So, you can directly compare fluxes between taxa, even if they are present at very different abundances. \n",
        "\n",
        "However, the metabolite names may not be very informative. That's why we have our annotations! For instance, to figure out what `ac[e]` is (air conditioning?), we can do the following:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DphXa9hw1yxM"
      },
      "source": [
        "anns = growth_results.annotations\n",
        "anns[anns.metabolite == \"ac[e]\"]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CImtzqRJbbGj"
      },
      "source": [
        "Ohhh, it's acetate. Yeah that makes more sense 🕵️‍♀️. For the AGORA models you can also use the official VMH knowledge base at https://vmh.life maintained by Dr. Thiele's, lab which will give you rich information on metabolites and reactions. For instance, you can find out a lot more about acetate at: https://www.vmh.life/#metabolite/ac. \n",
        "\n",
        "# 📊 Visualizations\n",
        "\n",
        "Ok, we have seen that we generate a lot of output data from the growth simulations. But how do we make sense of it all? \n",
        "\n",
        "We will use the standard visualizations included in MICOM. These tools take in the growth results we obtained before and create visualizations in standalone HTML files that bundle the plots and raw data and can be viewed directly in your browser."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G1JbbKrLcVye"
      },
      "source": [
        "Normally, we could call `viz.view()` afterwards and it would open it in our web browser. However, this will not work in Colab. However, the plot function create the file `growth_rates_[DATE].html` in your `materials` folder. To open it simply download that file and view it in your browser. We can see that there are many things going on, but it's not super clear. Let's continue.\n",
        "\n",
        "## Growth niches\n",
        "\n",
        "Two really important questions are 'what dietary nutrients are consumed by the microbiota and what metabolites do the microbiota produce?' We provided nutrients in our medium, but we don't actually know yet what was eaten by the microbiota. Let's check that out using the `plot_exchanges_per_sample` function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5zh17r6xf15I"
      },
      "source": [
        "from micom.viz import *\n",
        "\n",
        "plot_exchanges_per_sample(growth_results)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qXnbUCCs2yVG"
      },
      "source": [
        "You can see that even though all samples have been subjected to the same environmental conditions there are some differences in the set of consumed metabolites. Why do you think that is?\n",
        "\n",
        "## Metabolic connections to a phenotype\n",
        "\n",
        "That is all nice, but how does that relate to recurrent *C. diff* infections? To answer that question we can use the `plot_fit` function. This will run a logistic regression with an L1 penalty on coefficients, with our diesase status as the response variable and the normalizeed fluxes as independent variables. In general, import fluxes are not as predictive because, well, they are more relevant to the bacteria than us. What we usually care about (from a host-health perspective) are the production fluxes of metabolites. These are the total production fluxes into the extracellular (lumenal) space, which includes the set of metabolites that are bioavailable to the host. \n",
        "\n",
        "We will start by reading the metadata annotations for our samples that tell us which samples belong to which health status."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "X9hqoO4go0h1"
      },
      "source": [
        "metadata = pd.read_table(\"metadata.tsv\").rename(columns={\"id\": \"sample_id\"})\n",
        "metadata.index = metadata.sample_id\n",
        "metadata"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nx8ZGEKOuKSu"
      },
      "source": [
        "\n",
        "Because OSQP has a somewhat lower solver accuracy, we will be conservative for what we consider to be 'substantial flux' and will filter out fluxes smaller than 0.001 mmol/l."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ebsw0XD-l2Lj"
      },
      "source": [
        "from micom.viz import *\n",
        "\n",
        "pl = plot_fit(growth_results, metadata.disease_stat, atol=1e-3, flux_type=\"production\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4Y_XfHkB4sO8"
      },
      "source": [
        "This will again create a file `fit_[DATE].html` that you can open. You will see the production fluxes most predictive of the phenotypes of interest and you can compare them across the groups. In the coefficient plot negative coefficients mean the particular metabolite is produced at a higher flux in the first group (recurrent C. diff in our case), whereas positive coefficients mean it is produced at a higher flux in the second group (healthy). You can plot the fluxes for the individual metabolites in the lower left. \n",
        "\n",
        "Which ones do you think can be trusted on which ones can not?\n",
        "\n",
        "Whatever we conclude here is based on *very* little data so it would be purely hypothetical. However, it may point us to future avenues to explore."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Sppu1VKEKohS"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}