micom.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "micom.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/isb_course_2020/blob/master/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 2020 ISB Microbiome Course. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2020/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/). 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. However, MICOM will require a solver for quadratic programming problems and all the best ones are commercial (boo) even though they often have free academic licenses 😌. We will use a brand new open source QP solver named OSQP, but this will require us to pull in development versions of certain packages.\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/isb_course_2020 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. OSQP itself can be installed right along with it."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M_TeC5yrst3h"
      },
      "source": [
        "!pip install -q osqp cobra micom\n",
        "\n",
        "print(\"Done! 🎉 \")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l9ssgi8ts_fP"
      },
      "source": [
        "## Enabling OSQP support\n",
        "\n",
        "For today we will also install some development versions that enable full support for OSQP in MICOM. This will hopefully not be necessary anymore soon. Also this would be unnecessary if we would have access to CPLEX og Gurobi (both commercial solvers with free academic licenses)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p5zNnjaCotQ_"
      },
      "source": [
        "!rm -r /tmp /var/tmp\n",
        "!pip install --force-reinstall --no-deps -q \\\n",
        "  https://github.com/cdiener/cobrapy/archive/feature/osqp_support.zip \\\n",
        "  https://github.com/cdiener/optlang/archive/feature/osqp.zip\n",
        "  \n",
        "print(\"Done! 🎉 \")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6oJrxxz6tV9T"
      },
      "source": [
        "## Enable QIIME 2 interactions\n",
        "\n",
        "Finally we will need to install packages to read the \"biom\" format which is a file format QIIME 2 uses to save tables. This is only necessary if you want to read QIIME 2 FeatureTable artifacts."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5rZX7SK_toLp"
      },
      "source": [
        "!pip install -q numpy Cython\n",
        "!pip install -q biom-format\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. To see how those tables look we can import MICOM and look at an example table:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UV9SObSQkSZh"
      },
      "source": [
        "from micom.data import test_data\n",
        "\n",
        "test_data().head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AEk7yfd1lbYp"
      },
      "source": [
        "The `file` column is not required when using a taxonomy database like we will do here. The `id` column specifies identifiers for the taxa and should be expressive and not include spaces or special characters. Each row needs to contain the the abundance of a single taxon in a single sample. \n",
        "\n",
        "Oh no, that's not what we have generated in the previous step. We only have 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, 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 [MICOM data repository](https://doi.org/10.5281/zenodo.3755182) which contains a whole lot of prebuilt databases.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "eXGHZPYwlEa5"
      },
      "source": [
        "!wget -O agora103_genus.qza https://zenodo.org/record/3755182/files/agora103_genus.qza?download=1"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_57iya0D3L6-"
      },
      "source": [
        "Okay. We've got everything we need now. The data from the prior analysis can be found in the `treasure_chest` folder, so we can use those files. We will remove all taxa that make up less than 2.5% of the community to keep the models small and speed up this tutorial."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "o0vBAiiqqPLC"
      },
      "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",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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 are you to run into those 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. The GreenGenes database is pretty old and many taxonomic names have been superceded by now. So we will stick with the \"lax\" option of only matching on genus ranks.\n",
        "\n",
        "We can also look at the generated MICOM taxonomy."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2eS_g-ffNUAt"
      },
      "source": [
        "tax"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R5zwLVFcNTbq"
      },
      "source": [
        "One helpful thing to do is to merge in our metadata, so we'll have it at hand for the following steps."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "X9hqoO4go0h1"
      },
      "source": [
        "import pandas as pd\n",
        "\n",
        "metadata = pd.read_table(\"metadata.tsv\").rename(columns={\"id\": \"sample_id\"})\n",
        "tax = pd.merge(tax, metadata, on=\"sample_id\")\n",
        "tax"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RTbYBR8cJfup"
      },
      "source": [
        "With the taxonomic metadata, we can finally build our community-level models.\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 have to specify where to write the models. We will also run that in parallel over two threads. It should take around 10 minutes to finish."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XlCk7GyB7B5e"
      },
      "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=2.5e-2, threads=2)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kwya6vbZZSUo"
      },
      "source": [
        "For different data a warning may pop up if less than 50% of the abundances can be matched to the database. If this happens, you can still continue, but be aware that such a sparse model may not accurately represent your sample. 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",
        "We won't see any warnings here. So, we will go ahead for now. Let's also take a look what we got back from the `build` process."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r9qwglr88Ise"
      },
      "source": [
        "manifest"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y4KAJkhIdspQ"
      },
      "source": [
        "This will tell you many taxa were found in the database and what fraction of the total abundance was represented by the database. Looks okay here. \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://doi.org/10.5281/zenodo.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 again 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 (indicated by the dilution column). \n",
        "\n",
        "Okay let's go right ahead and simulate growth. This will take a little while and give us time to dive into some details 🏊"
      ]
    },
    {
      "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)\n",
        "\n",
        "# We'll save the results to a file\n",
        "pickle.dump(growth_results, open(\"growth.pickle\", \"wb\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "tZiPmdO4S-OE"
      },
      "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.\n",
        "\n",
        "The first things we might want to look at are the growth rates for each taxon."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EaplMHFLcMT7"
      },
      "source": [
        "from micom.viz import *\n",
        "\n",
        "viz = plot_growth(growth_results)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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": "Mkgydl_jdqc7"
      },
      "source": [
        "plot_exchanges_per_sample(growth_results)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lRDCXRDidupd"
      },
      "source": [
        "We can have a look at the results after downloading `materials/sample_exchanges_[DATE].html`. It would be even better if we could visualize which taxa compete for similar resources. We can create a niche plot by using `plot_exchanges_per_taxon`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NlZrfv38esj8"
      },
      "source": [
        "plot_exchanges_per_taxon(growth_results, perplexity=4, direction=\"import\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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 close together. Thus, 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 enough data.\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",
        "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.01 mmol/l."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ebsw0XD-l2Lj"
      },
      "source": [
        "from micom.viz import *\n",
        "\n",
        "manifest.index = manifest.sample_id\n",
        "pheno = manifest.disease_stat\n",
        "\n",
        "pl = plot_fit(growth_results, pheno, atol=1e-3, flux_type=\"production\")\n"
      ],
      "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).\n",
        "\n",
        "Some observations to help you intepret the results:\n",
        "\n",
        "- [maltotriose](https://en.wikipedia.org/wiki/Maltotriose) is a breakdown product of alpha-amylase\n",
        "- [laurate](https://en.wikipedia.org/wiki/Lauric_acid) increases serum cholesterol more than any other fatty acid\n",
        "- *C. diff.* toxins contain multipe aspartate residues and toxin B [has aspartate protease activity](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4885049/) \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": "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",
      "metadata": {
        "id": "bllnCa1O28YZ"
      },
      "source": [
        "manifest"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VjUustQG2_bX"
      },
      "source": [
        "Let's look further into the surprising *C. diff.* individual that looked very similar to the healthy subjects. We will apply the same diet as before."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uh7QFi9V3HIX"
      },
      "source": [
        "from micom import load_pickle\n",
        "from micom.qiime_formats import load_qiime_medium\n",
        "\n",
        "medium = load_qiime_medium(\"western_diet_gut.qza\")\n",
        "medium.index = medium.reaction\n",
        "\n",
        "com = load_pickle(\"models/ERR1883248.pickle\")\n",
        "com.medium = medium.flux\n",
        "com"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GFcWJKuO3bnT"
      },
      "source": [
        "This is a MICOM community object. MICOM community models are full [COBRApy](https://opencobra.github.io/cobrapy/) models (with sprinkles on top) and there is a whole bunch of stuff we could do with it."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CUaUHsKm2fk7"
      },
      "source": [
        "## Microbe-microbe interactions\n",
        "\n",
        "Let's dive a bit more into competition and cooperation between taxa. We can start by simulating taxa knockouts using `com.knockout_taxa`. For that we will remove each of the taxa from the model (one-at-a-time) and see how that affects the growth rates of all other taxa. If other taxa grow faster after the knockout, they were competing with the knocked-out taxon. However, if they grow slower, they were cooperating with the knocked-out taxon. \n",
        "\n",
        "See the docs for more info: https://micom-dev.github.io/micom/taxa_knockouts.html.\n",
        "\n",
        "Bonus points if you can visualize your results. How would you deal with vastly different scales of growth rates between taxa?\n",
        "\n",
        "> Oh geez. I once saw somebody plot a heatmap with Seaborn using a Pandas DataFrame. I wish I could remember where that was... 🤔\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5IlJ-YpUCQzp"
      },
      "source": [
        "# Your code here"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p81Fsdp-CPTF"
      },
      "source": [
        "## The flux capacitator\n",
        "\n",
        "Here is how you would run the cooperative tradeoff for a single model. We can follow that up with pFBA to get all fluxes in the system."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tMAz3FoSMZLf"
      },
      "source": [
        "sol = com.cooperative_tradeoff(fraction=0.5, fluxes=True, pfba=True)\n",
        "sol"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J3j9NZC1cFR2"
      },
      "source": [
        "The returned solution contains all fluxes in `sol.fluxes`. An `NaN` entry denotes that this reaction is not present in the organism."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "X26njphUcJ2k"
      },
      "source": [
        "sol.fluxes.head()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Aeg_X2KQucKK"
      },
      "source": [
        "The acetate export reaction has the name `EX_ac(e)`. Identify the primary acetate producer in the system. Don't forget about the accuracy of 10<sup>-3</sup>."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aWC4kiq10zlO"
      },
      "source": [
        "# Your code here"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LGGBSgwVwZiC"
      },
      "source": [
        "Look up at least one other reaction at https://www.vmh.life/#microbes/reactions, and get it's predicted fluxes. Does the reaction take place? If yes in which organism? Can you identify the most active fluxes in the community?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PFwS8VIbuIY6"
      },
      "source": [
        "# Your code here"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "8_1jesZTHYra"
      },
      "source": [
        "from micom.workflows import tradeoff\n",
        "import micom\n",
        "\n",
        "tradeoff_results = tradeoff(manifest, \"models\", medium, threads=2)\n",
        "tradeoff_results.to_csv(\"tradeoff.csv\", index=False)\n",
        "\n",
        "plot_tradeoff(tradeoff_results, tolerance=1e-4)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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."
      ]
    }
  ]
}