ipynb/tutorials/Phylogenetic-RPCA-moving-pictures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial \n",
    "\n",
    "**Note**: This guide assumes you have installed QIIME 2 using one of the procedures in the [install documents](https://docs.qiime2.org/2019.1/install/) and have installed [gemelli](https://github.com/biocore/gemelli).\n",
    "\n",
    "\n",
    "## Introduction \n",
    "\n",
    "In this tutorial you will learn how to interpret and perform Phylogenetic Robust Aitchison PCA through QIIME 2. This tutorial builds from the Robust Aitchison PCA tutorial, which should be completed first.\n",
    "\n",
    "In order to incorporate phylogenetic information we integrate the FastUniFrac method (see [here](https://www.nature.com/articles/ismej200997)) with the robust centered log ratio transform (rclr). The phylogenetic tree is denoted as $P\\langle\\Upsilon, \\Lambda\\rangle$ where the nodes of the tree are $\\Upsilon=v_{1}, v_{2}, \\ldots v_{a}$ and the branch weights $ A=\\epsilon_{1}, \\epsilon_{2}, \\ldots \\epsilon_{b}$. The total weight is defined as the sum of the branch lengths $W\\left(\\epsilon_{i}\\right)$ which is vectorized and defined as $V_{l}$. In Fast UniFrac the distance between sample $x$ and $x^{\\prime}$ is given as \n",
    "\n",
    "$$\n",
    "d_{U}\\left(x, x^{\\prime}\\right)=\\frac{\\sum_{i=1}^{m} V_{l} \\cdot\\left(x \\oplus x^{\\prime}\\right)}{\\sum_{i=1}^{m} V_{l} \\cdot\\left(x \\vee x^{\\prime}\\right)}\n",
    "$$\n",
    "\n",
    "Updating the rclr in the previous tutorial to include the FastUniFrac transformation can be given as\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    " y_{aj }= \\\n",
    "\\log( V_{l} \\cdot x_{aj}) -\\frac{1}{ \\vert  \\Omega _{V_{l} \\cdot x_{a.}} \\vert } \\sum _{k \\in  \\Omega _{V_{l} \\cdot x_{a.}}}^{}x_{k} - \\frac{1}{ \\vert  \\Omega _{x_{.j}} \\vert } \\sum _{a \\in  \\Omega _{x_{.j}}}^{}x_{k} \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The phylogenetic-RCLR can then be used with the Robust Aitchison PCA algorithm. To demonstrate this in action we will run an example dataset below, where the output can be viewed as a compositional biplot through emperor. \n",
    "\n",
    "## Example \n",
    "\n",
    "\n",
    "In this example we will use Phylogenetic Robust Aitchison PCA via gemelli on the “Moving Pictures” tutorial, if you have not yet completed the tutorial it can be found [here](https://docs.qiime2.org/2019.1/tutorials/moving-pictures/). The dataset consists of human microbiome samples from two individuals at four body sites at five timepoints, the first of which immediately followed antibiotic usage ([Caporaso et al. 2011](https://www.ncbi.nlm.nih.gov/pubmed/21624126)). If you have completed this tutorial run the following command and skip the download section.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "##### Table [view](https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2019.1%2Fdata%2Ftutorials%2Fmoving-pictures%2Ftable.qza) | [download](https://docs.qiime2.org/2019.1/data/tutorials/moving-pictures/table.qza)\n",
    "**save as:** table.qza \n",
    "\n",
    "##### Sample Metadata [download](https://data.qiime2.org/2019.1/tutorials/moving-pictures/sample_metadata.tsv)\n",
    "**save as:** sample-metadata.tsv\n",
    "\n",
    "##### Feature Metadata  [view](https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2019.1%2Fdata%2Ftutorials%2Fmoving-pictures%2Ftaxonomy.qza) | [download](https://docs.qiime2.org/2019.1/data/tutorials/moving-pictures/taxonomy.qza)\n",
    "**save as:** taxonomy.qza\n",
    "\n",
    "##### Rooted Phylogeny  [view](https://view.qiime2.org/?src=https%3A%2F%2Fdocs.qiime2.org%2F2019.1%2Fdata%2Ftutorials%2Fmoving-pictures%2Frooted-tree.qza) | [download](https://docs.qiime2.org/2019.1/data/tutorials/moving-pictures/rooted-tree.qza)\n",
    "**save as:** rooted-tree.qza\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cd qiime2-moving-pictures-tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using table.qza, of the type raw count table (FeatureTable[Frequency]), and rooted-tree.qza of type Phylogeny[Rooted], we will generate our beta diversity ordination file. The first tutorial covers many of the  parameters included in gemelli that we may want to consider. Phylogenetic RPCA included one additional parameter called `--p-min-depth`. This parameter will remove internal nodes of the tree that have less than the input number of descendants. This can be a useful method to prune large trees with high redundancy. \n",
    "\n",
    "**Note:** it is _not_ recommended to bin your features by taxonomic assignment (i.e. by genus level). \n",
    "**Note:** it is _not_ recommended to rarefy your data before using gemelli. \n",
    "\n",
    "Now that we understand the acceptable parameters, we are ready to run gemelli.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mQIIME is caching your current deployment for improved performance. This may take a few moments and should only happen once per deployment.\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime dev refresh-cache"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved PCoAResults % Properties('biplot') to: qiime2-moving-pictures-tutorial/phylo-ordination.qza\u001b[0m\n",
      "\u001b[32mSaved DistanceMatrix to: qiime2-moving-pictures-tutorial/phylo-distance.qza\u001b[0m\n",
      "\u001b[32mSaved Phylogeny[Rooted] to: qiime2-moving-pictures-tutorial/phylo-tree.qza\u001b[0m\n",
      "\u001b[32mSaved FeatureTable[Frequency] to: qiime2-moving-pictures-tutorial/phylo-table.qza\u001b[0m\n",
      "\u001b[32mSaved FeatureData[Taxonomy] to: qiime2-moving-pictures-tutorial/phylo-taxonomy.qza\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "!qiime gemelli phylogenetic-rpca-with-taxonomy \\\n",
    "    --i-table qiime2-moving-pictures-tutorial/table.qza \\\n",
    "    --i-phylogeny qiime2-moving-pictures-tutorial/rooted-tree.qza \\\n",
    "    --m-taxonomy-file qiime2-moving-pictures-tutorial/taxonomy.qza \\\n",
    "    --p-min-feature-count 10 \\\n",
    "    --p-min-sample-count 500 \\\n",
    "    --o-biplot qiime2-moving-pictures-tutorial/phylo-ordination.qza \\\n",
    "    --o-distance-matrix qiime2-moving-pictures-tutorial/phylo-distance.qza \\\n",
    "    --o-counts-by-node-tree qiime2-moving-pictures-tutorial/phylo-tree.qza \\\n",
    "    --o-counts-by-node qiime2-moving-pictures-tutorial/phylo-table.qza \\\n",
    "    --o-t2t-taxonomy qiime2-moving-pictures-tutorial/phylo-taxonomy.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the ordination and distance that RPCA produces, phylogenetic RPCA produces a tree that was used, the table expanded to include the internal nodes, and a taxonomy file with lowest common ancestor taxonomy filled. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our ordination file, with type (PCoAResults % Properties(['biplot'])), we are ready to visualize the results. In this tutorial we can view both the biplot and the phylogeny at the same time with [empress](https://github.com/biocore/empress). In this case we will include metadata for our features (optional) and our samples (required). \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved Visualization to: qiime2-moving-pictures-tutorial/phylo-empress.qzv\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime empress community-plot\\\n",
    "    --i-tree qiime2-moving-pictures-tutorial/phylo-tree.qza\\\n",
    "    --i-feature-table qiime2-moving-pictures-tutorial/phylo-table.qza\\\n",
    "    --i-pcoa qiime2-moving-pictures-tutorial/phylo-ordination.qza\\\n",
    "    --m-sample-metadata-file qiime2-moving-pictures-tutorial/sample-metadata.tsv\\\n",
    "    --m-feature-metadata-file qiime2-moving-pictures-tutorial/phylo-taxonomy.qza\\\n",
    "    --p-filter-missing-features\\\n",
    "    --p-number-of-features 50\\\n",
    "    --o-visualization qiime2-moving-pictures-tutorial/phylo-empress.qzv\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Biplots are exploratory visualization tools that allow us to represent the features (i.e. taxonomy, OTUs, or internal nodes of the phylogeny)  that strongly influence the principal component axis as arrows. The interpretation of the compositional biplot differs slightly from classical biplot interpretation. The previous tutorial describes how to interpret the biplot in more depth.The important features with regard to sample clusters are not a single arrow but by the log ratio between features represented by arrows pointing in different directions.  To effectively use Emperor we fist will color the samples by BodySite. \n",
    "\n",
    "Furthermore, unlike the RPCA output, Phylo-RPCA includes internal nodes in the arrows with taxonomy labeled by LCA. Following the similar methodology in the previous tutorial we can observed that BodySite seems to explain the clusters well. \n",
    "\n",
    "![](etc/img21.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this visualization we noticed that BodySite seems to explain the clusters well. We can run [PERMANOVA](https://docs.qiime2.org/2019.1/plugins/available/diversity/beta-group-significance/) on the distances to get a statistical significance for this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved Visualization to: qiime2-moving-pictures-tutorial/phylo-BodySite_significance.qzv\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime diversity beta-group-significance \\\n",
    "    --i-distance-matrix qiime2-moving-pictures-tutorial/phylo-distance.qza \\\n",
    "    --m-metadata-file qiime2-moving-pictures-tutorial/sample-metadata.tsv \\\n",
    "    --m-metadata-column BodySite \\\n",
    "    --p-method permanova \\\n",
    "    --o-visualization qiime2-moving-pictures-tutorial/phylo-BodySite_significance.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed we can now see that the clusters we saw in the biplot were significant by viewing the phylo-BodySite_significance.qzv at [view.qiime2](https://view.qiime2.org).\n",
    "\n",
    "![](etc/img19.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the biplot and PERMANOVA results that gut is very different from the skin samples. \n",
    "\n",
    "Next we can color the internal nodes and sOTUs by their LCA and directly predicted taxonomy respectively. \n",
    "\n",
    "![](etc/img22.png)\n",
    "\n",
    "By double clicking on the arrows associated with the different body sites we can explore where they fall in the phylogeny. For example Faecalibacterium (g) associated with gut\n",
    "\n",
    "![](etc/img23.png)\n",
    "\n",
    "\n",
    "and Streptococcus (g) associated with the skin/tongue body site samples.\n",
    "\n",
    "![](etc/img24.png)\n",
    "\n",
    "\n",
    "Next we can use [qurro](https://github.com/biocore/qurro) to explore log-ratios of the phylogenetic partitions of microbes highlighted by gemelli. For more about why log-ratios are useful you may want to read [\"Establishing microbial composition measurement standards with reference frames\"](https://www.nature.com/articles/s41467-019-10656-5).\n",
    "\n",
    "**Note: It is important that log-ratios are not taken between nodes that shared descendants**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved Visualization to: qiime2-moving-pictures-tutorial/phylo-qurro_plot.qzv\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime qurro loading-plot \\\n",
    "    --i-ranks qiime2-moving-pictures-tutorial/phylo-ordination.qza \\\n",
    "    --i-table qiime2-moving-pictures-tutorial/phylo-table.qza \\\n",
    "    --m-sample-metadata-file qiime2-moving-pictures-tutorial/sample-metadata.tsv \\\n",
    "    --m-feature-metadata-file qiime2-moving-pictures-tutorial/phylo-taxonomy.qza \\\n",
    "    --o-visualization qiime2-moving-pictures-tutorial/phylo-qurro_plot.qzv\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two taxa groups whose arrows seem to be directly opposed with relation to the BodySite grouping is Faecalibacterium (associated with gut) and Streptococcus (associated with skin and oral). We can use Qurro to explore this relationship. To make a log-ratio we can filter by taxa who contain Faecalibacterium in the numerator and Streptococcus in the denominator of the log-ratio. Those features will then be summed according to thier taxonomic labels and used in the log-ratio. In Qurro the axis one loadings (or another axis) from gemelli are highlighted by if they are contained in the numerator or denominator. The log-ratio plot is contained on the left and can be visualized as a scatter or box-plot. From this it is clear these taxa can separate our BodySite groupings. The tsv file can be exported and a t-test by BodySite on the log-ratos could confirm this observation.\n",
    "\n",
    "![](etc/img26.png)\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
}