micom.ipynb

{
  "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_2022/blob/main/micom.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 second session of the 2022 ISB Microbiome Course. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2022/micom). \n",
        "\n",
        "You can save your own local copy of this notebook by using `File > Save a copy in Drive`. You may be promted 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/), which is installed along with MICOM. 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 alow 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",
        "First let's start by downloading the materials again and switching to the folder."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ckON4xr3_bW5"
      },
      "outputs": [],
      "source": [
        "!git clone https://github.com/gibbons-lab/isb_course_2022 materials\n",
        "%cd materials"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "is6fmBUeorwv"
      },
      "source": [
        "## Basic Installation\n",
        "\n",
        "Installing MICOM is straight-forward in Python. OSQP itself will be installed automatically along with it."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "M_TeC5yrst3h"
      },
      "outputs": [],
      "source": [
        "!pip install -q micom\n",
        "\n",
        "print(\"Done! 🎉 \")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6oJrxxz6tV9T"
      },
      "source": [
        "## Enable QIIME 2 interactions\n",
        "\n",
        "Before we start, we also need to install packages to read the \"biom\" file format used by QIIME 2 to save tables. This is only necessary if you want to read QIIME 2 FeatureTable artifacts (like the ones we constructed yesterday)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5rZX7SK_toLp"
      },
      "outputs": [],
      "source": [
        "!pip install -q numpy Cython\n",
        "!pip install -q biom-format\n",
        "\n",
        "print(\"Done! 🎉 \")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oPy1f-WLI0lZ"
      },
      "source": [
        "Okay, all done. So let's get started building some models 🦺🛠d😁.\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 across both interfaces:\n",
        "![micom overview](https://github.com/micom-dev/q2-micom/raw/706f583a060b91c12c0cec7acea2354fdd0dd320/docs/assets/overview.png)\n",
        "\n",
        "The process of building a metabolic model in MICOM begins with constructing a combined abundance/taxonomy table, referred to hereafter as a taxonomy table. Let's load a sample taxonomy table to see what it looks like: \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UV9SObSQkSZh"
      },
      "outputs": [],
      "source": [
        "from micom.data import test_data\n",
        "\n",
        "test_data().head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AEk7yfd1lbYp"
      },
      "source": [
        "In this taxonomy table, we see four identical strains of _E. coli_ (1 through 4), across two samples (sample_1 and sample_2). We can see that each row represents a single taxon in a single sample, and the `abundance` column identifies the abundance of that taxon in the sample. \n",
        "\n",
        "The `id` column specifies identifiers for the taxa and should be expressive and not include spaces or special characters. Since we are using a taxonomy database to build our models (more on that soon), we don't need a `file` column.\n",
        "\n",
        "You might notice that this dataframe looks very different from what we generated in yesterday's tutorial, where we ended up with separate QIIME 2 artifacts 😱 \n",
        "\n",
        "No worries, we can deal with that.\n",
        "\n",
        "## Importing data from QIIME 2\n",
        "\n",
        "MICOM can read QIIME 2 artifacts. You don't even need to have QIIME 2 installed for that! But before we do so, let's resolve one issue. We discussed that MICOM summarizes genome-scale models into pangenome-scale models as a first step, but our data are on the ASV level...so how will we know what to summarize?\n",
        "\n",
        "Basically, a specific model database can be used to quickly summarize pangenome-scale models for use within MICOM. So, before we read our data we have to decide which model database to use. We will go with the [AGORA database](https://pubmed.ncbi.nlm.nih.gov/27893703/), which is a curated database of more than 800 bacterial strains that commonly live in the human gut. In particular, we will use a version of this database summarized on the genus rank which can be downloaded from the [MICOM data repository](https://doi.org/10.5281/zenodo.3755182), which contains a whole lot of prebuilt databases. This database is available from the materials folder that we previously cloned. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_57iya0D3L6-"
      },
      "source": [
        "Now we're all set to start building models! The data we previously collected can be found in the `treasure_chest` folder, so we can use those files to build our taxonomy for MICOM. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "o0vBAiiqqPLC"
      },
      "outputs": [],
      "source": [
        "from micom.taxonomy import qiime_to_micom\n",
        "\n",
        "tax = qiime_to_micom(\n",
        "    \"treasure_chest/dada2/table.qza\", \n",
        "    \"treasure_chest/taxa.qza\", \n",
        "    collapse_on=\"genus\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-TQ6Zp7wouk4"
      },
      "source": [
        "Notice the `collapse_on` argument. That will specify the rank on which to sumarize and can be a list of several ranks. When matching taxonomy you can either match by the particular rank of interest (for example, just comparing genus names here), or you could compare the entire taxonomy, which will require all taxonomic ranks prior to the target rank to match. For that you cloud specify `collapse_on=[\"kingdom\", \"phylum\", \"class\", \"order\", \"family\", \"genus\"]`. \n",
        "\n",
        "Taxonomic names will often not match 100% between databases. For instance, the genus name \"Prevotella\" in one database may be \"Prevotella_6\" in another. The more ranks you use for matching the more likely you are to run into these issues. However, the more taxonomic ranks you use to match the more confident you can be that your observed taxon really is the same taxon as the one in the model database.\n",
        "\n",
        "The resulting table will contain the same abundances but it will include more ranks if `collapse_on` is a list. All ranks present in the taxonomy will be used when matching to the database. We will stick with the \"lax\" option of only matching on genus ranks.\n",
        "\n",
        "Let's now take a look at the taxonomy table we generated:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2eS_g-ffNUAt"
      },
      "outputs": [],
      "source": [
        "tax"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R5zwLVFcNTbq"
      },
      "source": [
        "That looks more like the example! Again, we have a row for each taxon in each sample, so we're good to go. \n",
        "\n",
        "One helpful thing to do is to merge in our metadata, so we'll have it at hand for the following steps. In our case, the metadata will include the ethnic group, region, and subsistence type of each of the study participants. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X9hqoO4go0h1"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "\n",
        "metadata = pd.read_table(\"data/metadata.tsv\").rename(columns={\"id\": \"sample_id\"})\n",
        "tax = pd.merge(tax, metadata, on=\"sample_id\")\n",
        "tax"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RTbYBR8cJfup"
      },
      "source": [
        "With our taxonomy table ready to go, and our metadata merged, its finally time to get to the model building! 🎉\n",
        "\n",
        "## Building community models\n",
        "\n",
        "With the data we have now, building our models is pretty easy. We just pass our taxonomy table and model database to MICOM. We will remove all taxa that make up less than 1% of the community to keep the models small and speed up this tutorial. We will also have to specify where to write the models. For simplicity, we'll run this process in parallel over two threads. It should take around 10 minutes to finish."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "kDbSN71SmCZr"
      },
      "outputs": [],
      "source": [
        "from micom.workflows import build\n",
        "from micom import Community\n",
        "import pandas as pd\n",
        "\n",
        "manifest = build(tax, \"agora103_genus.qza\", \"models\", solver=\"osqp\", \n",
        "                 cutoff=0.01, threads=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kwya6vbZZSUo"
      },
      "source": [
        "You'll see a warning pop up indicating that less than 50% of the abundances can be matched to the database for one of the samples. This can happen with some data, and may indicate the models may not be completely representative of the samples.Typically a fraction of 80% or more is considered great. We'll continue, but remember to keep an eye out for this in future projects!\n",
        "\n",
        "In lower-biomass 16S amplicon sequencing samples from stool, many reads can match to food components or to host mitochondria and these hits probably do not contribute much to bacterial community metabolism. These hits will be excluded from MICOM. \n",
        "\n",
        "Let's take a look what we got back from the `build` process."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "r9qwglr88Ise"
      },
      "outputs": [],
      "source": [
        "manifest"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y4KAJkhIdspQ"
      },
      "source": [
        "This will tell you how many taxa were found in the database and what fraction of the total abundance was represented by the database. For most samples, this looks okay (i.e., >70% of abundance represented). \n",
        "\n",
        "So we now have our community models and can leverage MICOM fully by simulating community growth - let's discuss what we want to look at. \n",
        "\n",
        "### Dietary Context\n",
        "\n",
        "Now that our models are ready to go, let's think about some of the insights we might gain from these samples. First and foremost, we can investigate the metabolomic response of the gut microbiome in a set of individuals from different underrepresented populations - in these different groups, do the microbes fall into different niche spaces?\n",
        "\n",
        "Additionally, we can use MICOM to take a mechanistic look at how the microbiome of these individuals changes based on __dietary input__. That is, if we switch the diet of individuals from indigenous populations from a matched diet (i.e., a diet based on a typical meal eaten by an individual from that population) to an unmatched diet (e.g., the average Austrian diet), how does the metabolomic profile of the microbiome shift?\n",
        "\n",
        "All that and more, coming up. Stay tuned!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ENR_ptvrTtp_"
      },
      "source": [
        "First, we'll download the matched diets for each of our communities that we constructed by hand. Since this is a .csv file, we can read it using pandas."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7099I-AuUE8J"
      },
      "outputs": [],
      "source": [
        "matched_medium = pd.read_csv('data/per_sample_media.csv')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "etrvjwBLkKdR"
      },
      "source": [
        "You'll notice that this dataframe is pretty long - that's because this medium is \"per-sample\", meaning each of the nine models has it's own medium represented in the file, dependent on the traditional food of the community they live in. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4s8R4WYUez4g"
      },
      "source": [
        "### Growing the models\n",
        "Great, now we have our media & our models, it's time to get growing. We'll grow each of the nine samples using the per-sample matched medium. This will take some time, so we'll use that time as an opportunity to discuss more in depth what these processes do, and what to look for in the results. First, let's run the `grow()` command using the matched medium. This will take the models we've built, and find an optimal solution to the fluxes based upon the medium that's been applied."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3WH8VVrVS4mv"
      },
      "source": [
        "If that takes too long or was aborted, we can read it in from the treasure chest."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IjDguZEcWGjG"
      },
      "outputs": [],
      "source": [
        "from micom.workflows import grow, save_results\n",
        "\n",
        "matched_growth = grow(manifest, \"models\", matched_medium, tradeoff=0.5, threads=2)\n",
        "\n",
        "# We'll save the results to a file\n",
        "save_results(matched_growth, \"treasure_chest/matched_growth.zip\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rHedHJxHWkjy"
      },
      "source": [
        "Again, if that takes too long or was aborted, we can read it in from the treasure chest."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "rcPNBkDpWGrQ"
      },
      "outputs": [],
      "source": [
        "from micom.workflows import load_results\n",
        "\n",
        "try:               # Will only run if the previous step failed\n",
        "  matched_growth\n",
        "except NameError:\n",
        "  matched_growth = load_results(\"treasure_chest/matched_growth.zip\")"
      ]
    },
    {
      "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",
        "### 📈 Growth Rates\n",
        "\n",
        "The growth rates are pretty straightforward."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lIz36dSQWyO4"
      },
      "outputs": [],
      "source": [
        "matched_growth.growth_rates"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F5BK7DDv0UfA"
      },
      "source": [
        "### ↔️ Exchange Fluxes \n",
        "\n",
        "More interesting are the exchange fluxes. These reactions represent the import and export of metabolites from the system Let's look at those now:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lQW2BBS10jdN"
      },
      "outputs": [],
      "source": [
        "matched_growth.exchanges"
      ]
    },
    {
      "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",
        "If you're curious what the abbreviation for each of these metabolites represents, that can be found in the annotations dataframe. For instance, let's find out what `\"but[e]\"` represents. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DphXa9hw1yxM"
      },
      "outputs": [],
      "source": [
        "anns = matched_growth.annotations\n",
        "anns[anns.metabolite == \"but[e]\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GVHLD2dm4a6B"
      },
      "source": [
        "Butyrate! Interesting, [that's an important metabolite](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8608412/)! All of these annotations and more information at are also available at https://vmh.life, maintained by Dr. Ines Thiele's lab. \n",
        "\n",
        "We want to compare these results from our matched medium against an unmatched medium, to illustrate how important it is to use the proper dietary inputs for the modeling process. Rather than repeat all the steps above, we can load in the results from the treasure chest. The following cell will load the growth results for the same models as above, but with an unmatched medium applied. In this case, the unmatched medium is an average Austrian diet. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "I14SxizN45ya"
      },
      "outputs": [],
      "source": [
        "from micom.workflows import load_results\n",
        "\n",
        "unmatched_growth = load_results(\"treasure_chest/unmatched_growth.zip\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CImtzqRJbbGj"
      },
      "source": [
        "\n",
        "# 📊 Visualizations\n",
        "\n",
        "Let's visualize our results. Because of the rich output of these models, it can be overwhelming to represent it all, but don't worry! There are tools in place for this already. \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.\n",
        "\n",
        "First, let's look at the growth rates of each taxon across samples. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "EaplMHFLcMT7"
      },
      "outputs": [],
      "source": [
        "from micom.viz import *\n",
        "\n",
        "viz_unmatched = plot_growth(unmatched_growth, filename=\"unmatched_growthrates.html\")\n",
        "viz_matched = plot_growth(matched_growth, filename=\"matched_growthrates.html\")"
      ]
    },
    {
      "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 creates the file `growth_rates_[DATE].html` in your `materials` folder. To open it, simply download that file and view it in your web 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",
        "Another thing we can look at is whether individual taxa inhabit different growth niches across different dietary contexts. Here we can use the `plot_exchanges_per_taxon` function to see how exchanges differ within and between taxa, within and across human populations. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NlZrfv38esj8"
      },
      "outputs": [],
      "source": [
        "plot_exchanges_per_taxon(unmatched_growth, perplexity=4, direction=\"import\", filename=\"unmatched_niche.html\")\n",
        "plot_exchanges_per_taxon(matched_growth, perplexity=4, direction=\"import\", filename=\"matched_niche.html\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qXnbUCCs2yVG"
      },
      "source": [
        "\n",
        "This function projects the full set of import or export fluxes onto a two dimensional plane, and arranges taxa so that more similar flux patterns lie nearer together. Taxa closer to one another compete for a more similar set of resources (and/or produce a more similar set of metabolites). The center of the plot signifies a more competitive nutrient space, whereas clusters on the outskirts denote more isolated niches.\n",
        "\n",
        "You can tune [TSNE parameters](https://distill.pub/2016/misread-tsne/), such as perplexity, to get a more meaningful grouping. We will lower the perplexity here since we don't have a lot of data points.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4Y_XfHkB4sO8"
      },
      "source": [
        "## Comparative Metabolomics\n",
        "\n",
        "Now let's compare the metabolomic exports between the two dietary contexts. We're interested to see how the metabolomic profile of the microbiome changes when the diet changes, as changes in diet can lead to changes in host health.  To look into this deeper, we'll concatenate the exchange results into one dataframe, transform the data and then plot the metabolite exports on a heatmap. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pAtec3I78DJJ"
      },
      "source": [
        "We can use the `production_rates` function in MICOM to calculate production rates from the growth results. We'll also run a centered log ratio transformation on the data, to account for the compositional nature of these data and compare all the fluxes against each other. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DJAsxjnFdajN"
      },
      "outputs": [],
      "source": [
        "from micom.measures import production_rates\n",
        "import numpy as np\n",
        "\n",
        "\n",
        "matched_prod = production_rates(matched_growth)\n",
        "matched_prod['diet'] = 'matched'\n",
        "unmatched_prod = production_rates(unmatched_growth)\n",
        "unmatched_prod['diet'] = 'unmatched'\n",
        "exchanges = pd.concat([matched_prod, unmatched_prod])  # merge the production rates\n",
        "\n",
        "exchanges = pd.pivot_table(                            # convert to a matrix of samples vs. metabolites\n",
        "    exchanges,                                         # that contains the production rates\n",
        "    index = ['diet','sample_id'], \n",
        "    columns = 'name', \n",
        "    values = 'flux'\n",
        ")\n",
        "exchanges = exchanges.T.fillna(0.0)                    # if a metabolite is not produced its flux is zero\n",
        "exchanges = exchanges.apply(                           # ...and a CLR transform again, normalizes the fluxes\n",
        "    lambda xs: np.log(xs + 0.001) - np.log(xs.mean() + 0.001),\n",
        "    axis=0)\n",
        "exchanges = exchanges.reindex(                         # sort by variance, highest variance fluxes first\n",
        "    exchanges.var(axis = 1).sort_values(ascending=False).index\n",
        ")\n",
        "\n",
        "exchanges"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nYuc7Wu38nYd"
      },
      "source": [
        "We can use seaborn to plot our heatmap:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YHnKeFuF3qAt"
      },
      "outputs": [],
      "source": [
        "import seaborn as sns \n",
        "import numpy as np\n",
        "\n",
        "sns.clustermap(\n",
        "    exchanges.head(50),  # take 50 highest fluxes\n",
        "    cmap = 'viridis',\n",
        "    yticklabels = True,  # show all metabolite names\n",
        "    figsize = (8, 12)    # size of the heatmap\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b1i2LStZiUst"
      },
      "source": [
        "We can see here that the dietary context applied to the metabolic models is important - there are significant differences in export fluxes between the matched and unmatched diets. Overall, we can see clustering of the matched and unmatched groups, and more defined patterns when using the matched diet as opposed to the unmatched diet. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l8maOr3w2bOo"
      },
      "source": [
        "# 🏫 Exercises"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TqdappAU20oT"
      },
      "source": [
        "Up to now, we have mostly used MICOM's \"high-level\" API, which is designed for working with several samples in parallel. However, MICOM also allows you to work with single models. We will choose a single sample now for further analysis.\n",
        "\n",
        "First, let's recall what samples we had: "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bllnCa1O28YZ"
      },
      "outputs": [],
      "source": [
        "manifest"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VjUustQG2_bX"
      },
      "source": [
        "## Diet-Specific Interactions\n",
        "We can use the `load_pickle` funciton in MICOM to load one of the models we're interested in, and look at differences in growth between the two diets we have available. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uh7QFi9V3HIX"
      },
      "outputs": [],
      "source": [
        "from micom import load_pickle\n",
        "from micom.qiime_formats import load_qiime_medium\n",
        "\n",
        "matched_medium = load_qiime_medium(\"data/himalaya.qza\")\n",
        "\n",
        "com = load_pickle(\"models/chepang1.pickle\")\n",
        "com.medium = matched_medium.flux\n",
        "sol = com.cooperative_tradeoff(fraction=0.5, fluxes=True, pfba=True)\n",
        "sol"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mkpqy7yz7cJV"
      },
      "source": [
        "Go ahead and try growing this model using the other available medium (i.e., the average Austrian diet):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8G655l0W7Zw-"
      },
      "outputs": [],
      "source": [
        "# Your code here"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5IlJ-YpUCQzp"
      },
      "outputs": [],
      "source": [
        "# Your code here"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CsqIRTbC7doD"
      },
      "source": [
        "# 🔵 Addendum\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hycoXNTi5xsH"
      },
      "source": [
        "## Choosing a tradeoff value\n",
        "\n",
        "Even if you don't have growth rates available you can still use your data to choose a decent tradeoff value. This can be done by choosing the largest tradeoff value that still allows growth for the majority of the taxa that you observed in the sample (if they are present at an appreciable abundance, they should be able to grow). This can be done with the `tradeoff` workflow in MICOM that will run cooperative tradeoff with varying tradeoff values, which can be visualized with the `plot_tradeoff` function."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8_1jesZTHYra"
      },
      "outputs": [],
      "source": [
        "from micom.workflows import tradeoff\n",
        "import micom\n",
        "\n",
        "tradeoff_results = tradeoff(manifest, \"models\", matched_medium, threads=2)\n",
        "tradeoff_results.to_csv(\"tradeoff.csv\", index=False)\n",
        "\n",
        "plot_tradeoff(tradeoff_results, tolerance=1e-4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y9703vhK6d6c"
      },
      "source": [
        "After opeing `tradeoff_[DATE].html` you will see that, for our example here, all tradeoff values work great. This is because we modeled very few taxa, which keeps the compettion down. If you would allow for fewer abundant taxa in the models, this would change drastically. For instance, here is an example from a colorectal cancer data set:\n",
        "\n",
        "[![tradeoff example](https://micom-dev.github.io/micom/_images/tradeoff.png)](https://micom-dev.github.io/micom/_static/tradeoff.html)\n",
        "\n",
        "You can see how not using the cooperative tradeoff would give you nonsense results where only 10% of all observed taxa grew. A tradeoff value of 0.6-0.8 would probably be a good choice for this particular data set."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3.10.7 ('micom')",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.10.7"
    },
    "vscode": {
      "interpreter": {
        "hash": "c991a7ed881363492957ff225bb30af9d5174cd8515a21cbef71fcaa303e4050"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}