diff --git a/micom_2023.ipynb b/micom_2023.ipynb
index 9ef1d5e..eeebb11 100644
--- a/micom_2023.ipynb
+++ b/micom_2023.ipynb
@@ -1,843 +1,881 @@
{
- "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 2023 ISB Microbiome Course. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2023/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-3. 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-6.\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_2023 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": "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 π¦Ίπ 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",
- "\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/dada/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 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": "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",
- "collapsed": true
- },
- "outputs": [],
- "source": [
- "import pandas as pd\n",
- "\n",
- "metadata = pd.read_table(\"data/metadata.tsv\").rename(columns={\"sample-id\": \"sample_id\"})\n",
- "tax = pd.merge(tax, metadata, on=\"sample_id\")\n",
- "tax"
- ]
- },
- {
- "cell_type": "markdown",
- "source": [
- "Ok, now we want to invade our samples with C. diff. The goal is to predict susceptibility to invasion and see how disease context can influence predicted engraftment. To do this we will introduce 10% C. diff to all the samples! You can read more about this approach and its applications [here](https://www.biorxiv.org/content/10.1101/2023.04.28.538771v2) "
- ],
- "metadata": {
- "id": "ZsFgo_CRjvsr"
- }
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "id": "vHyvSZJJ-iIi"
- },
- "outputs": [],
- "source": [
- "def invade(data,invader,invader_rel): #define a function that will handle community invasion task\n",
- " invaded=pd.DataFrame() # set up results dataframe\n",
- " for smp,df in data.groupby(by='sample_id'): # loop through data, one sample at a time\n",
- " df=df[df.relative>0.05].copy() # filter out genera below 5% (this is a high threshold, which will make simulation run faster)\n",
- " df=df[df.genus!='Clostridioides'] # some samples already have C. diff, so lets remove it and then re-introduce\n",
- " abund=df.abundance.sum()*invader_rel/(1-invader_rel) # calculate the abundance needed to achieve desired relative abundance\n",
- "\n",
- " info=df.iloc[0,:].copy() # get necessary sample info\n",
- " info.genus=invader # add invader name\n",
- " info.id=invader\n",
- " info.abundance=abund # add invader abundance\n",
- " df.loc[df.shape[0]+1,info.index]=info.values #append invader info\n",
- " df.relative=df.abundance.apply(lambda x:x/df.abundance.sum()) # re-calculate relative abundance\n",
- " invaded=pd.concat([invaded,df]) # append results to output dataframe\n",
- " return invaded\n",
- "\n",
- "invaded=invade(tax,'Clostridioides',0.1) #lets test the function using Clostridioides introducd at 10%\n",
- "invaded"
- ]
- },
- {
- "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 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"
- },
- "outputs": [],
- "source": [
- "from micom.workflows import build\n",
- "from micom import Community\n",
- "import pandas as pd\n",
- "\n",
- "manifest = build(invaded, \"agora103_genus.qza\", \"models\", solver=\"osqp\",\n",
- " cutoff=0.05, threads=2)"
- ]
- },
- {
- "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",
- "### 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 the invasion potential of C. diff. How does its ability to invade vary in samples from healthy donors versus individuals with dysbiotic gut microbiomes (pre-FMT)?\n",
- "\n",
- "Additionally, we can use MICOM to take a mechanistic look at what metabolic strategies are leveraged by C. diff (e.g., what niche(s) does it occupy) in these different contexts.\n",
- "\n",
- "All that and more, coming up. Stay tuned!"
- ]
- },
- {
- "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"
- },
- "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": "code",
- "execution_count": null,
- "metadata": {
- "id": "lQW2BBS10jdN"
- },
- "outputs": [],
- "source": [
- "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 `\"tre[e]\"` represents."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "id": "DphXa9hw1yxM"
- },
- "outputs": [],
- "source": [
- "anns = growth.annotations\n",
- "anns[anns.metabolite == \"tre[e]\"]"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "id": "GVHLD2dm4a6B"
- },
- "source": [
- "Trehalose! Interesting, [that's an important metabolite](https://pubmed.ncbi.nlm.nih.gov/34277467/) in the context of CDI! All of these annotations and more information at are also available at https://vmh.life, maintained by Dr. Ines Thiele's lab."
- ]
- },
- {
- "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",
- "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."
- ],
- "metadata": {
- "id": "_sjNIuCSXkyb"
- }
- },
- {
- "cell_type": "markdown",
- "source": [
- "We're interested in understanding the invasion potential of C. diff so lets extract the predicted C. diff growth rates. In addition to C. diff growth rate we can also look at what fraction of the community growth rate this represents.\n",
- "\n"
- ],
- "metadata": {
- "id": "qPTCyw-7RqFm"
- }
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "id": "9YlWytSp-yKI"
- },
- "outputs": [],
- "source": [
- "cdiff = growth_rates[growth_rates.taxon=='Clostridioides'].copy()\n",
- "cdiff"
- ]
- },
- {
- "cell_type": "markdown",
- "source": [
- "Now that we've extracted the C. diff specific growth rates lets take a look at how they compare between patients with healthy and disbiotic gut microbiomes"
- ],
- "metadata": {
- "id": "bQKYxSZTSPk1"
- }
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "id": "AWG3U094-_OI"
- },
- "outputs": [],
- "source": [
- "import seaborn as sns\n",
- "sns.boxplot(x='disease_state',y='growth_rate',data=cdiff)"
- ]
- },
- {
- "cell_type": "markdown",
- "source": [
- "Looks like C. diff is predicted to grow in all samples but its predicted growth rate is ~2x higher in the Pre-FMT samples. You can see there is also a decent amount of variation in the Pre-FMT results."
- ],
- "metadata": {
- "id": "MEuym_PVui-k"
- }
- },
- {
- "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",
- "source": [
- "One small take away from this analysis is the speration between C. diff pre- and post- FMT samples, suggesting that C. diff may leverage different metabolic strategies in these contexts. Lets take a closer look at this..."
- ],
- "metadata": {
- "id": "weoSaLDPZYAK"
- }
- },
- {
- "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 consumption profile of the microbiome changes when the disease state changes, as changes in microbiome context can lead to changes in host succeptibility to infection. To look into this deeper, we'll transform the microbiome import data and then plot the metabolite imports on a heatmap."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "id": "pAtec3I78DJJ"
- },
- "source": [
- "We can use the `consumption_rates` function in MICOM to calculate consumption rates from the growth results. This will tell us what the patient microbiomes are consuming and provide additional insight into available niches. 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. Importantly, here we consider the consumption rates for samples with no C. diff present to understand how the initial state of the patient microbiomes may influence invasion potential."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "id": "DJAsxjnFdajN"
- },
- "outputs": [],
- "source": [
- "from micom.workflows import load_results\n",
- "from micom.measures import consumption_rates\n",
- "import numpy as np\n",
- "no_cdiff_growth = load_results(\"treasure_chest/no-cdiff-growth.zip\") # Load growth results with no C. diff invasion\n",
- "exchanges = consumption_rates(no_cdiff_growth) # extract consumption 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 = ['disease_state','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 disease context is important - there are apparent differences in consumption rates between the healthy and pre-fmt microbiomes. These differences may be why C. diff can exploit the pre-FMT microbiomes and achieve higher predicted growth rates."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "id": "l8maOr3w2bOo"
- },
- "source": [
- "# π« Exercises"
- ]
- },
- {
- "cell_type": "markdown",
- "source": [
- "Time for you to try your hand at some analysis, lets take a closer looks at the metabolic strategies used by C. diff"
- ],
- "metadata": {
- "id": "XZt3ojssnCzz"
- }
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "id": "VjUustQG2_bX"
- },
- "source": [
- "## Metabolic strategies used by C. diff\n",
- "We've alread looked at the community wide consumption fluxes in the absencence of C. diff and found that they differ between disease contexts. What about the import fluxes of C. diff specifically? Can you develop a visualization to look at those?"
- ]
- },
- {
- "cell_type": "code",
- "source": [
- "exchanges=#your code here # extract exchanges from growth data\n",
- "cdiff_exchanges=#your code here # get cdiff specific exchanges\n",
- "cdiff_imports=#your code here # specifically look at imports\n",
- "cdiff_imports=#your code here # add metadata\n",
- "cdiff_imports = pd.pivot_table(#your code here) # convert to a matrix of samples vs. metabolites\n",
- " # that contains the production rates\n",
- "cdiff_imports = abs(cdiff_imports.T.fillna(0.0)) # fill nans with 0s\n",
- "\n",
- "annot=#your code here # optionally map reactions to metabolite names\n",
- "annot.index=annot.reaction\n",
- "cdiff_imports.index=cdiff_imports.index.map(annot.name.to_dict()) # not necessary but makes results more human readible\n",
- "\n",
- "cdiff_imports = #your coder here # ...and a CLR transform again, normalizes the fluxes\n",
- "\n",
- "cdiff_imports = #your code here # sort by variance, highest variance fluxes first\n",
- "\n",
- "cdiff_imports"
- ],
- "metadata": {
- "id": "zrUZSSp62uUt"
- },
- "execution_count": null,
- "outputs": []
- },
- {
- "cell_type": "code",
- "source": [
- "#Make a heatmap with top 50 fluxes with highest variance\n",
- "\n",
- "#your code here"
- ],
- "metadata": {
- "id": "1gB6X3DJanii"
- },
- "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",
- "source": [
- "manifest"
- ],
- "metadata": {
- "id": "I8NnkbnymxoE"
- },
- "execution_count": null,
- "outputs": []
- },
- {
- "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": "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",
- "[](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.9.7"
- },
- "vscode": {
- "interpreter": {
- "hash": "c991a7ed881363492957ff225bb30af9d5174cd8515a21cbef71fcaa303e4050"
- }
+ "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 2023 ISB Microbiome Course. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2023/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-3. 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-6.\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_2023 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": "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 π¦Ίπ 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",
+ "\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/dada/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 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": "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": {
+ "collapsed": true,
+ "id": "X9hqoO4go0h1",
+ "jupyter": {
+ "outputs_hidden": true
}
+ },
+ "outputs": [],
+ "source": [
+ "import pandas as pd\n",
+ "\n",
+ "metadata = pd.read_table(\"data/metadata.tsv\").rename(columns={\"sample-id\": \"sample_id\"})\n",
+ "tax = pd.merge(tax, metadata, on=\"sample_id\")\n",
+ "tax"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "ZsFgo_CRjvsr"
+ },
+ "source": [
+ "Ok, now we want to invade our samples with C. diff. The goal is to predict susceptibility to invasion and see how disease context can influence predicted engraftment. To do this we will introduce 10% C. diff to all the samples! You can read more about this approach and its applications [here](https://www.biorxiv.org/content/10.1101/2023.04.28.538771v2) "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "vHyvSZJJ-iIi"
+ },
+ "outputs": [],
+ "source": [
+ "def invade(data, invader, invader_rel):\n",
+ " \"\"\"Artificially invade an existing abundance table with another genus.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " data : pandas.DataFrame\n",
+ " A MICOM-style taxonomy table.\n",
+ " invader : str\n",
+ " The genus of the invading taxon.\n",
+ " invader_rel : float in (0, 1)\n",
+ " The relative abundance the invading taxon should have in\n",
+ " each sample.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " pandas.DataFrame\n",
+ " The table with the invaded samples.\n",
+ "\n",
+ " \"\"\"\n",
+ " invaded = pd.DataFrame() # set up results dataframe\n",
+ " for smp, df in data.groupby(by='sample_id'): # loop through data, one sample at a time\n",
+ " df = df[df.relative > 0.05].copy() # filter out genera below 5% (this is a high threshold, which will make simulation run faster)\n",
+ " df = df[df.genus != 'Clostridioides'] # some samples already have C. diff, so lets remove it and then re-introduce\n",
+ " abund = df.abundance.sum() * invader_rel / (1-invader_rel) # calculate the abundance needed to achieve desired relative abundance\n",
+ "\n",
+ " info = df.iloc[0, :].copy() # get necessary sample info\n",
+ " info.genus = invader # add invader name\n",
+ " info.id = invader\n",
+ " info.abundance = abund # add invader abundance\n",
+ " df.loc[df.shape[0] + 1, info.index] = info.values # append invader info\n",
+ " df.relative = df.abundance.apply( # re-calculate relative abundance\n",
+ " lambda x:x/df.abundance.sum()\n",
+ " )\n",
+ " invaded = pd.concat([invaded, df]) # append results to output dataframe\n",
+ " return invaded"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's test the function using *Clostridioides* introducd at 10% final relative abundance."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "invaded = invade(tax, 'Clostridioides', 0.1)\n",
+ "invaded[invaded.genus == 'Clostridioides']"
+ ]
+ },
+ {
+ "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 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"
+ },
+ "outputs": [],
+ "source": [
+ "from micom.workflows import build\n",
+ "from micom import Community\n",
+ "import pandas as pd\n",
+ "\n",
+ "manifest = build(invaded, \"agora103_genus.qza\", \"models\", solver=\"osqp\",\n",
+ " cutoff=0.05, threads=2)"
+ ]
+ },
+ {
+ "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",
+ "### 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 the invasion potential of C. diff. How does its ability to invade vary in samples from healthy donors versus individuals with dysbiotic gut microbiomes (pre-FMT)?\n",
+ "\n",
+ "Additionally, we can use MICOM to take a mechanistic look at what metabolic strategies are leveraged by C. diff (e.g., what niche(s) does it occupy) in these different contexts.\n",
+ "\n",
+ "All that and more, coming up. Stay tuned!"
+ ]
+ },
+ {
+ "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"
+ },
+ "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": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "lQW2BBS10jdN"
+ },
+ "outputs": [],
+ "source": [
+ "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 `\"tre[e]\"` represents."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "DphXa9hw1yxM"
+ },
+ "outputs": [],
+ "source": [
+ "anns = growth.annotations\n",
+ "anns[anns.metabolite == \"tre[e]\"]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "GVHLD2dm4a6B"
+ },
+ "source": [
+ "Trehalose! Interesting, [that's an important metabolite](https://pubmed.ncbi.nlm.nih.gov/34277467/) in the context of CDI! All of these annotations and more information at 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()` 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."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "qPTCyw-7RqFm"
+ },
+ "source": [
+ "We're interested in understanding the invasion potential of C. diff so lets extract the predicted C. diff growth rates. In addition to C. diff growth rate we can also look at what fraction of the community growth rate this represents.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "9YlWytSp-yKI"
+ },
+ "outputs": [],
+ "source": [
+ "cdiff = growth_rates[growth_rates.taxon=='Clostridioides'].copy()\n",
+ "cdiff"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "bQKYxSZTSPk1"
+ },
+ "source": [
+ "Now that we've extracted the C. diff specific growth rates lets take a look at how they compare between patients with healthy and disbiotic gut microbiomes"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "AWG3U094-_OI"
+ },
+ "outputs": [],
+ "source": [
+ "import seaborn as sns\n",
+ "sns.boxplot(x='disease_state',y='growth_rate',data=cdiff)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "MEuym_PVui-k"
+ },
+ "source": [
+ "Looks like C. diff is predicted to grow in all samples but its predicted growth rate is ~2x higher in the Pre-FMT samples. You can see there is also a decent amount of variation in the Pre-FMT results."
+ ]
+ },
+ {
+ "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\")\n"
+ ]
+ },
+ {
+ "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": [
+ "One small take away from this analysis is the speration between C. diff pre- and post- FMT samples, suggesting that C. diff may leverage different metabolic strategies in these contexts. Lets take a closer look at this..."
+ ]
+ },
+ {
+ "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 consumption profile of the microbiome changes when the disease state changes, as changes in microbiome context can lead to changes in host succeptibility to infection. To look into this deeper, we'll transform the microbiome import data and then plot the metabolite imports on a heatmap."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "pAtec3I78DJJ"
+ },
+ "source": [
+ "We can use the `consumption_rates` function in MICOM to calculate consumption rates from the growth results. This will tell us what the patient microbiomes are consuming and provide additional insight into available niches. 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. Importantly, here we consider the consumption rates for samples with no C. diff present to understand how the initial state of the patient microbiomes may influence invasion potential."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "DJAsxjnFdajN"
+ },
+ "outputs": [],
+ "source": [
+ "from micom.workflows import load_results\n",
+ "from micom.measures import consumption_rates\n",
+ "import numpy as np\n",
+ "\n",
+ "no_cdiff_growth = load_results(\"treasure_chest/no-cdiff-growth.zip\") # Load growth results with no C. diff invasion\n",
+ "exchanges = consumption_rates(no_cdiff_growth) # extract consumption 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 = ['disease_state', '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",
+ "\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 disease context is important - there are apparent differences in consumption rates between the healthy and pre-fmt microbiomes. These differences may be why C. diff can exploit the pre-FMT microbiomes and achieve higher predicted growth rates."
+ ]
+ },
+ {
+ "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, lets take a closer looks at the metabolic strategies used by C. diff"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "id": "VjUustQG2_bX"
+ },
+ "source": [
+ "## Metabolic strategies used by C. diff\n",
+ "We've alread looked at the community wide consumption fluxes in the absencence of C. diff and found that they differ between disease contexts. What about the import fluxes of C. diff specifically? Can you develop a visualization to look at those?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "id": "zrUZSSp62uUt"
+ },
+ "outputs": [],
+ "source": [
+ "exchanges = #your code here # extract exchanges from growth data\n",
+ "cdiff_exchanges = #your code here # get cdiff specific exchanges\n",
+ "cdiff_imports = #your code here # specifically look at imports\n",
+ "cdiff_imports = #your code here # add metadata\n",
+ "cdiff_imports = pd.pivot_table(#your code here) # convert to a matrix of samples vs. metabolites\n",
+ " # that contains the production rates\n",
+ "cdiff_imports = abs(cdiff_imports.T.fillna(0.0)) # fill nans with 0s\n",
+ "\n",
+ "annot = #your code here # optionally map reactions to metabolite names\n",
+ "annot.index = annot.reaction\n",
+ "cdiff_imports.index = cdiff_imports.index.map(annot.name.to_dict()) # not necessary but makes results more human readible\n",
+ "\n",
+ "cdiff_imports = #your coder here # and a log-transform\n",
+ "\n",
+ "cdiff_imports = #your code here # sort by variance, highest variance fluxes first\n",
+ "\n",
+ "cdiff_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": "I8NnkbnymxoE"
+ },
+ "outputs": [],
+ "source": [
+ "manifest"
+ ]
+ },
+ {
+ "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": "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",
+ "[](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.11.6"
},
- "nbformat": 4,
- "nbformat_minor": 0
+ "vscode": {
+ "interpreter": {
+ "hash": "c991a7ed881363492957ff225bb30af9d5174cd8515a21cbef71fcaa303e4050"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
}