ipynb/tutorials/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 [DEICODE](https://library.qiime2.org/plugins/q2-deicode).\n",
    "\n",
    "\n",
    "## Introduction \n",
    "\n",
    "In this tutorial you will learn how to interpret and perform Robust Aitchison PCA through QIIME 2. The focus of this tutorial is compositional beta diversity. There are many beta diversity metrics that have been proposed, all with varying benefits on varying data structures. However, presence/absence metric often prove to give better results than those that rely on abundances (i.e. unweighted vs. weighted UniFrac). One component of this phenomenon is the interpretation of relative abundances can provide spurious results (see [the differential abundance analysis introduction](https://docs.qiime2.org/2019.1/tutorials/gneiss/). One solution to this problem is to use a compositional distance metric such as Aitchison distance. \n",
    "\n",
    "\n",
    "As a toy example let’s build three taxa. These three taxa represent common distributions we see in microbiome datasets. Where the first taxon is increasing exponentially across samples, this is a trend that we would be interested in. However, taxon 2 and 3 have much higher counts and taxon 3 is randomly fluctuating across samples.  \n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/7/72ebdf6a3303ce0a5850ce52a46befac564cc26d.png)\n",
    "\n",
    "In our distances below we have Euclidean, Bray-Curtis, Jaccard, and Aitchison distances (from left to right). We can see that the abundance based metrics Euclidean and Bray-Curtis are heavily influenced by the abundance of taxon 3 and seem to randomly fluctuate. In the presence/absence metric, Jaccard, we see that the distance saturates to one very quickly. However, in the Aitchison distance we see a linear curve representing taxon 1. The reason the distance is linear is because Aitchison distance relies on log transforms (the log of an exponential taxon 1 is linear). \n",
    "\n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/b/bc002a51edcd3e34cba1874a6aa97d7d08b6c0b5.png)\n",
    "\n",
    "From this toy example, it is clear that Aitchison distance better accounts for the proportions. However, we made the unrealistic assumption in our toy example that there were no zero counts. In real microbiome datasets there are a large number of zeros (i.e. sparsity). Sparsity complicates log ratio transformations because the log-ratio of zero is undefined. To solve this pseudo counts, that can skew results, are commonly used (see [Naught all zeros in sequence count data are the same](https://www.biorxiv.org/content/10.1101/477794v1)). \n",
    "\n",
    "Robust Aitchison PCA solves this problem in two steps:\n",
    "\n",
    "**1.** Compostional preprocessing using the centered log ratio transform on only the non-zero values of the data (no pseudo count)\n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/4/43fe1323791b5cea419e0973b8983621dbf31a20.gif)\n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/1/13b8c6f415d6ab10c81dec1a27f1f24079be398f.gif)\n",
    "\n",
    "**2.** Dimensionality reduction through Robust PCA on only the non-zero values of the data ( [matrix completion]( https://arxiv.org/pdf/0906.2027.pdf)). \n",
    "\n",
    "![](//forum-qiime2-org.s3.dualstack.us-west-2.amazonaws.com/original/2X/a/a327d5600f68b96457c227c660f533e94ee68341.gif)\n",
    "\n",
    "\n",
    "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 Robust Aitchison PCA via DEICODE 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"
   ]
  },
  {
   "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]), we will generate our beta diversity ordination file. There are a few parameters to DEICODE that we may want to consider. The first is filtering cutoffs, these are `p-min-feature-count`, `p-min-sample-count`, and `p-min-feature-frequency`. Both p-min-feature-count and p-min-sample-count accept integer values and remove feature or samples, respectively, with sums below this cutoff. The feature cut-off is useful in the case that features with very low total counts among all samples represent contamination or chimeric sequences. The sample cut off is useful for the case that some sample received very few reads relative to other samples. The p-min-feature-frequency can be useful to remove features that only appear in a small portion of samples, which may be difficult to further asses using tools like [Qurro](https://github.com/biocore/qurro).\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 DEICODE. \n",
    "\n",
    "The other two parameters are --p-rank and --p-iterations. These parameters should rarely have to change from the default. However, the minimum value of --p-rank can be 1 and the maximum recommended value is 10. Similarly, the minimum value of --p-iterations is 1 and is recommended to be below 500. The default value for p-min-feature-frequency is zero, meaning it will not be used.\n",
    "\n",
    "The other main parameter of the DEICODE is the number of components to use (i.e. the rank). DEICODE relies on a low-rank assumption and therefore it is recommended to choose a value between 2 and 10. If you are unsure of what value to use for n-components the command `auto-rpca` will estimate it for you. In all other aspects the auto-rpca is the same as the `rpca` command.\n",
    "\n",
    "Now that we understand the acceptable parameters, we are ready to run DEICODE.  \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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved PCoAResults % Properties('biplot') to: qiime2-moving-pictures-tutorial/ordination.qza\u001b[0m\r\n",
      "\u001b[32mSaved DistanceMatrix to: qiime2-moving-pictures-tutorial/distance.qza\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime deicode auto-rpca \\\n",
    "    --i-table qiime2-moving-pictures-tutorial/table.qza \\\n",
    "    --p-min-feature-count 10 \\\n",
    "    --p-min-sample-count 500 \\\n",
    "    --o-biplot qiime2-moving-pictures-tutorial/ordination.qza \\\n",
    "    --o-distance-matrix qiime2-moving-pictures-tutorial/distance.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our ordination file, with type (PCoAResults % Properties(['biplot'])), we are ready to visualize the results. This can be done using the [Emperor](https://docs.qiime2.org/2019.1/plugins/available/emperor/) biplot functionality. In this case we will include metadata for our features (optional) and our samples (required). \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved Visualization to: qiime2-moving-pictures-tutorial/biplot.qzv\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime emperor biplot \\\n",
    "    --i-biplot qiime2-moving-pictures-tutorial/ordination.qza \\\n",
    "    --m-sample-metadata-file qiime2-moving-pictures-tutorial/sample-metadata.tsv \\\n",
    "    --m-feature-metadata-file qiime2-moving-pictures-tutorial/taxonomy.qza \\\n",
    "    --o-visualization qiime2-moving-pictures-tutorial/biplot.qzv \\\n",
    "    --p-number-of-features 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Biplots are exploratory visualization tools that allow us to represent the features (i.e. taxonomy or OTUs)  that strongly influence the principal component axis as arrows. The interpretation of the compositional biplot differs slightly from classical biplot interpretation (we can view the qzv file at [view.qiime2](https://view.qiime2.org). 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",
    "![](etc/img12.png)\n",
    "\n",
    "Then by scrolling down the color selections, colors for arrows can be chosen based on taxonomy.\n",
    "\n",
    "![](etc/img13.png)\n",
    "\n",
    "Finally by toggling the colors on the legend, we can add custom coloring by the phylum. \n",
    "\n",
    "![](etc/img14.png)\n",
    "\n",
    "For a more detailed description see this discussion on the [QIIME2 forum](https://forum.qiime2.org/t/deicode-installation-error/8702/7?u=cmartino)\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": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved Visualization to: qiime2-moving-pictures-tutorial/BodySite_significance.qzv\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime diversity beta-group-significance \\\n",
    "    --i-distance-matrix qiime2-moving-pictures-tutorial/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/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 BodySite_significance.qzv at [view.qiime2](https://view.qiime2.org).\n",
    "\n",
    "![](etc/img11.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the biplot and PERMANOVA results that gut is very different from the skin samples. Next we can use [qurro](https://github.com/biocore/qurro) to explore log-ratios of the microbes highlighted by DEICODE. 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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mSaved Visualization to: qiime2-moving-pictures-tutorial/qurro_plot.qzv\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!qiime qurro loading-plot \\\n",
    "    --i-ranks qiime2-moving-pictures-tutorial/ordination.qza \\\n",
    "    --i-table qiime2-moving-pictures-tutorial/table.qza \\\n",
    "    --m-sample-metadata-file qiime2-moving-pictures-tutorial/sample-metadata.tsv \\\n",
    "    --m-feature-metadata-file qiime2-moving-pictures-tutorial/taxonomy.qza \\\n",
    "    --o-visualization qiime2-moving-pictures-tutorial/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 Bacteroides (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 Bacteroides 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 DEICODE 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/img15.png)\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
}