Parkinson’s Mouse Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parkinson’s Mouse Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**_note:_ This guide assumes you have QIIME 2 installed (e.g. using [this procedure](https://docs.qiime2.org/2021.2/install/native/)). To execute the script properly, open this notebook in a Jupyter Notebook from within a conda QIIME 2 environment.**\n",
    "\n",
    "**_note:_ This tutorial is an adaptation of the same [tutorial](https://docs.qiime2.org/2021.2/tutorials/pd-mice/) that may be found on the official [QIIME 2 docs website](https://docs.qiime2.org/2021.2/). The original tutorial uses the QIIME 2 CLI interface.**\n",
    "\n",
    "This tutorial will demonstrate a “typical” QIIME 2 analysis of 16S rRNA gene amplicon data, using a set of fecal samples from [humanized mice](https://en.wikipedia.org/wiki/Humanized_mouse). The original study, [Sampson et al, 2016](https://pubmed.ncbi.nlm.nih.gov/27912057/), was designed to determine whether the fecal microbiome contributed to the development of Parkinson’s Disease (PD). Several observation studies showed a difference in the microbiome between PD patients and controls, although the organisms identified across studies were not consistent. However, this was sufficient evidence to suggest that there might be a relationship between PD and the fecal microbiome.\n",
    "\n",
    "To determine whether that relationship was incidental or actually disease associated, a second study was needed. A human cohort study was not feasible; the disease only affects about 1% of the population over 60 years old, PD takes a long time to develop and to be diagnosed, and it would be difficult to determine when to collect the samples. Therefore, a gnotobiotic mouse study was utilized to evaluate the role of the microbiome in the development of PD symptoms. Feces were collected from six donors with Parkinson’s disease and six age- and sex-matched neurologically health controls, and then transplanted into mice who were either predisposed to developing Parkinson’s disease due to a mutation (“aSyn”) or resistant wild type mice (“BDF1”). Mice from different donors were kept in separate cages, but mice from different genetic backgrounds were co-housed. The mice were followed for 7 weeks to see if they developed symptoms of Parkinson’s disease.\n",
    "\n",
    "We’ll look a subset of data from two human donors (one healthy and one with PD) whose samples were each transplanted into three separate cages of mice from the susceptible genotype. For this tutorial, a subset of the metadata has been prepared, and the sequences have been subsampled to approximately 5000 sequences per sample to allow the tutorial to run in a short time. The sequences for the full study are accessible at EBI with accession [PRJEB17694](https://www.ebi.ac.uk/ena/data/view/PRJEB17694); processed tables from the full study can be downloaded from the [Qiita](https://qiita.ucsd.edu/) database from study 10483."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial will explore the hypothesis that the genetic background of a humanized mouse influences the microbial community. However, we’ll also need to consider other confounders which might drive the shape of the microbiome instead of the mouse genotype."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qiime2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiime2.plugins import (demux, dada2, metadata, feature_table,\n",
    "                            fragment_insertion, diversity, \n",
    "                            longitudinal, phylogeny, sample_classifier,\n",
    "                            feature_classifier, taxa, composition, \n",
    "                            longitudinal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial assumes that you have QIIME 2 installed according to the [installation instructions](https://docs.qiime2.org/2020.6/install/).\n",
    "\n",
    "Before running the tutorial, you will need to make a directory for the tutorial data and navigate into that directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "workdir='/path/to/your/directory/qiime2-pdmice-tutorial/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir -p $workdir\n",
    "!cd $workdir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before starting any analysis, it’s important to be familiar with the metadata. In this study, the metadata file contains 7 columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:26:30--  https://data.qiime2.org/2019.10/tutorials/pd-mice/manifest\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://s3-us-west-2.amazonaws.com/qiime2-data/2019.10/tutorials/pd-mice/manifest [following]\n",
      "--2020-07-28 18:26:31--  https://s3-us-west-2.amazonaws.com/qiime2-data/2019.10/tutorials/pd-mice/manifest\n",
      "Resolving s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)... 52.218.176.0\n",
      "Connecting to s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)|52.218.176.0|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 4640 (4.5K) [binary/octet-stream]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//manifest.tsv’\n",
      "\n",
      "/home/user/Document 100%[===================>]   4.53K  --.-KB/s    in 0s      \n",
      "\n",
      "2020-07-28 18:26:32 (32.2 MB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//manifest.tsv’ saved [4640/4640]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"manifest.tsv\" \\\n",
    "  \"https://data.qiime2.org/2019.10/tutorials/pd-mice/manifest\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:26:32--  https://data.qiime2.org/2019.10/tutorials/pd-mice/sample_metadata.tsv\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://docs.google.com/spreadsheets/d/e/2PACX-1vRovoHFup2t5MxddVNLcJ96f5hkh5DJsx84uqCsTwAkz06i0ibrCI5Tn0xFF15xHDmVXfDYVb2o_YXD/pub?gid=1509704122&single=true&output=tsv [following]\n",
      "--2020-07-28 18:26:33--  https://docs.google.com/spreadsheets/d/e/2PACX-1vRovoHFup2t5MxddVNLcJ96f5hkh5DJsx84uqCsTwAkz06i0ibrCI5Tn0xFF15xHDmVXfDYVb2o_YXD/pub?gid=1509704122&single=true&output=tsv\n",
      "Resolving docs.google.com (docs.google.com)... 216.58.208.206, 2a00:1450:401b:806::200e\n",
      "Connecting to docs.google.com (docs.google.com)|216.58.208.206|:443... connected.\n",
      "HTTP request sent, awaiting response... 307 Temporary Redirect\n",
      "Location: https://doc-0o-6o-sheets.googleusercontent.com/pub/l5l039s6ni5uumqbsj9o11lmdc/s7i4ijh3iq3hltlpj9vs92b3ug/1595953590000/103995680502445084602/*/e@2PACX-1vRovoHFup2t5MxddVNLcJ96f5hkh5DJsx84uqCsTwAkz06i0ibrCI5Tn0xFF15xHDmVXfDYVb2o_YXD?gid=1509704122&single=true&output=tsv [following]\n",
      "Warning: wildcards not supported in HTTP.\n",
      "--2020-07-28 18:26:33--  https://doc-0o-6o-sheets.googleusercontent.com/pub/l5l039s6ni5uumqbsj9o11lmdc/s7i4ijh3iq3hltlpj9vs92b3ug/1595953590000/103995680502445084602/*/e@2PACX-1vRovoHFup2t5MxddVNLcJ96f5hkh5DJsx84uqCsTwAkz06i0ibrCI5Tn0xFF15xHDmVXfDYVb2o_YXD?gid=1509704122&single=true&output=tsv\n",
      "Resolving doc-0o-6o-sheets.googleusercontent.com (doc-0o-6o-sheets.googleusercontent.com)... 216.58.209.1, 2a00:1450:401b:803::2001\n",
      "Connecting to doc-0o-6o-sheets.googleusercontent.com (doc-0o-6o-sheets.googleusercontent.com)|216.58.209.1|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: unspecified [text/tab-separated-values]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//metadata.tsv’\n",
      "\n",
      "/home/user/Document     [ <=>                ]   4.33K  --.-KB/s    in 0s      \n",
      "\n",
      "2020-07-28 18:26:34 (28.3 MB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//metadata.tsv’ saved [4431]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"metadata.tsv\" \\\n",
    "  \"https://data.qiime2.org/2019.10/tutorials/pd-mice/sample_metadata.tsv\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing metadata as a qiime2 artifact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All data that is used as input to QIIME 2 is in form of QIIME 2 artifacts, which contain information about the type of data and the source of the data. The first thing we need to do is import these sequence data files into a QIIME 2 artifact."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In QIIME 2, all data is structured as an Artifact of a specific semantic type. Artifacts contain the data as well as information about the data, including a record of the original data and the tools used to process it. This allows for better tracking of how you actually got to where you are in your analysis. You can learn more about common QIIME 2 Artifacts and semantic types here.\n",
    "\n",
    "Our samples were amplified using the [EMP 515f-806r](https://www.earthmicrobiome.org/protocols-and-standards/16s/) primers and sequenced on an Illumina MiSeq with a 2x150bp kit. The hypervariable region covered by the primers we used is 290bp long, so with 150bp reads our sequences will be slightly too short to be able to do paired-end analysis downstream. Therefore, we’re going to work with single-end sequences. We will work with a version of the samples which have already been demultiplexed, for example, by the sequencing center. If you need to demultiplex your sequences, the [Moving Pictures tutorial](https://docs.qiime2.org/2020.6/tutorials/moving-pictures/)describes how to demultiplex sequences if they were sequenced using the Earth Microbiome Project protocol.\n",
    "\n",
    "We will import the sequences as `SampleData[SequencesWithQuality]`, which is the demultiplexed single-end sequence format. If we wanted to import paired-sequences, we would specify the semantic type `SampleData[PairedEndSequencesWithQuality]`. We will import the sequences using the sample [manifest format](https://docs.qiime2.org/2020.6/tutorials/importing/#manifest-file), a versatile way to import demultiplexed data in QIIME 2. We create a tab-separated sample manifest file that maps the sample name we want to use in QIIME 2 to the path of the sequence file. The benefit is that the demultiplexed sequence files can be named anything you want; there are not fixed assumptions about the conventions, and the file names do not dictate the final name. When QIIME 2 reads the file, it ignores any line prefixed with the # symbol. The first line that doesn’t contain a # is the header line and must be `sample-id<TAB>absolute-filepath`. The sample order after the header line does not matter. Read more about importing data into QIIME 2 artifacts [here](https://docs.qiime2.org/2020.6/tutorials/importing/) and more about sample metadata formatting requirements [here](https://docs.qiime2.org/2020.6/tutorials/metadata/).\n",
    "\n",
    "Let’s start by downloading the manifest and corresponding sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "manifest = qiime2.Metadata.load(workdir+'/manifest.tsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_metadata = qiime2.Metadata.load(workdir+'/metadata.tsv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:26:34--  https://data.qiime2.org/2019.10/tutorials/pd-mice/demultiplexed_seqs.zip\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://s3-us-west-2.amazonaws.com/qiime2-data/2019.10/tutorials/pd-mice/demultiplexed_seqs.zip [following]\n",
      "--2020-07-28 18:26:35--  https://s3-us-west-2.amazonaws.com/qiime2-data/2019.10/tutorials/pd-mice/demultiplexed_seqs.zip\n",
      "Resolving s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)... 52.218.176.0\n",
      "Connecting to s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)|52.218.176.0|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 21508775 (21M) [application/zip]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//demultiplexed_seqs.zip’\n",
      "\n",
      "/home/user/Document 100%[===================>]  20.51M  2.18MB/s    in 11s     \n",
      "\n",
      "2020-07-28 18:26:47 (1.84 MB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//demultiplexed_seqs.zip’ saved [21508775/21508775]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"demultiplexed_seqs.zip\" \\\n",
    "  \"https://data.qiime2.org/2019.10/tutorials/pd-mice/demultiplexed_seqs.zip\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Archive:  demultiplexed_seqs.zip\n",
      "replace demultiplexed_seqs/10483.recip.539.ASO.PD4.D7_4_L001_R1_001.fastq.gz? [y]es, [n]o, [A]ll, [N]one, [r]ename: ^C\n"
     ]
    }
   ],
   "source": [
    "!unzip demultiplexed_seqs.zip"
   ]
  },
  {
   "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>absolute-filepath</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sample-id</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>recip.220.WT.OB1.D7</th>\n",
       "      <td>$PWD/demultiplexed_seqs/10483.recip.220.WT.OB1...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>recip.290.ASO.OB2.D1</th>\n",
       "      <td>$PWD/demultiplexed_seqs/10483.recip.290.ASO.OB...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>recip.389.WT.HC2.D21</th>\n",
       "      <td>$PWD/demultiplexed_seqs/10483.recip.389.WT.HC2...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>recip.391.ASO.PD2.D14</th>\n",
       "      <td>$PWD/demultiplexed_seqs/10483.recip.391.ASO.PD...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>recip.391.ASO.PD2.D21</th>\n",
       "      <td>$PWD/demultiplexed_seqs/10483.recip.391.ASO.PD...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                       absolute-filepath\n",
       "sample-id                                                               \n",
       "recip.220.WT.OB1.D7    $PWD/demultiplexed_seqs/10483.recip.220.WT.OB1...\n",
       "recip.290.ASO.OB2.D1   $PWD/demultiplexed_seqs/10483.recip.290.ASO.OB...\n",
       "recip.389.WT.HC2.D21   $PWD/demultiplexed_seqs/10483.recip.389.WT.HC2...\n",
       "recip.391.ASO.PD2.D14  $PWD/demultiplexed_seqs/10483.recip.391.ASO.PD...\n",
       "recip.391.ASO.PD2.D21  $PWD/demultiplexed_seqs/10483.recip.391.ASO.PD..."
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "manifest.to_dataframe().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using this manifest format, a sample name can only appear in one line and can only map to one sequencing file per column (one column for single-end, two columns for paired-end). The absolute-filepath for each sample must be an [absolute path](https://en.wikipedia.org/wiki/Path_(computing)#Absolute_and_relative_paths), which specifies the “complete” location of the file. We do that here using the `$PWD` variable, which expands the current directory in absolute terms.\n",
    "\n",
    "We’ll use the manifest to import our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "demux_seqs = qiime2.Artifact.import_data(type='SampleData[SequencesWithQuality]',\n",
    "                                         view_type='SingleEndFastqManifestPhred33V2',\n",
    "                                         view=workdir+'/manifest.tsv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s check the sequences and the sequencing depth of the samples using the `qiime demux summarize command`. It provides information about the number of sequences in each sample, as well as the quality of the sequences.\n",
    "\n",
    "Before running the command, let’s review the help documentation to make sure we understand the arguments for the command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Usage: \u001b[34mqiime demux summarize\u001b[0m [OPTIONS]\r\n",
      "\r\n",
      "  Summarize counts per sample for all samples, and generate interactive\r\n",
      "  positional quality plots based on `n` randomly selected sequences.\r\n",
      "\r\n",
      "\u001b[1mInputs\u001b[0m:\r\n",
      "  \u001b[34m\u001b[4m--i-data\u001b[0m ARTIFACT \u001b[32mSampleData[SequencesWithQuality |\u001b[0m\r\n",
      "    \u001b[32mPairedEndSequencesWithQuality | JoinedSequencesWithQuality]\u001b[0m\r\n",
      "                       The demultiplexed sequences to be summarized.\r\n",
      "                                                                    \u001b[35m[required]\u001b[0m\r\n",
      "\u001b[1mParameters\u001b[0m:\r\n",
      "  \u001b[34m--p-n\u001b[0m INTEGER        The number of sequences that should be selected at\r\n",
      "                       random for quality score plots. The quality plots will\r\n",
      "                       present the average positional qualities across all of\r\n",
      "                       the sequences selected. If input sequences are paired\r\n",
      "                       end, plots will be generated for both forward and\r\n",
      "                       reverse reads for the same `n` sequences.\r\n",
      "                                                              \u001b[35m[default: 10000]\u001b[0m\r\n",
      "\u001b[1mOutputs\u001b[0m:\r\n",
      "  \u001b[34m\u001b[4m--o-visualization\u001b[0m VISUALIZATION\r\n",
      "                                                                    \u001b[35m[required]\u001b[0m\r\n",
      "\u001b[1mMiscellaneous\u001b[0m:\r\n",
      "  \u001b[34m--output-dir\u001b[0m PATH    Output unspecified results to a directory\r\n",
      "  \u001b[34m--verbose\u001b[0m / \u001b[34m--quiet\u001b[0m  Display verbose output to stdout and/or stderr during\r\n",
      "                       execution of this action. Or silence output if\r\n",
      "                       execution is successful (silence is golden).\r\n",
      "  \u001b[34m--citations\u001b[0m          Show citations and exit.\r\n",
      "  \u001b[34m--help\u001b[0m               Show this message and exit.\r\n"
     ]
    }
   ],
   "source": [
    "!qiime demux summarize --help"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the documentation, we should specify the file (Artifact) with the demultiplexed sequences for the `--i-`data argument, since this expects data of semantic type `SampleData[SequencesWithQuality]`. We’ll specify the location we want to save the visualization to by specifying the output path to `--o-visualization`.\n",
    "\n",
    "The help documentation is a good reference for any command, and the first place to look if you’re getting errors, especially errors about parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-rvam5w1j')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 2789d9c9-7620-44b1-af39-554ec387a248>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demux_summary = demux.visualizers.summarize(demux_seqs)\n",
    "demux_summary.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "1. After demultiplexing, which sample has the lowest sequencing depth?\n",
    "2. What is the median sequence length?\n",
    "3. What is the median quality score at position 125?\n",
    "4. If you are working on this tutorial alongside someone else, why does your plot look slightly different from your neighbors? If you aren’t working alongside someone else, try running this command a few times and compare the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sequence quality control and feature table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several ways to construct a feature table in QIIME 2. The first major choice to make is to work with Operational Taxonomic Units (OTUs) or Absolute Sequence Variants (ASVs). OTUs have been widely used in microbiome research since the mid 2010s, and assign sequences to clusters either based on a reference database or de novo assignment. QIIME 2 offers clustering through [q2-vsearch](https://docs.qiime2.org/2020.6/tutorials/otu-clustering/) and [q2-dbOTU] plug-ins, currently.\n",
    "\n",
    "ASVs are a more recent development and provide better resolution in features than traditional OTU-based methods. ASVs can separate features based on differences of a single nucleotide in sequences of 400 bp or more, a resolution not possibly even with 99% identity OTU clustering. QIIME 2 currently offers denoising via [DADA2](https://pubmed.ncbi.nlm.nih.gov/27214047/) (q2-dada2) and [Deblur](https://pubmed.ncbi.nlm.nih.gov/28289731/) (q2-deblur). The major differences in the algorithms and motivation for denoising are nicely described in [Nearing et al, 2018](https://pubmed.ncbi.nlm.nih.gov/30123705/).\n",
    "\n",
    "It is worth noting in either case that denoising to ASVs and clustering to OTUs are separate, but parallel steps. A choice should be made for a single pathway: either denoising or OTU based clustering; it is not recommended to combine the steps.\n",
    "\n",
    "In this tutorial, we’ll denoise with DADA2 (using single-end sequences). Please see the [Atacama Soil tutorial](https://docs.qiime2.org/2020.6/tutorials/atacama-soils/) for an example of using DADA2 on paired-end sequences. For those interested in using Deblur, you can refer to the [Moving Pictures](https://docs.qiime2.org/2020.6/tutorials/moving-pictures/) tutorial and [Alternative methods of read joining](https://docs.qiime2.org/2020.6/tutorials/read-joining/) tutorial for running Deblur on single- and paired-end sequences, respectively.\n",
    "\n",
    "The qiime `dada2.methods.denoise_single` method requires us to set the `trunc_len` parameter. This controls the length of the sequences and should be selected based on a drop in quality scores. In our dataset, the quality scores are relatively evenly distributed along the sequencing run, so we’ll use the full 150 bp sequences. However, the selection of the trim length is a relatively subjective measurement and relies on the decision making capacity of the analyst."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running external command line application(s). This may print messages to stdout and/or stderr.\n",
      "The command(s) being run are below. These commands cannot be manually re-run as they will depend on temporary files that no longer exist.\n",
      "\n",
      "Command: run_dada_single.R /tmp/qiime2-archive-ltlh76m5/53a431ce-a88f-42c5-a5c4-72568d1312d2/data /tmp/tmp6ipwgvsz/output.tsv.biom /tmp/tmp6ipwgvsz/track.tsv /tmp/tmp6ipwgvsz 150 0 2.0 2 Inf consensus 1.0 1 1000000 NULL 16\n",
      "\n"
     ]
    }
   ],
   "source": [
    "denoised_sequences = dada2.methods.denoise_single(demux_seqs,\n",
    "                                                  trunc_len = 150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also review the denoising statistics using the `metadata.visualizers.tabulate` command.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-8sqjruzf')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 831ca847-6755-4f3f-9b8f-eebc4a27126f>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filter_stats = metadata.visualizers.tabulate(denoised_sequences.denoising_stats.view(qiime2.Metadata))\n",
    "filter_stats.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature table summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we finish denoising the data, we can check the results by looking at the summary of the feature table. This will provide us with the counts associated with each sequence and each feature, as well as other useful plots and metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_viz = feature_table.visualizers.summarize(denoised_sequences.table)\n",
    "output_viz.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "1. How many total features remain after denoising?\n",
    "2. Which sample has the highest total count of features? How many sequences did that sample have prior to DADA2 denoising?\n",
    "3. How many samples have fewer than 4250 total features?\n",
    "4. Which features are observed in at least 47 samples?\n",
    "5. Which sample has the fewest features? How many does it have?\n",
    "\n",
    "If you open the denoising summary, can you find the step where the sample with the fewest sequences fails?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating a phylogenetic tree for diversity analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QIIME 2 analysis allows the use of phylogenetic trees for diversity metrics such as Faith’s Phylogenetic Diversity and UniFrac distance as well as feature-based analyses in Gneiss. The tree provides an inherent structure to the data, allowing us to consider an evolutionary relationship between organisms.\n",
    "\n",
    "QIIME 2 offers several ways to construct a phylogenetic tree. For this tutorial, we’re going to create a fragment insertion tree using the `q2-fragment-insertion` plugin. The authors of the fragment insertion plugin suggest that it can outperform traditional alignment based methods based on short Illumina reads by alignment against a reference tree built out of larger sequences. Our command, qiime `fragment-insertion` sepp will use the representative sequences (a `FeatureData[Sequence]` artifact) we generated during denoising to create a phylogenetic tree where the sequences have been inserted into the greengenes 13_8 99% identity reference tree backbone.\n",
    "\n",
    "First, we will download the reference database:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:28:34--  https://data.qiime2.org/2019.10/common/sepp-refs-gg-13-8.qza\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://s3-us-west-2.amazonaws.com/qiime2-data/2019.10/common/sepp-refs-gg-13-8.qza [following]\n",
      "--2020-07-28 18:28:35--  https://s3-us-west-2.amazonaws.com/qiime2-data/2019.10/common/sepp-refs-gg-13-8.qza\n",
      "Resolving s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)... 52.218.185.8\n",
      "Connecting to s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)|52.218.185.8|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 50161069 (48M) [binary/octet-stream]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//sepp-refs-gg-13-8.qza’\n",
      "\n",
      "/home/user/Document 100%[===================>]  47.84M  1.13MB/s    in 42s     \n",
      "\n",
      "2020-07-28 18:29:17 (1.15 MB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//sepp-refs-gg-13-8.qza’ saved [50161069/50161069]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"sepp-refs-gg-13-8.qza\" \\\n",
    "  \"https://data.qiime2.org/2019.10/common/sepp-refs-gg-13-8.qza\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Importing the reference database as a qiime2 artifact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "database_sepp = qiime2.Artifact.load(workdir+'/sepp-refs-gg-13-8.qza')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note\n",
    "\n",
    "This command is resource intensive - if your computation environment supports it, we suggest including an appropriately-set --p-threads parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'denoised_sequences' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-b8869bd7e2d2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mrooted_tree\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mphylogeny\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpipelines\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0malign_to_tree_mafft_fasttree\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdenoised_sequences\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrepresentative_sequences\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'denoised_sequences' is not defined"
     ]
    }
   ],
   "source": [
    "rooted_tree = phylogeny.pipelines.align_to_tree_mafft_fasttree(denoised_sequences.representative_sequences)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_placements = fragment_insertion.methods.sepp(denoised_sequences.representative_sequences, \n",
    "                                                  database_sepp,\n",
    "                                                  threads = 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alpha Rarefaction and Selecting a Rarefaction Depth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a feature table (observation matrix) of ASVs in each sample, and a phylogenetic tree representing those ASVs, so are almost ready to perform various analyses of microbial diversity. However, first we must normalize our data to account for uneven sequencing depth between samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although sequencing depth in a microbiome sample does not directly relate to the original biomass in a community, the relative sequencing depth has a large impact on observed communities [(Weiss et al, 2017)](https://pubmed.ncbi.nlm.nih.gov/28253908/). Therefore, for most diversity metrics, a normalization approach is needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Current best practices suggest the use of rarefaction, a normalization via sub-sampling without replacement. Rarefaction occurs in two steps: first, samples which are below the rarefaction depth are filtered out of the feature table. Then, all remaining samples are subsampled without replacement to get to the specified sequencing depth. It’s both important and sometimes challenging to select a rarefaction depth for diversity analyses. Several strategies exist to figure out an appropriate rarefaction depth - we will primarily consider alpha rarefaction in this tutorial, because it is a data-driven way to approach the problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We’ll use `diversity.visualizers.alpha_rarefaction` to subsample the ASV table at different depths (between `min_depth` and `max_depth`) and calculate the alpha diversity using one or more metrics . When we checked the feature table, we found that the sample with the fewest sequences in the denoised table has 85 features and that the sample with the most has 4996 features. We want to set a maximum depth close to the maximum number of sequences. We also know that if we look at a sequencing depth around 4250 sequences per sample, we’ll be looking at information from 22 samples. So, let’s set this as our maximum sequencing depth.\n",
    "\n",
    "At each sampling depth, 10 rarefied tables are usually calculated to provide an error estimate, although this can be adjusted using the `--p-iterations` parameter. We can check and see if there is a relationship between the alpha diversity and metadata by specifying the metadata file for the `metadata = sample_metadata` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized_data = diversity.visualizers.alpha_rarefaction(denoised_sequences.table, \n",
    "                                                          metadata = sample_metadata,\n",
    "                                                          min_depth = 10,\n",
    "                                                          max_depth = 4250)\n",
    "                                                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized_data_1 = diversity.visualizers.alpha_rarefaction(denoised_sequences.table,\n",
    "                                                          4250,\n",
    "                                                          metadata = sample_metadata,\n",
    "                                                          min_depth = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-qn74rbc9')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: f7700ee4-3e8c-4221-899e-ab9171100a58>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normalized_data_1.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The visualization file will display two plots. The upper plot will display the alpha diversity (observed OTUs or shannon) as a function of the sampling depth. This is used to determine whether the richness or evenness has saturated based on the sampling depth. The rarefaction curve should “level out” as you approach the maximum sampling depth. Failure to do so, especially with a diversity-only metric such as observed OTUs or Faith’s PD diversity, may indicate that the richness in the samples has not been fully saturated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second plot shows the number of samples in each metadata category group at each sampling depth. This is useful to determine the sampling depth where samples are lost, and whether this may be biased by metadata column group values. Remember that rarefaction is a two-step process and samples that do not meet the rarefaction depth are filtered out of the table. We can use the curves to look at the number of samples by different metadata columns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you’re still unsure of the rarefaction depth, you can also use the sample summary to look at which samples are lost by supplying sample metadata to the feature table summary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Question\n",
    "\n",
    "Start by opening the alpha rarefaction visualization.\n",
    "\n",
    "1. Are all metadata columns represented in the visualization? If not, which columns were excluded and why?\n",
    "\n",
    "2. Which metric shows saturation and stabilization of the diversity?\n",
    "\n",
    "3. Which mouse genetic background has higher diversity, based on the curve? Which has shallower sampling depth?\n",
    "\n",
    "Now, let’s check the feature table summary.\n",
    "\n",
    "4. What percentage of samples are lost if we set the rarefaction depth to 2500 sequences per sample?\n",
    "\n",
    "5. Which mice did the missing samples come from?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we’ve looked through the data, we need to select a rarefaction depth. In general, selecting a rarefaction depth is a subjective process that requires that analyst’s discretion. Selecting a rarefaction depth is an exercise in minimizing sequence loss while maximizing the sequences retained for diversity analysis. For high-biomass samples (fecal, oral, etc), a general best estimate is a rarefaction depth of no less than 1000 sequences per sample. In low biomass samples where sequencing is shallower, a lower rarefaction depth may be selected although it’s important to keep in mind that the diversity measurements on these samples will be quite noisy and the overall quality will be low."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checkpoint\n",
    "\n",
    "Based on the current rarefaction curve and sample summary, what sequencing depth would you pick? Why?\n",
    "\n",
    "In this case, we can retain 47 samples with a rarefaction depth of 2000 sequences/sample.\n",
    "\n",
    "Based on the sequencing depth and distribution of samples, we’ll use 2000 sequences/sample for this analysis. This will let us keep 47 of 48 high quality samples (discarding the one sample with sequencing depth below 1000 sequences/sample)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diversity analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step in hypothesis testing in microbial ecology is typically to look at within- (alpha) and between-sample (beta) diversity. We can calculate diversity metrics, apply appropriate statistical tests, and visualize the data using the `q2-diversity` plugin.\n",
    "\n",
    "We’ll start by using the qiime `diversity.pipelines.core_metrics_phylogenetic` method, which ratifies the input feature table, calculates several commonly used alpha- and beta-diversity metrics, and produces principal coordinate analysis (PCoA) visualizations in Emperor for the beta diversity metrics. By default, the metrics computed are:\n",
    "\n",
    "Alpha Diversity\n",
    "\n",
    "- Shannon’s diversity index\n",
    "\n",
    "- Observed Features\n",
    "\n",
    "- Faith’s phylogenetic diversity\n",
    "\n",
    "- Pielou’s evenness\n",
    "\n",
    "Beta Diversity\n",
    "\n",
    "- Jaccard distance\n",
    "\n",
    "- Bray-Curtis distance\n",
    "\n",
    "- Unweighted UniFrac distance\n",
    "\n",
    "- Weighted UniFrac distance\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note\n",
    "\n",
    "🏗👷 Some descriptions are changing in QIIME 2’s diversity tools. The phrase “observed otus” is being replaced with “observed features”, because “features” better describes the different ways in which users work with non-taxonomic features. This will affect both documentation and (in places) the names of command arguments. You will see both phrases, but these measures of diversity are identical under the hood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a very good discussion of diversity metrics and their meanings in a [forum post by Stephanie Orchanian](https://forum.qiime2.org/t/alpha-and-beta-diversity-explanations-and-commands/2282).\n",
    "\n",
    "The qiime `diversity.pipelines.core_metrics_phylogenetic` method wraps several other methods, and it’s worthwhile to note that the steps can also be executed independently.\n",
    "\n",
    "One important consideration for diversity calculations is the rarefaction depth. Above, we used the alpha rarefaction visualization and the sample summary visualization to pick a rarefaction depth. So, for these analyses, we’ll use a depth of 2000 sequences per sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'tree_placements' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-26-f86ac10522b1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m core_metrics_results = diversity.pipelines.core_metrics_phylogenetic(denoised_sequences.table,\n\u001b[0;32m----> 2\u001b[0;31m                                                                      \u001b[0mtree_placements\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtree\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m                                                                      \u001b[0msampling_depth\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2000\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m                                                                      metadata = sample_metadata)\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'tree_placements' is not defined"
     ]
    }
   ],
   "source": [
    "core_metrics_results = diversity.pipelines.core_metrics_phylogenetic(denoised_sequences.table,\n",
    "                                                                     tree_placements.tree,\n",
    "                                                                     sampling_depth = 2000,\n",
    "                                                                     metadata = sample_metadata)\n",
    "                                                                    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Where did we get the value `2000` from? Why did we pick that?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alpha diversity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Alpha diversity asks whether the distribution of features within a sample (or groups of samples) differs between different conditions. The comparison makes no assumptions about the features that are shared between samples; two samples can have the same alpha diversity and not share any features. The rarefied alpha diversity produced by `diversity` is a univariate, continuous value and can be tested using common non-parametric statistical tests.\n",
    "\n",
    "We can test our covariates of interest against Faith’s phylogenetic diversity and Pielou’s evenness value by running:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "faith_pd_statistics = diversity.visualizers.alpha_group_significance(core_metrics_results.faith_pd_vector,\n",
    "                                                                     sample_metadata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-s53w3dcr')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 46dabba6-d534-4864-9e9e-00d54ab928c8>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "faith_pd_statistics.visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "evenness_statistics = diversity.visualizers.alpha_group_significance(core_metrics_results.evenness_vector,\n",
    "                                                                     sample_metadata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "1. Is there a difference in evenness between genotype? Is there a difference in phylogenetic diversity between genotype?\n",
    "\n",
    "2. Based on the group significance test, is there a difference in phylogenetic diversity by genotype? Is there a difference based on the donor?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we had a continuous covariate that we thought was associated with the alpha diversity, we could test that using qiime diversity alpha-correlation. However, the only continuous variable in this dataset is the days_since_transplant.\n",
    "\n",
    "In some experiments, multiple interacting factors may impact alpha diversity together. If our alpha diversity estimates follow a normal distribution, we may use analysis of variance (ANOVA) to test whether multiple effects significantly impact alpha diversity. This test is present in the `longitudinal` plugin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "faith_pd_anova = longitudinal.visualizers.anova(metadata = sample_metadata.merge(\n",
    "                                                            core_metrics_results.faith_pd_vector.view(qiime2.Metadata)),\n",
    "                                                 formula = 'faith_pd ~ genotype * donor_status')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-dq879tis')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 94a1958c-0ef7-464a-82fc-0b6746df7250>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "faith_pd_anova.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beta diversity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we’ll compare the structure of the microbiome communities using beta diversity. Start by making a visual inspection of the principle coordinates plots (PCoA) plots that were generated by `emperor` and `core_metrics_results/weighted_unifrac_emperor.qzv`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "1. Open the unweighted UniFrac emperor plot (core-metrics-results/unweighted_unifrac_emperor.qzv) first. Can you find separation in the data? If so, can you find a metadata factor that reflects the separation? What if you used weighted UniFrac distance (`core-metrics-results/weighted_unifrac_emperor.qzv`)?\n",
    "\n",
    "2. One of the major concerns in mouse studies is that sometimes differences in communities are due to natural variation in cages. Do you see clustering by cage?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let’s analyze the statistical trends using PERMANOVA. [PERMANOVA](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1442-9993.2001.01070.pp.x) tests the hypothesis that samples within a group are more similar to each other than they are to samples in another group. To put it another way, it tests whether the within-group distances from each group are different from the between group distance. We expect samples that are similar to have smaller distances from each other, so if our hypothesis that one group is different from another is true, we’d expect the within-group distances to be smaller than the between group distance.\n",
    "\n",
    "Let’s use the `beta_group_significance` command to test whether the donor identity (which we qualitatively identified as a major separator in PCoA space) is associated with significant differences in weighted and unweighted UniFrac distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'diversity' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-11594dc19f10>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0munweighted_unifrac_donor_significance\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdiversity\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m                                         \u001b[0mvisualizers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m                                         beta_group_significance(core_metrics_results.\\\n\u001b[1;32m      4\u001b[0m                                                                 \u001b[0munweighted_unifrac_distance_matrix\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m                                                                 sample_metadata.get_column('donor'))\n",
      "\u001b[0;31mNameError\u001b[0m: name 'diversity' is not defined"
     ]
    }
   ],
   "source": [
    "unweighted_unifrac_donor_significance = diversity.\\\n",
    "                                        visualizers.\\\n",
    "                                        beta_group_significance(core_metrics_results.\\\n",
    "                                                                unweighted_unifrac_distance_matrix,\n",
    "                                                                sample_metadata.get_column('donor'))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-_vxeixvf')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: dcaa774b-7d42-4462-afc4-b7832d4ad0b9>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unweighted_unifrac_donor_significance.visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "unexpected character after line continuation character (<ipython-input-3-71ef858d8d3b>, line 3)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-3-71ef858d8d3b>\"\u001b[0;36m, line \u001b[0;32m3\u001b[0m\n\u001b[0;31m    beta_group_significance(core_metrics_results.\\weighted_unifrac_distance_matrix,\u001b[0m\n\u001b[0m                                                                                   ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m unexpected character after line continuation character\n"
     ]
    }
   ],
   "source": [
    "weighted_unifrac_donor_significance = diversity.\\\n",
    "                                        visualizers.\\\n",
    "                                        beta_group_significance(core_metrics_results.\\\n",
    "                                                                weighted_unifrac_distance_matrix,\n",
    "                                                                sample_metadata.get_column('donor'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-5svyyowy')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 2a686e6d-aa4f-4bf4-9ae1-afed86f1f560>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weighted_unifrac_donor_significance.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s also check whether there’s a relationship between the cage in which a mouse lives and the beta diversity, since “cage effect” is often an important technical effect to consider. Since we have several cages, we’ll use the `pairwise` parameter that will let us check whether there are individual differences between the cages driving the difference. This may be useful, since if we check the metadata, we may find that cage is nested by donor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "unweighted_unifrac_cage_significance = diversity.\\\n",
    "                                       visualizers.\\\n",
    "                                       beta_group_significance(core_metrics_results.\\\n",
    "                                                               unweighted_unifrac_distance_matrix,\n",
    "                                                               sample_metadata.get_column('cage_id'),\n",
    "                                                               pairwise = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-mvbjp0gd')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 103f85aa-1fbc-428a-b783-a9163bf04865>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unweighted_unifrac_cage_significance.visualization\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "weighted_unifrac_cage_significance = diversity.\\\n",
    "                                    visualizers.\\\n",
    "                                    beta_group_significance(core_metrics_results.\\\n",
    "                                                            weighted_unifrac_distance_matrix,\n",
    "                                                            sample_metadata.get_column('cage_id'),\n",
    "                                                            pairwise = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-2zx3ukgt')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: fe147b6e-adeb-42ff-ab9f-33c0be157465>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weighted_unifrac_cage_significance.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "1. Is there a significant effect of donor?\n",
    "\n",
    "2. From the metadata, we know that cage C31, C32, and C42 all house mice transplanted from one donor, and that cages C43, C44, and C49 are from the other. Is there a significant difference in the microbial communities between samples collected in cage C31 and C32? How about between C31 and C43? Do the results look the way you expect, based on the boxplots for donor?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A significant difference in PERMANOVA can reflect a large difference between the group or differences in variances within a group. This means that we might see a statistically significant difference even if it’s caused by variation within one group. Distance boxplots can help give a visual sense of this, but it’s nice to use a statistical test to confirm this. We can use the permdisp to help rule out differences due to a high degree of dispersion (within-group variance) in one of the groups of interest.\n",
    "\n",
    "We can specify that we want to use permdisp using the `method` flag in `qiime diversity beta_group_significance`. Let’s explore dispersion based on `cage_id` to check whether are cage-related differences are due to large within-cage variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "unweighted_unifrac_cage_significance_disp = diversity.\\\n",
    "                                            visualizers.\\\n",
    "                                            beta_group_significance(core_metrics_results.weighted_unifrac_distance_matrix,\n",
    "                                                                    sample_metadata.get_column('cage_id'),\n",
    "                                                                    pairwise = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-e8guazlg')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: ce3bf25d-27a4-4248-827b-83e6ae6fd0d2>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unweighted_unifrac_cage_significance_disp.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Is there a significant difference in variance for any of the cages?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the `adonis` action to look at a multivariate model. The adonis action uses a PERMANOVA test, but a different implementation that permits multiple effects to be tested simultaneously (similar to how we used ANOVA earlier for multivariate effects on alpha diversity). Let’s look at the intersection between donor and genotype."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "1. Is there a significant effect of donor?\n",
    "\n",
    "2. From the metadata, we know that cage C31, C32, and C42 all house mice transplanted from one donor, and that cages C43, C44, and C49 are from the other. Is there a significant difference in the microbial communities between samples collected in cage C31 and C32? How about between C31 and C43? Do the results look the way you expect, based on the boxplots for donor?\n",
    "\n",
    "3. If you adjust for donor in the adonis model, do you retain an effect of genotype? What percentage of the variation does genotype explain?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taxonomic classification\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Up until now we have been performing diversity analyses directly on ASVs; in other words, we have assessed the similarity between samples based purely on the unique sequence variants that were observed in each sample. In most experiments we would like to get a sense of what microbial taxa are present — to identify ASVs and give them “names”. To do this, we’ll use the `feature_classifier` plugin to classify ASVs taxonomically.\n",
    "\n",
    "For this analysis, we’ll use a pre-trained naive Bayes machine-learning classifier that was trained to differentiate taxa present in the 99% Greengenes 13_8 reference set trimmed to 250 bp of the V4 hypervariable region (corresponding to the 515F-806R primers). [This classifier works](https://microbiomejournal.biomedcentral.com/articles/10.1186/s40168-018-0470-z) by identifying k-mers that are diagnostic for particular taxonomic groups, and using that information to predict the taxonomic affiliation of each ASV. We can download the pre-trained classifier here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:49:40--  https://data.qiime2.org/2020.2/common/gg-13-8-99-515-806-nb-classifier.qza\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://s3-us-west-2.amazonaws.com/qiime2-data/2020.2/common/gg-13-8-99-515-806-nb-classifier.qza [following]\n",
      "--2020-07-28 18:49:41--  https://s3-us-west-2.amazonaws.com/qiime2-data/2020.2/common/gg-13-8-99-515-806-nb-classifier.qza\n",
      "Resolving s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)... 52.218.176.16\n",
      "Connecting to s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)|52.218.176.16|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 28373581 (27M) [application/x-www-form-urlencoded]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//nb_classifier.qza’\n",
      "\n",
      "/home/user/Document 100%[===================>]  27.06M  1.46MB/s    in 20s     \n",
      "\n",
      "2020-07-28 18:50:02 (1.36 MB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//nb_classifier.qza’ saved [28373581/28373581]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"nb_classifier.qza\" \\\n",
    "  \"https://data.qiime2.org/2020.2/common/gg-13-8-99-515-806-nb-classifier.qza\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "gg_13_8_99_515_806 = qiime2.Artifact.load(workdir+'nb_classifier.qza')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It’s worth noting that Naive Bayes classifiers perform best when they’re trained for the specific hypervariable region amplified. You can train a classifier specific for your dataset based on the [training classifiers tutorial](https://docs.qiime2.org/2020.6/tutorials/feature-classifier/) or download classifiers for other datasets from the [QIIME 2 resource page](https://docs.qiime2.org/2020.6/data-resources/). Classifiers can be re-used for consistent versions of the underlying packages, database, and region of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "taxonomy = feature_classifier.methods.classify_sklearn(denoised_sequences.representative_sequences,\n",
    "                                                       gg_13_8_99_515_806)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let’s review the taxonomy associated with the sequences using the `metadata.visualizers.tabulate`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-ss7uuu_c')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 67bfb016-0fc6-42fc-8dd9-3aed372c6187>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seq_taxonomy = metadata.visualizers.tabulate(taxonomy.classification.view(qiime2.Metadata))\n",
    "seq_taxonomy.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s also tabulate the representative sequences `(FeatureData[Sequence])`. Tabulating the representative sequences will allow us to see the sequence assigned to the identifier and interactively blast the sequence against the NCBI database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "dada2_rep_set_visualization = feature_table.visualizers.tabulate_seqs(denoised_sequences.representative_sequences)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "1. Find the feature, `07f183edd4e4d8aef1dcb2ab24dd7745`. What is the taxonomic classification of this sequence? What’s the confidence for the assignment?\n",
    "\n",
    "2. How many features are classified as `g__Akkermansia`?\n",
    "\n",
    "3 Use the tabulated representative sequences to look up these features. If you blast them against NCBI, do you get the same taxonomic identifier as you obtained with q2-feature-classifier?\n",
    "\n",
    "### Note\n",
    "\n",
    "You might notice that some features do not have taxonomic assignments, which for the Greengenes database is indicated by a blank string at the level (e.g., \"`g__`\"). These indicate that there is not enough information for the Greengenes database to differentiate members of that clade, either due to ambiguity in the database or because the gene region being sequenced doesn’t provide the resolution to distinguish members of that clade. This is distinct from cases where `feature_classifier` cannot reliably classify the ASV to a deeper level: in those cases, an incomplete taxonomy string will be provided. Hence, you may see two different types of “underclassification” in your data: e.g., `k__Bacteria`; `p__Firmicutes; c__Clostridia; o__Clostridiales; f__Christensenellaceae; g__; s__` (genus and species annotations are missing in Greengenes) as well as `k__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Christensenellaceae` (`feature_classifier` could not confidently classify that ASV at genus level).\n",
    "\n",
    "### Note\n",
    "\n",
    "You may also notice that more than one ASV has the same taxonomic assignment. This is normal — unique ASVs do not necessarily map to unique taxonomic groups! We visualize the frequency of each taxonomic group in barplots (as described below) or use the `taxa` plugin to `collapse` our feature table based on taxonomic affiliation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taxonomy barchart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we saw a difference in diversity in this dataset, we may want to look at the taxonomic composition of these samples. To visualize this, we will build a taxonomic barchart of the samples we analyzed in the diversity dataset.\n",
    "\n",
    "Before doing this, we will first filter out any samples with fewer features than our rarefaction threshold (`2000`). We can filter samples using the `feature_table` plugin with the `filter_samples` method. This lets us filter our table based on a variety of criteria such as the number of counts (frequency, `min_frequency` and `max_frequency`), number of features (`min_features` and `max_features`), or sample metadata (`where`). See the [filtering tutorial](https://docs.qiime2.org/2020.6/tutorials/filtering/) for more details and examples.\n",
    "\n",
    "For this example, we need to filter out samples with fewer sequences than our rarefaction depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiime2.plugins.taxa.visualizers import barplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_min_2k = feature_table.methods.filter_samples(denoised_sequences.table,\n",
    "                                                     min_frequency = 2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Now, let’s use the filtered table to build an interactive barplot of the taxonomy in each sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "taxa_barplot_seq = taxa.visualizers.barplot(table_min_2k.filtered_table,\n",
    "                                            taxonomy.classification,\n",
    "                                            sample_metadata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Visualize the data at level 2 (phylum level) and sort the samples by donor, then by genotype. Can you observe a consistent difference in phylum between the donors? Does this surprising you? Why or why not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differential abundance with ANCOM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many microbiome investigators are interested in testing whether individual ASVs or taxa are more or less abundant in different sample groups. This is known as differential abundance. Microbiome data present several challenges for performing differential abundance using convential methods. Microbiome abundance data are inherently sparse (have a lot of zeros) and compositional (everything adds up to 1). Because of this, traditional statistical methods that you may be familiar with, such as ANOVA or t-tests, are not appropriate for performing differential abundance tests of microbiome data and lead to a high false-positive rate. ANCOM is a compositionally aware alternative that allows to test for differentially abundant features. If you’re unfamiliar with the technique, it’s worthwhile to review the [ANCOM paper](https://pubmed.ncbi.nlm.nih.gov/26028277/) to better understand the method.\n",
    "\n",
    "Before we begin, we will filter out low abundance/low prevalence ASVs. Filtering can provide better resolution and limit false discovery rate (FDR) penalty on features that are too far below the noise threshhold to be applicable to a statistical test. A feature that shows up with 10 counts may be a real feature that is present only in that sample, may be a feature that’s present in several samples but only got amplified and sequenced in one sample because PCR is a somewhat stochastic process, or it may be noise. It’s not possible to tell, so feature-based analysis may be better after filtering low abundance features. However, filtering also shifts the composition of a sample, further disrupting the relationship. Here, the filtering is performed as a trade off between the model, computational efficiency, and statistical practicality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "table_2k_abund = feature_table.methods.filter_features(table_min_2k.filtered_table,\n",
    "                                                      min_frequency = 50,\n",
    "                                                      min_samples = 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ANCOM fundamentally operates on a `FeatureTable[Frequency]`, which contains the frequencies of features in each sample. However, ANCOM cannot tolerate zeros (because compositional methods typically use a log-transform or a ratio and you can’t take the log or divide by zeros). To remove the zeros from our table, we add a pseudocount to the `FeatureTable[Frequency]` Artifact, creating a `FeatureTable[Composition]` in its place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "table2k_abund_comp = composition.methods.add_pseudocount(table_2k_abund.filtered_table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s use ANCOM to check whether there is a difference in the mice based on their donor and then by their genetic background. The test will calculate the number of ratios between pairs of ASVs that are significantly different with FDR-corrected p < 0.05."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-1huq71gn')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 300b560a-e33a-46b3-819e-051bd23b4a43>"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ancom_donor = composition.visualizers.ancom(table2k_abund_comp.composition_table,\n",
    "                                            sample_metadata.get_column('donor'))\n",
    "ancom_donor.visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-pvmt6bk8')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 92d2ffd7-7a60-48d2-bd8b-e16543acbb32>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ancom_genotype = composition.visualizers.ancom(table2k_abund_comp.composition_table,\n",
    "                                               sample_metadata.get_column('genotype'))\n",
    "ancom_genotype.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you open the ANCOM visualizations, you’ll see a [volcano plot](https://en.wikipedia.org/wiki/Volcano_plot_(statistics)) on top, which relates the ANCOM W statistic to the CLR (center log transform) for the groups. The W statistic is the number of ANCOM subhypotheses that have passed for each individual taxon, indicating that the ratios of that taxon’s relative abundance to the relative abundances of W other taxa were detected to be significantly different (typically FDR-adjusted p < 0.05). Because differential abundance in ANCOM is based on the ratio between tests, it does not produce a traditional p-value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Open the ANCOM visualizations for the donor and genotype and the taxonomy visualization artifact.\n",
    "\n",
    "1. Are there more differentially abundant features between the donors or the mouse genotype? Did you expect this result based on the beta diversity?\n",
    "\n",
    "2. Are there any features that are differentially abundant in both the donors and by genotype?\n",
    "\n",
    "3. How many differentially abundant features are there between the two genotypes? Using the percentile abundances and volcano plot as a guide, can you tell if they are more abundant in wild type or susceptible mice?\n",
    "\n",
    "4. Use taxonomy metadata visualization and search sequence identifiers for the significantly different features by genotype. What genera do they belong to?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taxonomic classification again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to [increase taxonomic classification accuracy](https://www.biorxiv.org/content/10.1101/406611v2) by showing the taxonomic classifier what a typical animal stool sample looks like before attempting classification. To do that we will have to retrain the naive Bayes classifier. Fortunately, a representation of a typical stool sample that is derived from [Qiita](https://qiita.ucsd.edu/) data is available from the [readytowear collection](https://github.com/BenKaehler/readytowear).\n",
    "\n",
    "If you feel that these samples are not typical stool samples, it is possible to, for instance, assemble data on just mouse or just human (or just human and mouse) stool samples using [q2_clawback](https://library.qiime2.org/plugins/q2-clawback/7/). We will not attempt that here because it takes a while to run, but details are available in the [tutorial](https://forum.qiime2.org/t/using-q2-clawback-to-assemble-taxonomic-weights/5859).\n",
    "\n",
    "Start by downloading the stool data, along with the 99% Greengene 13_8 reference data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:52:28--  https://data.qiime2.org/2020.6/tutorials/pd-mice/ref_seqs_v4.qza\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://qiime2-data.s3-us-west-2.amazonaws.com/2020.6/tutorials/pd-mice/ref_seqs_v4.qza [following]\n",
      "--2020-07-28 18:52:29--  https://qiime2-data.s3-us-west-2.amazonaws.com/2020.6/tutorials/pd-mice/ref_seqs_v4.qza\n",
      "Resolving qiime2-data.s3-us-west-2.amazonaws.com (qiime2-data.s3-us-west-2.amazonaws.com)... 52.218.236.49\n",
      "Connecting to qiime2-data.s3-us-west-2.amazonaws.com (qiime2-data.s3-us-west-2.amazonaws.com)|52.218.236.49|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 9319426 (8.9M) [binary/octet-stream]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//ref_seqs_v4.qza’\n",
      "\n",
      "/home/user/Document 100%[===================>]   8.89M  1.54MB/s    in 6.5s    \n",
      "\n",
      "2020-07-28 18:52:37 (1.36 MB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//ref_seqs_v4.qza’ saved [9319426/9319426]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"ref_seqs_v4.qza\" \\\n",
    "  \"https://data.qiime2.org/2020.6/tutorials/pd-mice/ref_seqs_v4.qza\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:52:37--  https://data.qiime2.org/2020.6/tutorials/pd-mice/ref_tax.qza\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://qiime2-data.s3-us-west-2.amazonaws.com/2020.6/tutorials/pd-mice/ref_tax.qza [following]\n",
      "--2020-07-28 18:52:38--  https://qiime2-data.s3-us-west-2.amazonaws.com/2020.6/tutorials/pd-mice/ref_tax.qza\n",
      "Resolving qiime2-data.s3-us-west-2.amazonaws.com (qiime2-data.s3-us-west-2.amazonaws.com)... 52.218.212.113\n",
      "Connecting to qiime2-data.s3-us-west-2.amazonaws.com (qiime2-data.s3-us-west-2.amazonaws.com)|52.218.212.113|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 2604632 (2.5M) [binary/octet-stream]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//ref_tax.qza’\n",
      "\n",
      "/home/user/Document 100%[===================>]   2.48M  1.33MB/s    in 1.9s    \n",
      "\n",
      "2020-07-28 18:52:41 (1.33 MB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//ref_tax.qza’ saved [2604632/2604632]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"ref_tax.qza\" \\\n",
    "  \"https://data.qiime2.org/2020.6/tutorials/pd-mice/ref_tax.qza\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2020-07-28 18:52:41--  https://data.qiime2.org/2020.6/tutorials/pd-mice/animal_distal_gut.qza\n",
      "Resolving data.qiime2.org (data.qiime2.org)... 52.35.38.247\n",
      "Connecting to data.qiime2.org (data.qiime2.org)|52.35.38.247|:443... connected.\n",
      "HTTP request sent, awaiting response... 302 FOUND\n",
      "Location: https://qiime2-data.s3-us-west-2.amazonaws.com/2020.6/tutorials/pd-mice/animal_distal_gut.qza [following]\n",
      "--2020-07-28 18:52:42--  https://qiime2-data.s3-us-west-2.amazonaws.com/2020.6/tutorials/pd-mice/animal_distal_gut.qza\n",
      "Resolving qiime2-data.s3-us-west-2.amazonaws.com (qiime2-data.s3-us-west-2.amazonaws.com)... 52.218.212.113\n",
      "Connecting to qiime2-data.s3-us-west-2.amazonaws.com (qiime2-data.s3-us-west-2.amazonaws.com)|52.218.212.113|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 229276 (224K) [binary/octet-stream]\n",
      "Saving to: ‘/home/user/Documents/qiime2-mouse-tutorial//animal_distal_gut.qza’\n",
      "\n",
      "/home/user/Document 100%[===================>] 223.90K   257KB/s    in 0.9s    \n",
      "\n",
      "2020-07-28 18:52:43 (257 KB/s) - ‘/home/user/Documents/qiime2-mouse-tutorial//animal_distal_gut.qza’ saved [229276/229276]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget -O $workdir/\"animal_distal_gut.qza\" \\\n",
    "  \"https://data.qiime2.org/2020.6/tutorials/pd-mice/animal_distal_gut.qza\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_seqs_v4 = qiime2.Artifact.load(workdir+'ref_seqs_v4.qza')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "ref_tax = qiime2.Artifact.load(workdir+'ref_tax.qza')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "animal_distal_gut = qiime2.Artifact.load(workdir+'animal_distal_gut.qza')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next retrain the classifier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/user/miniconda2/envs/qiime2-2020.2/lib/python3.6/site-packages/q2_feature_classifier/classifier.py:102: UserWarning: The TaxonomicClassifier artifact that results from this method was trained using scikit-learn version 0.22.1. It cannot be used with other versions of scikit-learn. (While the classifier may complete successfully, the results will be unreliable.)\n",
      "  warnings.warn(warning, UserWarning)\n"
     ]
    }
   ],
   "source": [
    "classifier = feature_classifier.actions.fit_classifier_naive_bayes(ref_seqs_v4,\n",
    "                                                                ref_tax,\n",
    "                                                                animal_distal_gut)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the new classifier in exactly the same way as the standard classifier that we downloaded above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "bespoke_taxonomy = feature_classifier.actions.classify_sklearn(denoised_sequences.representative_sequences,\n",
    "                                                              classifier.classifier)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/user/miniconda2/envs/qiime2-2020.2/lib/python3.6/site-packages/q2_feature_classifier/classifier.py:102: UserWarning: The TaxonomicClassifier artifact that results from this method was trained using scikit-learn version 0.22.1. It cannot be used with other versions of scikit-learn. (While the classifier may complete successfully, the results will be unreliable.)\n",
      "  warnings.warn(warning, UserWarning)\n"
     ]
    }
   ],
   "source": [
    "classifier = feature_classifier.actions.fit_classifier_naive_bayes(ref_seqs_v4,\n",
    "                                                                ref_tax,\n",
    "                                                                animal_distal_gut)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Open up the old `taxonomy` visualization and compare it to the `bespoke_taxonomy` visualization.\n",
    "\n",
    "1. Search for “ovatus” in both. Is there an ASV in the new taxonomy that wasn’t present in the original?\n",
    "\n",
    "2. Revisit the `ancom_donor` visualization. Can you find that ASV?\n",
    "\n",
    "When analyzing ANCOM results, it is possible to trace the ASVs that we found using the taxonomies that we have created. It is also possible to run ANCOM directly on taxonomic groups that we have discovered in our samples by counting features according to taxonomic classification. This has the advantage of pooling feature counts across taxonomically similar ASVs, for instance allowing exact species substitution between samples. The output is also more readable. On the down side, it has all the inaccuracies that come with automated taxonomic classification.\n",
    "\n",
    "We will run through the pipeline twice, once with our original taxonomy and once with the new taxonomy, for the purpose of comparison. First using the original taxonomy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "uniform_table = taxa.actions.collapse(table_min_2k.filtered_table,\n",
    "                                      taxonomy.classification,\n",
    "                                      level = 7)\n",
    "                                      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_uniform_table = feature_table.actions.filter_features(uniform_table.collapsed_table,\n",
    "                                                               min_frequency = 50,\n",
    "                                                               min_samples = 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "cfu_table = composition.actions.add_pseudocount(filtered_uniform_table.filtered_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-xxfqctsw')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 16f48422-e76e-47f9-960b-ce746b0d4cac>"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ancom_donor_uniform = composition.actions.ancom(cfu_table.composition_table,\n",
    "                                               sample_metadata.get_column('donor'))\n",
    "ancom_donor_uniform.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now redo with the new taxonomy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "collapsed_bespoke_table = taxa.actions.collapse(table_min_2k.filtered_table,\n",
    "                                               bespoke_taxonomy.classification,\n",
    "                                                level = 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_bespoke_table = feature_table.actions.filter_features(uniform_table.collapsed_table,\n",
    "                                                               min_frequency = 50,\n",
    "                                                               min_samples = 4)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "cfb_table = composition.actions.add_pseudocount(filtered_bespoke_table.filtered_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-budld4y8')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 759547ca-ee8e-4ca6-8cca-42603ab9f165>"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ancom_donor_bespoke = composition.actions.ancom(cfb_table.composition_table,\n",
    "                                               sample_metadata.get_column('donor'))\n",
    "ancom_donor_bespoke.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Compare final ANCOM visualizations. They are fairly similar, which is good.\n",
    "\n",
    "1. Is Bacteroides ovatus present in the ANCOM results derived from our original taxonomy?\n",
    "\n",
    "2. Is B. ovatus present in the new ANCOM results?\n",
    "\n",
    "3. Why is that?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Longitudinal analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This study includes a longitudinal component; samples from each mouse were collected 7, 14, 21, and 49 days post fecal transplant. We can use the `longitudinal` plug-in to explore the hypothesis that a mouse’s genetic background affected the change in the microbial community of each mouse. For this longitudinal analysis, we’re going to focus on beta diversity. Alpha diversity changes wildly in infants, but it’s often stable in adults over short time periods. We’re dealing with an adult fecal community over a relatively short time period, and there is no difference in alpha diversity with time. The [longitudinal analysis tutorial](https://docs.qiime2.org/2020.6/tutorials/longitudinal/) is an excellent resource for exploring longitudinal analyses of microbiome samples.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCoA-based analyses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can start by exploring temporal change in the PCoA using the animations tab.\n",
    "\n",
    "### Question\n",
    "\n",
    "1. Open the unweighted UniFrac emperor plot and color the samples by mouse id. Click on the “animations” tab and animate using the day_post_transplant as your gradient and mouse_id as your trajectory. Do you observe any clear temporal trends based on the PCoA?\n",
    "\n",
    "2. What happens if you color by day_post_transplant? Do you see a difference based on the day? Hint: Trying changing the colormap to a sequential colormap like viridis.\n",
    "\n",
    "A volatility plot will let us look at patterns of variation along principle coordinate axes starting from the same point. This can be helpful since inter-individual variation can be large and this visualizations lets us focus instead on magnitude of change in each group and in each individual.\n",
    "\n",
    "Let’s use the `longitudinal` plugin to look at how samples from an individual move along each PC. The `sample_metadata` column can take several types, including a metadata file (like our `metadata.tsv`) as well as a `sample_metadata[AlphaDiversity]`, `sample_metadata[Distance]` (files “viewable” as metadata), or a `PCoA` artifact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'longitudinal' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-797e9643bbe3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m pc_vol = longitudinal.actions.volatility(sample_metadata.merge(core_metrics_results.\\\n\u001b[0m\u001b[1;32m      2\u001b[0m                                                                \u001b[0munweighted_unifrac_pcoa_results\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m                                                                view(qiime2.Metadata)),\n\u001b[1;32m      4\u001b[0m                                                                \u001b[0mstate_column\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'days_post_transplant'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m                                                                \u001b[0mindividual_id_column\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'mouse_id'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'longitudinal' is not defined"
     ]
    }
   ],
   "source": [
    "pc_vol = longitudinal.actions.volatility(sample_metadata.merge(core_metrics_results.\\\n",
    "                                                               unweighted_unifrac_pcoa_results.\\\n",
    "                                                               view(qiime2.Metadata)),\n",
    "                                                               state_column = 'days_post_transplant',\n",
    "                                                               individual_id_column = 'mouse_id',\n",
    "                                                               default_group_column = 'donor_status',\n",
    "                                                               default_metric = 'Axis 2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-afc9gtvr')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: c2a9a417-8e76-49bc-b64d-07f19efb4776>"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pc_vol.visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Using the controls, look at variation in cage along PCs 1, 2, and 3. What kind of patterns do you see with time along each axis?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Distance-based analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let’s try looking directly at the pairwise distances between samples. Here, we’ll test the hypothesis that genotype affects the magnitude of change in the distance from the first sample collected from each mouse (7 days post transplant). We assume that given the rate of turnover in a microbial community, we might expect to see a change in the community over time. However, here we’ll ask if these changes are associated with host genotype.\n",
    "\n",
    "We’ll start this analysis by looking at how much the microbial community of each mouse changes from 7 days post transplant. We use the `baseline` parameter to specify a static time point against which all other time points are compared; if we remove this parameter from the command, we look instead at the rate of change for each individual between each time point. See the [longitudinal analysis tutorial](https://docs.qiime2.org/2020.6/tutorials/longitudinal/) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "from_first_unifrac = longitudinal.actions.first_distances(core_metrics_results.\\\n",
    "                                                          unweighted_unifrac_distance_matrix,\n",
    "                                                          metadata = sample_metadata,\n",
    "                                                          state_column = 'days_post_transplant',\n",
    "                                                          individual_id_column = 'mouse_id',\n",
    "                                                          baseline = 7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can again use volatility analysis to visualize the change in beta diversity based on distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'longitudinal' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-8-2618c245d918>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m from_first_unifrac_vol = longitudinal.actions.volatility(sample_metadata.\\\n\u001b[0m\u001b[1;32m      2\u001b[0m                                                          \u001b[0mmerge\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfrom_first_unifrac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfirst_distances\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mview\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mqiime2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMetadata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m                                                          \u001b[0mstate_column\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'days_post_transplant'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m                                                          \u001b[0mindividual_id_column\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'mouse_id'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m                                                          \u001b[0mdefault_metric\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'Distance'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'longitudinal' is not defined"
     ]
    }
   ],
   "source": [
    "from_first_unifrac_vol = longitudinal.actions.volatility(sample_metadata.\\\n",
    "                                                         merge(from_first_unifrac.first_distances.\\\n",
    "                                                               view(qiime2.Metadata)),\n",
    "                                                         state_column = 'days_post_transplant',\n",
    "                                                         individual_id_column = 'mouse_id',\n",
    "                                                         default_metric = 'Distance',\n",
    "                                                         default_group_column = 'donor_status')\n",
    "                                                    \n",
    "from_first_unifrac_vol.visualization                            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "Based on the volatility plot, does one donor change more over time than the other? What about by genotype? Cage?\n",
    "\n",
    "\n",
    "A linear mixed effects (LME) model lets us test whether there’s a relationship between a dependent variable and one or more independent variables in an experiment using repeated measures. Since we’re interested in genotype, we should use this as an independent predictor.\n",
    "\n",
    "For our experiment, we’re currently interested in the change in distance from the initial timepoint, so we’ll use this as our outcome variable (given by `metric`).\n",
    "\n",
    "The `linear_mixed_effects` action also requires a state column (`state_column`) which designates the time component in the metadata, and an individual identifier (`individual_id_column`). Which columns should we use in our data?\n",
    "\n",
    "We can build a model either using the `formula` parameter or the `group_columns` parameter. For this analysis, we’re interested in whether genotype affects the longitudinal change in the microbial community. However, we also know from our cross sectional analysis that donor plays a large role in shaping the fecal community. So, we should also probably include that in this analysis. We may also want to consider cage effect in our experiment, since this is a common confounder in rodent studies. However, the original experimental design here was clever: although cages were grouped by donor (mice are coprophagic), they were of mixed genotype. This partial randomization helps limit some of the cage effects we might otherwise see.\n",
    "\n",
    "Based on the experimental design, what group columns should we choose?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/user/miniconda2/envs/qiime2-2020.2/lib/python3.6/site-packages/statsmodels/base/model.py:568: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n",
      "/home/user/miniconda2/envs/qiime2-2020.2/lib/python3.6/site-packages/statsmodels/regression/mixed_linear_model.py:2114: ConvergenceWarning: Retrying MixedLM optimization with lbfgs\n",
      "  ConvergenceWarning)\n",
      "/home/user/miniconda2/envs/qiime2-2020.2/lib/python3.6/site-packages/statsmodels/regression/mixed_linear_model.py:2149: ConvergenceWarning: The MLE may be on the boundary of the parameter space.\n",
      "  warnings.warn(msg, ConvergenceWarning)\n",
      "/home/user/miniconda2/envs/qiime2-2020.2/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 Distance\" 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/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-hqd35563')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: fe618b88-6ab5-4b5d-b7b8-81e38d287084>"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from_first_unifrac_lme = longitudinal.actions.linear_mixed_effects(sample_metadata.merge(from_first_unifrac.first_distances.view(qiime2.Metadata)),\n",
    "                                                                   state_column = 'days_post_transplant',\n",
    "                                                                   individual_id_column = 'mouse_id',\n",
    "                                                                   metric = 'Distance',\n",
    "                                                                   group_columns = 'genotype,donor')\n",
    "\n",
    "from_first_unifrac_lme.visualization                "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let’s look at the results of the models\n",
    "\n",
    "### Question\n",
    "\n",
    "1. Is there a significant association between the genotype and temporal change?\n",
    "\n",
    "2. Which genotype is more stable (has lower variation)?\n",
    "\n",
    "3. Is there a temporal change associated with the donor? Did you expect or not expect this based on the volatility plot results?\n",
    "\n",
    "4. Can you find an interaction between the donor and genotype?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note\n",
    "\n",
    "Importantly, LME models also allow us to distinguish between two types of independent variables: fixed effects (e.g., experimental treatments) and random effects (e.g., biological factors that cannot be controlled in the experiment). By default, the `linear_mixed_effects` action in `longitudinal` uses the `individual_id_column` as a random effect, since we can expect that biological differences between individual subjects will impact the baseline values of the dependent variable we are testing. The rate of change — slope — is another inter-individual effect that we often might want to consider as a random effect in longitudinal experiments. See the [longitudinal analysis](https://docs.qiime2.org/2020.6/tutorials/longitudinal/) tutorial for more details and discussion of LME tests and effect types."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine-learning classifiers for predicting sample characteristics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative (or complementary) approach to the methods we have used in this tutorial for testing if and how samples are different from one another, we can utilize [machine-learning methods](https://docs.qiime2.org/2020.6/tutorials/sample-classifier/) to determine how predictive microbiome composition is of other characteristics about a sample. For example, we may use machine-learning classifiers to predict a patient’s susceptibility to disease, or predict the treatment group that a sample belongs to. Additionally, many machine-learning methods report which features are most important for predicting sample characteristics, making this a useful approach for determining which features (ASVs, species, etc) are associated with a particular treatment, disease state, or other category of interest. All of this and much more can be found in the `sample_classifier` plugin. Here we will use this plugin to predict each mouse’s genotype and donor status based on their ASV composition using a Random Forest classifier (this pipeline can access many different machine-learning methods via the estimator parameter, but Random Forest classifiers are used by default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'sample_classifier' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-3dcca86f6dd9>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m sample_classifier_results = sample_classifier.actions.classify_samples(denoised_sequences.table,\n\u001b[0m\u001b[1;32m      2\u001b[0m                                                                        \u001b[0msample_metadata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m                                                                        \u001b[0mget_column\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'genotype_and_donor_status'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m                                                                        \u001b[0mrandom_state\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m666\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m                                                                        n_jobs = 1)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'sample_classifier' is not defined"
     ]
    }
   ],
   "source": [
    "sample_classifier_results = sample_classifier.actions.classify_samples(denoised_sequences.table,\n",
    "                                                                       sample_metadata.\\\n",
    "                                                                       get_column('genotype_and_donor_status'),\n",
    "                                                                       random_state = 666,\n",
    "                                                                       n_jobs = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><img onload=\"(function(div, url){\n",
       "if (typeof require !== 'undefined') {\n",
       "    var baseURL = require.toUrl('').split('/').slice(0, -2).join('/');\n",
       "} else {\n",
       "    var baseURL = JSON.parse(\n",
       "        document.getElementById('jupyter-config-data').innerHTML\n",
       "    ).baseUrl.slice(0, -1);\n",
       "}\n",
       "url = baseURL + url;\n",
       "fetch(url).then(function(res) {\n",
       "    if (res.status === 404) {\n",
       "        div.innerHTML = 'Install QIIME 2 Jupyter extension with:<br />' +\n",
       "                        '<code>jupyter serverextension enable --py qiime2' +\n",
       "                        ' --sys-prefix</code><br />then restart your server.' +\n",
       "                        '<br /><br />(Interactive output not available on ' +\n",
       "                        'static notebook viewer services like nbviewer.)';\n",
       "    } else if (res.status === 409) {\n",
       "        div.innerHTML = 'Visualization no longer in scope. Re-run this cell' +\n",
       "                        ' to see the visualization.';\n",
       "    } else if (res.ok) {\n",
       "        url = res.url;\n",
       "        div.innerHTML = '<iframe src=\\'' + url + '\\' style=\\'' +\n",
       "                        'width: 100%; height: 700px; border: 0;\\'>' +\n",
       "                        '</iframe><hr />Open in a: <a href=\\'' + url + '\\'' +\n",
       "                        ' target=\\'_blank\\'>new window</a>'\n",
       "    } else {\n",
       "        div.innerHTML = 'Something has gone wrong. Check notebook server for' +\n",
       "                        ' errors.';\n",
       "    }\n",
       "});\n",
       "})(this.parentElement, '/qiime2/redirect?location=/tmp/qiime2-archive-v2oetses')\" src=\"data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==\" /></div>"
      ],
      "text/plain": [
       "<visualization: Visualization uuid: 43b6ca07-6850-45d7-a0b4-545933b1222d>"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sample_classifier_results.accuracy_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This pipeline generates a number of output artifacts and visualizations. You can read more about these in the [sample classifier](https://docs.qiime2.org/2020.6/tutorials/sample-classifier/) tutorial but right now let’s just focus on `sample_classifier_results` This visualization tells you how well your sample classifier performed via a [confusion matrix](https://en.wikipedia.org/wiki/Confusion_matrix) and accompanying table of accuracy scores. This tells you how frequently each sample type is classified to each sample class, including the correct class label. Overall error rates are also reported in the table below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "How did we do? Just for fun, try predicting some of the other metadata columns to see how easily `cage_id` and other columns can be predicted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like we did pretty well! So we can see what features are most predictive of each sample class (donor and genotype groups). The importance scores are stored in the ./sample-classifier-results/feature_importance.qza artifact (pro tip: this can be view with the qiime metadata tabulate command we covered earlier). Here we will generate a heatmap showing the mean abundance of the 100 most important ASVs in each genotype and donor group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_table_heatmap = sample_classifier.actions.heatmap(denoised_sequences.table,\n",
    "                                                           sample_classifier_results.feature_importance,\n",
    "                                                           sample_metadata.get_column('genotype_and_donor_status'),\n",
    "                                                           group_samples = True,\n",
    "                                                           feature_count = 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "\n",
    "What features appear to differentiate genotypes? What about donors? Are any ASVs specific to a single sample group?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synthesis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the results of the analysis, we can say that there is a difference in the microbial communities of these mice based on their donor and genetic background. (This recapitulates the results of the original analysis.)\n",
    "\n",
    "We found that the donor is the primary driver of alpha diversity.\n",
    "\n",
    "But, we saw differences by donor and genotype based on beta diversity. Using the PCoA emperor plots, we can see clear separation between the mice from the two donors. After adjusting for the donor, we saw a significant difference between the genotypes.\n",
    "\n",
    "Although there wasn’t a clear pattern in the barchart at the phylum level between donors or genotypes, we were still able to find ASVs which differentiated the genotypes using ANCOM and Random Forest classification. There was no overlap between these ASVs in the donor and genetic background, supporting the hypothesis that the difference due to genotype is separate from the difference due to donor.\n",
    "\n",
    "The volatility plots and temporal analysis showed that the microbiome in different genetic backgrounds changed differently over time.\n",
    "\n",
    "This suggests that there is a genotype-specific effect on the microbiome of mice receiving fecal transplants.\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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}