ipynb/tutorials/Phylogenetic-IBD-Tutorial-QIIME2-API.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phylogenetic Compositional Tensor Factorization (phylo-CTF) Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial builds upon the CTF tutorial and explains how you can also perform CTF weighted by phylogenetic information. If you would like to better understand CTF please first complete that tutorial. In this tutorial we will use Phylogenetic Compositional Tensor Factorization (phylo-CTF) which we provide in the software package [gemelli](https://github.com/biocore/gemelli). Phylo-CTF can account for repeated measures, compositionality, and sparsity in microbiome data.\n",
    "\n",
    "In this tutorial we use _gemelli_ to perform phylo-CTF on a time series dataset comparing Crohn's and control subjects over a period of 25 weeks published in [Vázquez-Baeza et al](https://gut.bmj.com/content/67/9/1743). First we will download the processed data originally from [here](https://qiita.ucsd.edu/study/description/2538#). This data can be downloaded with the following links:\n",
    "\n",
    "* **Table** (table.qza) | [download](https://github.com/biocore/gemelli/tree/master/ipynb/tutorials/IBD-2538/data/table.qza)\n",
    "* **Rarefied Table** (rarefied-table.qza) | [download](https://github.com/biocore/gemelli/tree/master/ipynb/tutorials/IBD-2538/data/rarefied-table.qza)\n",
    "* **Sample Metadata** (metadata.tsv) | [download](https://github.com/biocore/gemelli/tree/master/ipynb/tutorials/IBD-2538/data/metadata.tsv)\n",
    "* **Feature Metadata** (taxonomy.qza) | [download](https://github.com/biocore/gemelli/tree/master/ipynb/tutorials/IBD-2538/data/taxonomy.qza)\n",
    "* **Tree** (sepp-insertion-tree.qza) | [download](https://github.com/biocore/gemelli/tree/master/ipynb/tutorials/IBD-2538/data/sepp-insertion-tree.qza)\n",
    "\n",
    "**Note**: This tutorial assumes you have installed [QIIME2](https://qiime2.org/) using one of the procedures in the [install documents](https://docs.qiime2.org/2020.2/install/). This tutorial also assumed you have installed, [Qurro](https://github.com/biocore/qurro) and [gemelli](https://github.com/biocore/gemelli).\n",
    "\n",
    "First, we will make a tutorial directory and download the data above and move the files to the `IBD-2538/data` directory:\n",
    "\n",
    "```bash\n",
    "mkdir IBD-2538\n",
    "```\n",
    "```bash\n",
    "# move downloaded data here\n",
    "mkdir IBD-2538/data\n",
    "```\n",
    "\n",
    "Next we will import our data with the QIIME2 Python API. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import warnings\n",
    "import qiime2 as q2\n",
    "# hide pandas Future/Deprecation Warning(s) for tutorial\n",
    "warnings.filterwarnings(\"ignore\", category=DeprecationWarning) \n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "# import table(s)\n",
    "table = q2.Artifact.load('IBD-2538/data/table.qza')\n",
    "rarefied_table = q2.Artifact.load('IBD-2538/data/rarefied-table.qza')\n",
    "# import metadata\n",
    "metadata = q2.Metadata.load('IBD-2538/data/metadata.tsv')\n",
    "# import tree\n",
    "tree = q2.Artifact.load('IBD-2538/data/sepp-insertion-tree.qza')\n",
    "# import taxonomy\n",
    "taxonomy = q2.Artifact.load('IBD-2538/data/taxonomy.qza')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to account for the correlation among samples from the same subject we will employ phylogenetic compositional tensor factorization (CTF). CTF builds on the ability to account for compositionality and sparsity using the robust center log-ratio transform covered in the RPCA tutorial (found [here](https://forum.qiime2.org/t/robust-aitchison-pca-beta-diversity-with-deicode)) but restructures and factors the data as a tensor. Here we will run CTF through _gemelli_ and explore/interpret the different results. Additionally, phylogenetic CTF incorporates internal nodes and branch lengths of a tree into the factorization. \n",
    "\n",
    "There are two forms of phylogenetic CTF bieng (1) with taxonomy information and (2) without. Here we will run phylogenetic-ctf-with-taxonomy but it can also be run without taxonomy called just phylogenetic-ctf.\n",
    "\n",
    "To run phylo-CTF we only need to run one command (gemelli phylogenetic-ctf-with-taxonomy). The required input requirements are:\n",
    "\n",
    "1. table\n",
    "    - The table is of type `FeatureTable[Frequency]` which is a table where the rows are features (e.g. ASVs/microbes), the columns are samples, and the entries are the number of sequences for each sample-feature pair.\n",
    "2. tree\n",
    "    - This is a phylogeny of type `Phylogeny[Rooted]` where all the features in the `table` are represented in the tree. \n",
    "3. feature-metadata-file\n",
    "    - This is a metadata file (e.g. tsv, or `FeatureTable[Taxonomy]` .qza) where the rows are matched to the table features and the columns are feature metadata such as taxonomy, gene pathway, etc... \n",
    "4. sample-metadata\n",
    "    - This is a QIIME2 formatted [metadata](https://docs.qiime2.org/2020.2/tutorials/metadata/) (e.g. tsv format) where the rows are samples matched to the (1) table and the columns are different sample data (e.g. time point). \n",
    "5. individual-id-column\n",
    "    - This is the name of the column in the (2) metadata that indicates the individual subject/site (e.g. subject ID) that was sampled repeatedly.\n",
    "6. state-column\n",
    "    - This is the name of the column in the (2) metadata that indicates the numeric repeated measure (e.g., Time in months/days) or non-numeric category (i.e. decade/body-site). \n",
    "7. output-dir\n",
    "    - The desired location of the output. We will cover each output independently below.  \n",
    "\n",
    "\n",
    "In this tutorial our individual-id-column is `host_subject_id` and our state-column is different time points denoted as `timepoint` in the sample metadata. Now we are ready to run phylo-CTF:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Results (name = value)\n",
       "--------------------------------------------------------------------------------------------------------------------\n",
       "subject_biplot           = <artifact: PCoAResults % Properties('biplot') uuid: 153a5d84-648a-4ef5-a52c-c8db20a487e8>\n",
       "state_biplot             = <artifact: PCoAResults % Properties('biplot') uuid: 668f9315-9d79-4964-acd5-6f81afb43d5f>\n",
       "distance_matrix          = <artifact: DistanceMatrix uuid: b6ff2ad6-e869-4b3a-aeaa-cf091e983662>\n",
       "state_subject_ordination = <artifact: SampleData[SampleTrajectory] uuid: 5016ae27-ad16-4a30-a806-30b21f9ec115>\n",
       "state_feature_ordination = <artifact: FeatureData[FeatureTrajectory] uuid: dab556b3-06fe-4101-ba9b-ebedbf417795>\n",
       "counts_by_node_tree      = <artifact: Phylogeny[Rooted] uuid: 4ddbace6-2eeb-40b7-b99f-a91924341aa6>\n",
       "counts_by_node           = <artifact: FeatureTable[Frequency] uuid: bc5eb5fa-5313-4f47-84b5-b1cb04c99a4b>\n",
       "t2t_taxonomy             = <artifact: FeatureData[Taxonomy] uuid: 2af2010a-c636-48fe-a62b-5b1c97c82711>\n",
       "subject_table            = <artifact: FeatureTable[Frequency] uuid: 086272d8-3256-4242-bbfb-c04c57bf1fe1>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiime2.plugins.gemelli.actions import phylogenetic_ctf_with_taxonomy\n",
    "from qiime2.plugins.fragment_insertion.methods import filter_features\n",
    "\n",
    "# make a dir. for results\n",
    "#os.mkdir('IBD-2538/ctf-results')\n",
    "# first ensure all the table features are in the tree\n",
    "table = filter_features(table, tree).filtered_table\n",
    "# now run phylo-CTF\n",
    "ctf_results = phylogenetic_ctf_with_taxonomy(table, tree, taxonomy.view(q2.Metadata),\n",
    "                                             metadata, 'host_subject_id', 'timepoint',\n",
    "                                             min_feature_frequency=10)\n",
    "\n",
    "ctf_results\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We will now cover the output files:\n",
    "* subject_biplot\n",
    "* state_biplot\n",
    "* distance_matrix\n",
    "* state_subject_ordination\n",
    "* state_feature_ordination\n",
    "* counts_by_node_tree\n",
    "* counts_by_node\n",
    "* t2t_taxonomy\n",
    "* subject_table\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we will explore the `subject_biplot` which is an ordination where dots represent _subjects_ not _samples_ and arrows represent features (e.g. ASVs). First, we will need to aggregate the metadata by subject (i.e. collapsing the metadata of all samples from a given subject). This can be done by hand or using DataFrames in python (with pandas) or R like so:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ibd</th>\n",
       "      <th>active_disease</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>#SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>s1000100</th>\n",
       "      <td>Control</td>\n",
       "      <td>quiescent</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s1000200</th>\n",
       "      <td>Control</td>\n",
       "      <td>quiescent</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s1000300</th>\n",
       "      <td>Control</td>\n",
       "      <td>quiescent</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s1000500</th>\n",
       "      <td>Control</td>\n",
       "      <td>quiescent</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s1000600</th>\n",
       "      <td>Control</td>\n",
       "      <td>quiescent</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               ibd active_disease\n",
       "#SampleID                        \n",
       "s1000100   Control      quiescent\n",
       "s1000200   Control      quiescent\n",
       "s1000300   Control      quiescent\n",
       "s1000500   Control      quiescent\n",
       "s1000600   Control      quiescent"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from qiime2 import Metadata\n",
    "import numpy as np\n",
    "from biom import Table\n",
    "from skbio import OrdinationResults\n",
    "\n",
    "# first we import the metdata into pandas\n",
    "mf = pd.read_csv('IBD-2538/data/metadata.tsv', sep='\\t',index_col=0)\n",
    "# next we aggregate by subjects (i.e. 'host_subject_id') \n",
    "# and keep the first instance of 'diagnosis_full' by subject.\n",
    "mf = mf.groupby('host_subject_id').agg({'ibd':'first','active_disease':'first'})\n",
    "# now we save the metadata in QIIME2 format.\n",
    "mf.index.name = '#SampleID'\n",
    "mf.to_csv('IBD-2538/data/subject-metadata.tsv', sep='\\t')\n",
    "mf.head(5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Taxon</th>\n",
       "      <th>PC1</th>\n",
       "      <th>PC2</th>\n",
       "      <th>PC3</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>featureid</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGCGCAGGCGGTCTGTTAAGTCAGATGTGAAAGGTTAGGGCTCAACCCTGAACGTG</th>\n",
       "      <td>k__Bacteria; p__Firmicutes; c__Clostridia; o__...</td>\n",
       "      <td>0.042088</td>\n",
       "      <td>-0.033160</td>\n",
       "      <td>-0.039901</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGCGCAGGCGGACCGGCAAGTTGGAAGTGAAAACTATGGGCTCAACCCATAAATTG</th>\n",
       "      <td>k__Bacteria; p__Firmicutes; c__Clostridia; o__...</td>\n",
       "      <td>-0.024657</td>\n",
       "      <td>0.060054</td>\n",
       "      <td>-0.010885</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGCGCAGGCGGACCGGCAAGTTGGAAGTGAAATCCATGGGCTCAACCCGTGAATTG</th>\n",
       "      <td>k__Bacteria; p__Firmicutes; c__Clostridia; o__...</td>\n",
       "      <td>-0.004121</td>\n",
       "      <td>-0.061348</td>\n",
       "      <td>0.006902</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGCGCAGGCGGATTGGCAAGTTGGGAGTGAAATCTATGGGCTCAACCCATAAATTG</th>\n",
       "      <td>k__Bacteria; p__Firmicutes; c__Clostridia; o__...</td>\n",
       "      <td>-0.027500</td>\n",
       "      <td>-0.019843</td>\n",
       "      <td>0.036397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGCGCAGGCGGGAAGACAAGTTGGAAGTGAAAACCATGGGCTCAACCCATGAATTG</th>\n",
       "      <td>k__Bacteria; p__Firmicutes; c__Clostridia; o__...</td>\n",
       "      <td>-0.001026</td>\n",
       "      <td>0.030095</td>\n",
       "      <td>-0.000666</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                                                                Taxon  \\\n",
       "featureid                                                                                               \n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGC...  k__Bacteria; p__Firmicutes; c__Clostridia; o__...   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC...  k__Bacteria; p__Firmicutes; c__Clostridia; o__...   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC...  k__Bacteria; p__Firmicutes; c__Clostridia; o__...   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC...  k__Bacteria; p__Firmicutes; c__Clostridia; o__...   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC...  k__Bacteria; p__Firmicutes; c__Clostridia; o__...   \n",
       "\n",
       "                                                         PC1       PC2  \\\n",
       "featureid                                                                \n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGC...  0.042088 -0.033160   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC... -0.024657  0.060054   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC... -0.004121 -0.061348   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC... -0.027500 -0.019843   \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC... -0.001026  0.030095   \n",
       "\n",
       "                                                         PC3  \n",
       "featureid                                                     \n",
       "AACATAGGGGGCAAGCGTTGTCCGGAATCACTGGGCGTAAAGGGCGC... -0.039901  \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC... -0.010885  \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC...  0.006902  \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC...  0.036397  \n",
       "AACGTAGGGTGCAAGCGTTGTCCGGAATTACTGGGTGTAAAGGGAGC... -0.000666  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# combine feature metadata\n",
    "phylo_ctf_taxonomy = ctf_results.t2t_taxonomy.view(q2.Metadata).to_dataframe()\n",
    "phylo_ctf_feature_loadings = ctf_results.subject_biplot.view(OrdinationResults).features.rename({0:'PC1',\n",
    "                                                                                                 1:\"PC2\",\n",
    "                                                                                                 2:\"PC3\"},\n",
    "                                                                                                axis=1)\n",
    "phylo_ctf_taxonomy_and_loadings = pd.concat([phylo_ctf_taxonomy, phylo_ctf_feature_loadings], axis=1)\n",
    "phylo_ctf_taxonomy_and_loadings.index.name = 'featureid'\n",
    "phylo_ctf_taxonomy_and_loadings.head(5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unlike conventional CTF, the arrows in the phylo-CTF biplot can be both internal nodes and features of the table. In order to label the taxonomy we use the lowest common ancestor of that internal node. To do this we feed the `t2t_taxonomy` output from phylo-CTF into the biplot. Additionally, we can use Empress to generate plots combining both the tree and biplot view at one time. This view can help us understand what phylogenetic partitions separate our samples along a PC axis. \n",
    "\n",
    "Phylo-CTF output contains a table that allows us to use the internal nodes of the tree as possible log-ratios (i.e. `counts_by_node`). In this table the internal nodes contain the sum of all the counts up to that node. `Note,` it is important not to take the log-ratio of two internal nodes in the tree so we will visualize the loadings on the tree using [Empress](https://github.com/biocore/empress).\n",
    "\n",
    "The log-ratio can then be taken with [Qurro](https://github.com/biocore/qurro) to explore the feature loading partitions (arrows) in this biplot as a log-ratio of the original table counts. This allows us to relate these low-dimensional representations back to our original data. Additionally, log-ratios provide a nice set of data points for additional analysis such as LME models. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiime2.plugins.empress.actions import community_plot\n",
    "\n",
    "joint_tree_biplot = community_plot(ctf_results.counts_by_node_tree,\n",
    "                                   ctf_results.subject_table,\n",
    "                                   q2.Metadata(mf),\n",
    "                                   pcoa = ctf_results.subject_biplot,\n",
    "                                   feature_metadata = q2.Metadata(phylo_ctf_taxonomy_and_loadings),\n",
    "                                   number_of_features=50,\n",
    "                                   ignore_missing_samples=False,\n",
    "                                   filter_extra_samples=False,\n",
    "                                   filter_missing_features=True)\n",
    "joint_tree_biplot.visualization.save('IBD-2538/ctf-results/empress.qzv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IBD-2538/ctf-results/qurro.qzv'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiime2.plugins.qurro.actions import loading_plot\n",
    "\n",
    "# run Qurro\n",
    "qurro_plot = loading_plot(ctf_results.subject_biplot, ctf_results.counts_by_node,\n",
    "                          metadata,\n",
    "                          feature_metadata=ctf_results.t2t_taxonomy.view(q2.Metadata))\n",
    "# save visual\n",
    "qurro_plot.visualization.save('IBD-2538/ctf-results/qurro.qzv')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this visualization we can see that the Crohn's subjects clearly separate from the healthy controls in the ordination on the right. We can also observe, by adding a barplot of the PC2 loadings (see the [empress tutorial](https://github.com/biocore/empress) for more info on how to make these plots). \n",
    "\n",
    "![image.png](etc/ctf_empress_plot_one2.png)\n",
    "\n",
    "One additional benefit to having the phylogeny is the ability to explore ratios between phylogenetic paritions. For example the nodes n531 and n135/n119 which represent Faecalibacterium spp. and Akkermansia spp. respectively. \n",
    "\n",
    "\n",
    "![image.png](etc/ctf_empress_plot_two2.png)\n",
    "\n",
    "![image.png](etc/ctf_empress_plot_three2.png)\n",
    "\n",
    "We can then obtain the log-ratio based between these nodes from the qurro visualization.\n",
    "\n",
    "\n",
    "![image.png](etc/ctf_phylo_qurro.png)\n",
    "\n",
    "\n",
    "We can then export the log-ratio output and plot the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>timepoint</th>\n",
       "      <th>host_subject_id</th>\n",
       "      <th>ibd</th>\n",
       "      <th>Current_Natural_Log_Ratio</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>#SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2538.1000102</th>\n",
       "      <td>25</td>\n",
       "      <td>s1000100</td>\n",
       "      <td>Control</td>\n",
       "      <td>4.663439</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2538.1000104000004</th>\n",
       "      <td>6</td>\n",
       "      <td>s1000100</td>\n",
       "      <td>Control</td>\n",
       "      <td>0.470004</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                    timepoint host_subject_id      ibd  \\\n",
       "#SampleID                                                \n",
       "2538.1000102        25         s1000100        Control   \n",
       "2538.1000104000004  6          s1000100        Control   \n",
       "\n",
       "                    Current_Natural_Log_Ratio  \n",
       "#SampleID                                      \n",
       "2538.1000102        4.663439                   \n",
       "2538.1000104000004  0.470004                   "
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# import log-ratio data\n",
    "metadata_one = pd.read_csv('IBD-2538/data/metadata.tsv',\n",
    "                           sep='\\t', index_col=0)\n",
    "# import rest of the metadata\n",
    "metadata_two = pd.read_csv('IBD-2538/ctf-results/sample_plot_data.tsv',\n",
    "                           sep='\\t', index_col=0)[['Current_Natural_Log_Ratio']]\n",
    "# merge the data\n",
    "log_ratio_metdata = pd.concat([metadata_two, metadata_one], axis=1)\n",
    "# ensure no duplicate columns\n",
    "log_ratio_metdata = log_ratio_metdata.dropna(subset=['Current_Natural_Log_Ratio'])\n",
    "log_ratio_metdata.index = log_ratio_metdata.index.astype(str)\n",
    "# export in QIIME2 format\n",
    "log_ratio_metdata = log_ratio_metdata[['timepoint','host_subject_id',\n",
    "                                       'ibd','Current_Natural_Log_Ratio']]\n",
    "\n",
    "log_ratio_metdata.index.name = '#SampleID'\n",
    "log_ratio_metdata.to_csv('IBD-2538/ctf-results/merged_sample_plot_data.tsv', sep='\\t')\n",
    "log_ratio_metdata.head(2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "As you can see above the metadata now has the added column of `Current_Natural_Log_Ratio` from Qurro. So now we will continue to explore this log-ratio by first plotting it explicitly over time with q2-longitudinal.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'IBD-2538/ctf-results/log_ratio_plot.qzv'"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiime2.plugins.longitudinal.actions import (volatility, linear_mixed_effects)\n",
    "\n",
    "# make a time series plot of log-ratio\n",
    "temporal_plot = volatility(q2.Metadata(log_ratio_metdata),\n",
    "                           'timepoint',\n",
    "                           individual_id_column='host_subject_id',\n",
    "                           default_group_column='ibd',\n",
    "                           default_metric='Current_Natural_Log_Ratio')\n",
    "temporal_plot.visualization.save('IBD-2538/ctf-results/log_ratio_plot.qzv')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This demonstrates that we can recreate the separation by IBD that we saw in both the `subject_biplot` & `state_subject_ordination`, allowing us to associate specific taxa (in the numerator or denominator) with a particular phenotype.\n",
    "\n",
    "![image.png](etc/ctf_phylo-lr-vol.png)\n",
    "\n",
    "We can test the statistical power of this log-ratio to differentiate samples by IBD status using a linear mixed effects (LME) through q2-longitudinal. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cmartino/miniconda3/envs/qiime2-2021.2.songbird.mmvec/lib/python3.6/site-packages/q2_longitudinal/_longitudinal.py:291: UserWarning: This is only a warning, and the results of this action are still valid. The column name \"predicted Current_Natural_Log_Ratio\" already exists in your metadata file. Any \"raw\" metadata that can be downloaded from the resulting visualization will contain overwritten values for this metadata column, not the original values.\n",
      "  warnings.warn(warning, UserWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'IBD-2538/ctf-results/lme_log_ratio.qzv'"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run LME model on log-ratio\n",
    "lme_plot = linear_mixed_effects(q2.Metadata(log_ratio_metdata),\n",
    "                                'timepoint',\n",
    "                                individual_id_column='host_subject_id',\n",
    "                                group_columns='ibd',\n",
    "                                metric='Current_Natural_Log_Ratio')\n",
    "lme_plot.visualization.save('IBD-2538/ctf-results/lme_log_ratio.qzv')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this LME model we can see that indeed the IBD grouping is significant across time. \n",
    "\n",
    "![image.png](etc/ctf_lme-logratio-phylo.png)\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}