micom_2024.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "C0vqP4LJ9y6K"
   },
   "source": [
    "# 🧫🦠 Modeling microbiota-wide metabolism with MICOM\n",
    "\n",
    "This notebook will accompany the second session of the 2024 ISB Microbiome Course. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2024/micom).\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 certify 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 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_2024 materials\n",
    "%cd materials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "is6fmBUeorwv"
   },
   "source": [
    "## Basic Installation\n",
    "\n",
    "Installing MICOM is straightforward in Python. OSQP itself will be installed automatically along with it."
   ]
  },
  {
   "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": "M_TeC5yrst3h"
   },
   "outputs": [],
   "source": [
    "!pip install -q micom Cython 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 🦺🛠😁.\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/dada/table.qza\",\n",
    "    \"treasure_chest/taxa.qza\",\n",
    "    collapse_on=\"genus\"\n",
    ")\n",
    "tax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-TQ6Zp7wouk4"
   },
   "source": [
    "Notice the `collapse_on` argument. That will specify the rank on which to summarize 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": "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 sample id, disease state, and demographic information of each of the study participants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "X9hqoO4go0h1"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "metadata = pd.read_csv(\"data/metadata.tsv\",sep='\\t').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, it's 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 passed our taxonomy table and model database to MICOM. We will remove all taxa that make up less than 5% 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",
    "scrolled": true
   },
   "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.05, threads=2)"
   ]
  },
  {
   "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",
    "### Microbiome 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 want to investigate how the activity (e.g., how the metabolic production and consumption patterns) of microbiomes from patients with parkinson's disease might differ from those of healthy controls. Our goal here is to assess whether the microbiome might be playing a role in disease.\n",
    "\n",
    "Alright, lets get to it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "etrvjwBLkKdR"
   },
   "source": [
    "First we need to import our dietary context. For simplicity we will be using a formulation that represents an average western diet, but if information about host diet is known other formulations can be used (e.g., vegetarian or vegan diet). Additional dietary formulations can be found [here]( https://github.com/micom-dev/media/tree/main/media)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aJwXAR4PPAkA"
   },
   "outputs": [],
   "source": [
    "from micom.qiime_formats import load_qiime_medium\n",
    "medium = load_qiime_medium(\"western_diet_gut_agora.qza\")\n",
    "medium"
   ]
  },
  {
   "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. 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. This will take the models we've built, and find an optimal solution to the fluxes based upon the medium that's been applied. Should take about 10 minutes to complete."
   ]
  },
  {
   "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",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from micom.workflows import grow, save_results\n",
    "\n",
    "growth = grow(manifest, \"models\",medium, tradeoff=0.8, threads=2)\n",
    "\n",
    "# We'll save the results to a file\n",
    "save_results(growth, \"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",
    "  growth\n",
    "except NameError:\n",
    "  growth = load_results(\"treasure_chest/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": "2r_XUm7U-HSm"
   },
   "outputs": [],
   "source": [
    "growth_rates=pd.merge(growth.growth_rates,metadata,on='sample_id')\n",
    "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": "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 from 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 = growth.annotations\n",
    "anns[anns.metabolite == \"but[e]\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GVHLD2dm4a6B"
   },
   "source": [
    "Butyrate! Interesting, [that's an important metabolite](https://www.nature.com/articles/s41392-024-01743-1) in the context of the gut brain axis! All of these annotations and more information are also available at https://vmh.life, maintained by Drs. Ronan Fleming's and Ines Thiele's labs."
   ]
  },
  {
   "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 = plot_growth(growth, filename=\"growthrates.html\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_sjNIuCSXkyb"
   },
   "source": [
    "Normally, we could call `viz.view()` afterward 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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qPTCyw-7RqFm"
   },
   "source": [
    "We're interested in understanding the role of microbiome activity in parkinson's disease. Prior work has suggested short-chain fatty acids (SCFAs) and genera that produce them, like  Faecalibacterium, Roseburia, and Blautia may play a role in mitigating the onset and severity of disease. Let's take a look at the growth rates of these genera and how they vary across cases and controls.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9YlWytSp-yKI"
   },
   "outputs": [],
   "source": [
    "fcali = growth_rates[growth_rates.taxon=='g__Faecalibacterium'].copy()\n",
    "fcali"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bQKYxSZTSPk1"
   },
   "source": [
    "Now that we've extracted the Faecalibacterium specific growth rates let's take a look at how they compare between patients with and without parkinson's disease."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "AWG3U094-_OI"
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "sns.boxplot(x='parkinson_disease',y='growth_rate',data=fcali)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also at other genera to see if there are any significant differences in predicted growth rates between those with and without parkinson's disease"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "genera=growth_rates.groupby(by=['taxon']).count()['abundance'].nlargest(10).index\n",
    "sns.boxplot(x='growth_rate',y='taxon',hue='parkinson_disease',data=growth_rates[growth_rates.taxon.isin(genera)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like there are indeed some differences that are more substantial among some of the most prevalent genera. However, the differences we see are possibly not what we would have expected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "G1JbbKrLcVye"
   },
   "source": [
    "## Growth niches\n",
    "\n",
    "Another thing we can look at is whether individual taxa inhabit different growth niches across different disease 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(growth, perplexity=4, direction=\"import\", filename=\"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": "weoSaLDPZYAK"
   },
   "source": [
    "We see that most genera don't show much separation among our samples suggesting that the growth niches occupied across our sample may not be very different. However, this doesn't seem to be true for all genera! We do see some separation for Bacteroides!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4Y_XfHkB4sO8"
   },
   "source": [
    "## Comparative flux analysis\n",
    "\n",
    "Now let's compare the metabolite imports between the two disease contexts. We're interested to see how the production profile of the microbiome changes when the disease state changes. To look into this deeper, we'll transform the microbiome import 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. This will tell us what the patient microbiomes are producing and provide additional insight into the role microbial metabolites might play in disease. To visualize the results we'll  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": "markdown",
   "metadata": {},
   "source": [
    "First, let's take a quick at the SFCA butyrate and which genera produce it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "exports=growth.exchanges #Get samples exchanges \n",
    "exports=exports[exports.direction=='export'] #Filter for  exports\n",
    "but=exports[exports.reaction=='EX_but(e)'] #Get butyrate specific exports\n",
    "but['logflux']=np.log(but['flux']+0.001) #Log transform flux and account for possibility of zeros\n",
    "#Lets looks at top 3 genera with most non-zero butyrate fluxes\n",
    "genera=but[but.logflux>-6].groupby(by='taxon').count()['sample_id'].sort_values(ascending=False).index[:3]\n",
    "plt.figure(figsize=(5,3)) #Optional function call to specify figure size\n",
    "sns.boxplot(y='taxon',x='logflux',data=but,order=genera)\n",
    "plt.title('Butyrate') #Add title"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets see how the genera specific fluxes vary with disease context!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "but_meta=pd.merge(but,metadata,on='sample_id')\n",
    "plt.figure(figsize=(5,3))\n",
    "sns.boxplot(y='taxon',x='logflux',data=but_meta,hue='parkinson_disease',order=genera)\n",
    "plt.title('Butyrate')\n",
    "plt.legend(bbox_to_anchor=(1,1),title=\"Parkison's disease\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets try visualizing this result as a heatmap so we can get a better idea of the variability in individual samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logbut=but_meta.pivot_table(columns='taxon',index=['parkinson_disease','sample_id'],values='flux').fillna(0).apply(                           # and we log-transform the fluxes\n",
    "    lambda xs: np.log(xs + 0.001),\n",
    "    axis=0) #Convert data table from long to wide format and log transform\n",
    "\n",
    "#Create palette to color heatmap columns based on disease status\n",
    "df=logbut.reset_index().copy()\n",
    "i=0\n",
    "cmap={}\n",
    "colors=sns.color_palette('colorblind')\n",
    "for group in sorted(df['parkinson_disease'].unique()):\n",
    "    cmap[group] = colors[i]\n",
    "    i+=3\n",
    "#Use color dictionary to map color to speces column (provides color of each data point)\n",
    "col_colors=df['parkinson_disease'].copy()\n",
    "col_colors=col_colors.map(cmap)\n",
    "col_colors.index=logbut.index\n",
    "\n",
    "#Create clustered heatmap \n",
    "fig=sns.clustermap(logbut[genera].T,figsize=(5,2.5),cmap='viridis',xticklabels='',col_colors=col_colors)\n",
    "fig.ax_heatmap.set_xlabel('')\n",
    "fig.ax_heatmap.set_ylabel('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this visualization we can see how butyrate flux varies with both genera and disease context\n",
    "\n",
    "Now lets take a look at the global metabolite fluxes!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again we can see how burytate production varies across disease contexts!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from micom.measures import consumption_rates, production_rates\n",
    "\n",
    "#Extract production rates\n",
    "exchanges = consumption_rates(growth)    \n",
    "#Add metadata\n",
    "exchanges=pd.merge(exchanges,metadata,on='sample_id') \n",
    "#Log transform\n",
    "exchanges['logflux']=np.log(exchanges.flux+0.001)\n",
    "#Visualize\n",
    "sns.boxplot(x='parkinson_disease',y='logflux',data=exchanges[exchanges.metabolite=='but[e]'])\n",
    "plt.title('Butyrate')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like globally, butyrate flux is predicted to be a bit higher on average for those without parkinson's disease"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Now lets see how things look across all exports!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exchanges = production_rates(growth)                   # extract production rates\n",
    "exchanges=pd.merge(exchanges,metadata,on='sample_id')  # add metadata\n",
    "exchanges = pd.pivot_table(                            # convert to a matrix of samples vs. metabolites\n",
    "    exchanges,                                         # that contains the production rates\n",
    "    index = ['parkinson_disease', '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 we log-transform the fluxes\n",
    "    lambda xs: np.log(xs + 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",
    "exchanges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Exports\n",
    "\n",
    "#Create palette to color heatmap columns based on disease status\n",
    "df=exchanges.T.reset_index().copy()\n",
    "i=0\n",
    "cmap={}\n",
    "colors=sns.color_palette('colorblind')\n",
    "for group in sorted(df['parkinson_disease'].unique()):\n",
    "    cmap[group] = colors[i]\n",
    "    i+=3\n",
    "#Use color dictionary to map color to species column (provides the color of each data point)\n",
    "col_colors=df['parkinson_disease'].copy()\n",
    "col_colors=col_colors.map(cmap)\n",
    "col_colors.index=exchanges.columns\n",
    "\n",
    "fig=sns.clustermap(\n",
    "    exchanges.head(25),  # take 50 highest fluxes\n",
    "    cmap='viridis',\n",
    "    yticklabels=True,    # show all metabolite names\n",
    "    figsize=(12, 10),\n",
    "    col_colors=col_colors,# size of the heatmap\n",
    "    xticklabels=''\n",
    ")\n",
    "fig.ax_heatmap.set_xlabel('')\n",
    "fig.ax_heatmap.set_ylabel('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Globally, it appears there is some indication that the metabolic activity of the microbiome varies with disease context, however, given the limited sample size and the approximations that had to be made to generate these results it is hard to say with certainty if the differences in predictions we see are plausible or significant "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "l8maOr3w2bOo"
   },
   "source": [
    "# 🏫 Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XZt3ojssnCzz"
   },
   "source": [
    "Time for you to try your hand at some analysis, let's take a closer look at the metabolic strategies used by genera in the community"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VjUustQG2_bX"
   },
   "source": [
    "## Metabolic strategies used by a genus of your choice!\n",
    "We've already looked at the community wide production fluxes and found that they differ slightly between disease contexts. What about the import fluxes of the most prevalent genera? Can you develop a visualization to look at those?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's ensure you don't pick a genus that is only present in 1 sample. We want to see how imports vary across disease contexts. Use the code below to pick a genera that is present in at least 7 samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Count how many samples each genera is present in\n",
    "growth.growth_rates.groupby(by='taxon').count()['abundance'].nlargest(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "zrUZSSp62uUt"
   },
   "outputs": [],
   "source": [
    "exchanges = #your code here                                          # extract exchanges from growth data\n",
    "genus_exchanges = #your code here                                    # get genus specific exchanges\n",
    "genus_imports = #your code here                                      # specifically look at imports\n",
    "genus_imports = #your code here                                      # add metadata\n",
    "genus_imports = #your code here                                      # convert to a matrix of samples vs. metabolites\n",
    "                                                                     # that contains the production rates\n",
    "genus_imports = abs(genus_imports.T.fillna(0.0))                     # fill nans with 0s\n",
    "\n",
    "annot = #your code here                                              # optionally map reactions to metabolite names\n",
    "annot.index = #your code here\n",
    "genus_imports.index = genus_imports.index.map(annot.name.to_dict())  # not necessary but makes results more human readible\n",
    "\n",
    "genus_imports = #your code here                                      # and a log-transform\n",
    "\n",
    "genus_imports=#your code here]                                       # sort by variance, highest variance fluxes first\n",
    "\n",
    "genus_imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1gB6X3DJanii"
   },
   "outputs": [],
   "source": [
    "#Make a heatmap with top 50 fluxes with highest variance\n",
    "\n",
    "#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\", medium, threads=2)\n",
    "tradeoff_results.to_csv(\"tradeoff.csv\", index=False)\n",
    "\n",
    "plot_tradeoff(tradeoff_results, tolerance=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from micom.viz.tradeoff import plot_tradeoff\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": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.19"
  },
  "vscode": {
   "interpreter": {
    "hash": "c991a7ed881363492957ff225bb30af9d5174cd8515a21cbef71fcaa303e4050"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}