16S_2024.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/Gibbons-Lab/isb_course_2024/blob/main/16S_2024.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "2inlS7SNn0-o",
    "tags": []
   },
   "source": [
    "# 🦠 Amplicon Sequencing Data Analysis with QIIME 2\n",
    "\n",
    "### 🔗 Important Links\n",
    "- [Course Github](https://github.com/Gibbons-Lab/isb_course_2024): Github repository for the 2024 ISB Virtual Microbiome Series\n",
    "\n",
    "    - [Day 1 Slides](https://gibbons-lab.github.io/isb_course_2024/16S): Presentation slides that accompany this notebook\n",
    "\n",
    "    - [Day 1 Solutions](https://colab.research.google.com/github/Gibbons-Lab/isb_course_2024/blob/main/16S_2024_solutions.ipynb): Pre-run notebook \n",
    "\n",
    "    - [Treasure chest](https://github.com/Gibbons-Lab/isb_course_2024/tree/main/treasure_chest): All files generated in the notebooks\n",
    "    \n",
    "- [QIIME 2 view](http://view.qiime2.org): A browser-based interface for viewing QIIME artifacts and visualizations (tables and plots)\n",
    "- [QIIME 2 plugins](https://docs.qiime2.org/2024.5/plugins/): more information on the QIIME 2 plugins we're using today"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lcJVyVD80fAG"
   },
   "source": [
    "### 🛑 STOP! Before you run anything ...\n",
    "\n",
    "1. Save your own local copy of this notebook by using `File > Save a copy in Drive`. At some point you may be prompted to trust the notebook. We promise that it is safe 🤞\n",
    "\n",
    "2. Keep in mind, **the code in this notebook must be run IN ORDER** If you get lost or encounter errors, all outputs can be found in `materials/treasure chest` or on the link above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Disclaimer:** The Google Colab notebook environment will interpret any command as Python code by default. If we want to run bash commands we will have to prefix them by `!`. So any command you see with a leading `!` is a bash command and if you wanted to run it in your terminal you would omit the `!`. For example, if in the Colab notebook you ran `!wget` you would just run `wget` in your terminal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UVawr3gG0fAG",
    "tags": []
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CuWOn36b0fAH"
   },
   "source": [
    "QIIME 2 is usually installed by following the [official installation instructions](https://docs.qiime2.org/2024.5/install/). However, because we are using Google Colab and there are some caveats to using conda here, we will have to hack around the installation a little bit. But no worries, we provide a setup script below which does all this work for us. 😌\n",
    "\n",
    "So...let's start by pulling a local copy of the project repository down from GitHub. This repository, called `materials`, contains all the relevant data and other resources we'll need for this course. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "9Cyg41Tcn1Ua",
    "outputId": "a207fd86-51e3-4e07-8ae9-a6cbafc4ef18"
   },
   "outputs": [],
   "source": [
    "!git clone https://github.com/gibbons-lab/isb_course_2024 materials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MQ9XVvOToEBR"
   },
   "source": [
    "To view the directory, you can click on the folder 📁 icon on the left. To run code from this directory, we will navigate to it via command line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "GvNHs-EroFNM",
    "outputId": "ce246af8-4475-47c8-d71f-7eb3b431c914"
   },
   "outputs": [],
   "source": [
    "%cd materials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IHvKnTbdypyp"
   },
   "source": [
    "Notice here we use ```%``` instead of ```!``` to run our command line function. This makes the path change to our directory _permanent_. Using the ```!``` operator only switches the interpreter to expect command line prompts temporarily.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nJ56MqA50fAJ"
   },
   "source": [
    "Now that we have all our materials, we're _almost_ ready to get started - but not quite. Remember QIIME? We'll need to install that before getting into the actual analysis. Don't worry - this will only set up in the Colab notebook, not on your local machine.\n",
    "\n",
    "**Let's run the following cell, to install and setup QIIME2**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 298
    },
    "id": "Inkvkt_DoSQT",
    "outputId": "278c4207-ca12-4ecf-d562-2335990d6170"
   },
   "outputs": [],
   "source": [
    "%run setup_qiime2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Pajd6O0OpE_f"
   },
   "source": [
    "⬆️ This will take some time (usually 10 to 15 minutes), so we'll switch back over to the [presentation](https://gibbons-lab.github.io/isb_course_2024/16S) while we wait.\n",
    "\n",
    "If you want to learn more about QIIME2, we recommend you check out the [documentation](https://docs.qiime2.org/). This will also explain how to install QIIME2 on your local machine 🖥"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RccdNYBeqQVn",
    "tags": []
   },
   "source": [
    "---\n",
    "# Data Import"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QIIME2 can take us from raw sequences to ecological insight!\n",
    "![our workflow](https://github.com/Gibbons-Lab/isb_course_2024/raw/main/docs/16S/assets/steps.png)\n",
    "But we first need to get our raw data into QIIME2. To do this, we will need:\n",
    "1. **Raw sequencing files** (.fastq): FASTQ files contain both sequence and quality scores, as opposed to FASTA, which only have sequence. We want those quality scores so that we can trim or discard low-quality reads.\n",
    "\n",
    "2. **Metadata file** (.tsv): This file contains the sample id, along with all non-sequencing data for each sample. Here, our metadata includes disease status, age, sex, BMI, geographical location, and prescription drug use, as well as some information on diet and GI symptoms.\n",
    "\n",
    "3. **Manifest file** (.tsv): This file contains the sample id and the path to its sequencing file. This is how QIIME matches a sequence to its correct metadata. 📝"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bPExCVcl0fAK"
   },
   "source": [
    "Let's begin by taking a look at our data. In the `data` folder, you'll find 10 FASTQ files, a file manifest, and a metadata file.\n",
    "\n",
    "If we read in `manifest.tsv`, we can see that it is a table containing the _name_ and _filepath_ of all our samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 362
    },
    "id": "_pFaCnO6iqqW",
    "outputId": "160bf3f3-5f7c-4460-ef75-9599327044f8"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "manifest = pd.read_csv('data/manifest.tsv', sep = '\\t')\n",
    "manifest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vAIoc31bgTY_"
   },
   "source": [
    "We can also check out `metadata.tsv`, which will give us more context on our samples 🔬"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 712
    },
    "id": "sAfnoRnxgZOO",
    "outputId": "25ad5351-711f-41a9-82e8-fa6acec4406e",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "metadata = pd.read_csv('data/metadata.tsv', sep='\\t')\n",
    "metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_ssDjv0I0fAL"
   },
   "source": [
    "Looks good, all 10 FASTQ files are accounted for, five healthy and five with Parkinson's Disease. We can use the manifest file to import our files into QIIME."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qc10jMUYsXyh",
    "tags": []
   },
   "source": [
    "## FASTQ to Artifact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HLF-zdrJXq7q"
   },
   "source": [
    "To use sequencing data in QIIME, we first need to turn the FASTQ files containing our data into QIIME artifacts. Using the manifest we just checked out, let's run our first command:\n",
    "\n",
    "-- as a reminder, adding ```!``` before the command tells the notebook this is a **bash** command, rather than python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "QmLebzqFX_0N",
    "outputId": "022d22f5-5a0e-4b91-c125-2cedc826f03e"
   },
   "outputs": [],
   "source": [
    "!qiime tools import \\\n",
    "  --type 'SampleData[SequencesWithQuality]' \\\n",
    "  --input-path data/manifest.tsv \\\n",
    "  --output-path sequences.qza \\\n",
    "  --input-format SingleEndFastqManifestPhred33V2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "N2Pmtg-aaqvr"
   },
   "source": [
    "Let's take a look a the command. QIIME commands take following format:\n",
    "\n",
    "```\n",
    "qiime plugin action --i-argument1 ... --o-argument2 ...\n",
    "```\n",
    "In the previous command, we are calling the ```tools``` plugin within QIIME2 to import our data. The following arguments designate where the manifest is, as well as where the output should be saved. We also tell QIIME what sort of input to expect.\n",
    "\n",
    "Argument types usually begin with a letter denoting their meaning:\n",
    "\n",
    "- `--i-...` = input files\n",
    "- `--o-...` = output files\n",
    "- `--p-...` = parameters\n",
    "- `--m-...` = metadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Wb__5zxD0fAM",
    "tags": []
   },
   "source": [
    "## Artifact to Visualization 🔎\n",
    "\n",
    "Before we move on, let's use QIIME to visualize our sequencing data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "nVPsDHTNvNUw",
    "outputId": "eb7034d2-910b-4825-f668-18ff83978859"
   },
   "outputs": [],
   "source": [
    "!qiime demux summarize \\\n",
    "--i-data sequences.qza \\\n",
    "--o-visualization qualities.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "o29TMwWOy9L1"
   },
   "source": [
    ".qzv files like the one we just produced are visualization. You can view the plot by **downloading the file and opening it using http://view.qiime2.org**. To download the file click on the folder symbol to the left, open the `materials` folder, and choose download from the dot menu next to the `qualities.qzv` file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "va8RVKqP0fAM",
    "tags": []
   },
   "source": [
    "---\n",
    "# Quality Filtering: From Sequence to ASV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Lxb_lE8l0fAN"
   },
   "source": [
    "\n",
    "Before we can use our sequencing data, we need to \"denoise\" it. To do this, we'll use a *plugin* called ```DADA2```. This involves a few things:\n",
    "\n",
    "1. filter and trim the reads\n",
    "2. find the most likely set of unique sequences in the sample (ASVs)\n",
    "3. remove chimeras\n",
    "4. count the abundances of each ASV\n",
    "\n",
    "\n",
    "This command will take a little time - let's run it, and head back to the presentation to discuss what's happening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "HclmzHqyacLU",
    "outputId": "537ecb15-7847-4aa8-f38c-f49a75b8d85b"
   },
   "outputs": [],
   "source": [
    "!qiime dada2 denoise-single \\\n",
    "    --i-demultiplexed-seqs sequences.qza \\\n",
    "    --p-trunc-len 150 \\\n",
    "    --p-n-threads 2 \\\n",
    "    --output-dir dada --verbose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sEs5k1Za-54R"
   },
   "source": [
    "\n",
    "If this step takes too long or fails, you can also copy the results from the treasure chest with the following command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-Sba9UTNeoji"
   },
   "outputs": [],
   "source": [
    "# obscure magic that will only copy if the previous command failed\n",
    "![ -d dada ] || cp -r treasure_chest/dada ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we open up the `dada` folder 📁, we will see several outputs.\n",
    "- `representative_sequences.qza` is our *ASV abundance table* 🧬\n",
    "- `denoising-stats.qza` contains the *denoising statistics* 📊\n",
    "\n",
    "We'll want to look at the denoising statistics to make sure the denoising went well. It is normal to lose between 5-25% of reads to chimeras, but major loss could be due to running too many PCR cycles.\n",
    "\n",
    "As mentioned before, we can't view QIIME *artifacts* directly, but we can convert them to *visualizations* using the action `tabulate` from the plugin `metadata`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "z-iW0UBd_Iiu",
    "outputId": "c0c41e43-5e54-4290-89e3-9e12a27701fa"
   },
   "outputs": [],
   "source": [
    "!qiime metadata tabulate \\\n",
    "    --m-input-file dada/denoising_stats.qza \\\n",
    "    --o-visualization dada/denoising-stats.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TWqEgfzRKA5t"
   },
   "source": [
    "Like before, we can download the .qzv file and visualize the results using the [QIIME2 Viewer](https://view.qiime2.org/).\n",
    "\n",
    "It's important to understand what this output tells us. For instance, what percent of reads in our data pass the filtering step? What percent of reads were non-chimeric? Differences in these metrics between samples can affect diversity metrics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information on paired-end quality filtering, check out the [2023 course](https://github.com/Gibbons-Lab/isb_course_2023)! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mQKTRQXJ0fAO",
    "tags": []
   },
   "source": [
    "---\n",
    "# Phylogenetics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have an ASV table, we can calculate diversity metrics! But first, head back to the [presentation](https://gibbons-lab.github.io/isb_course_2024/16S) for an introduction to diversity!\n",
    "\n",
    "Note: we're getting into analyses that some may prefer to perform in RStudio or Jupyter Notebooks, without the use of the QIIME framework. For example, diversity analyses can be performed using `scikit-bio` in Python or `vegan` in R. In fact, QIIME2 is using `scikit-bio` and `vegan` under the hood for its diversity calculations!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bPoyj-6c0fAU",
    "tags": []
   },
   "source": [
    "### From ASVs to Phylogenetic Trees with QIIME2 `phylogeny`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Bw8D0aML0fAU"
   },
   "source": [
    "Some *diversity metrics* are calculated using **phylogenetic distance**, so let's start by building a phylogenetic tree for our sequences using the following command. This time, we call the `phylogeny` plugin in QIIME2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Ym7blrPpKQuJ",
    "outputId": "c8b24f51-63af-4f52-b4cb-d1fef9d72a00"
   },
   "outputs": [],
   "source": [
    "!qiime phylogeny align-to-tree-mafft-fasttree \\\n",
    "    --i-sequences dada/representative_sequences.qza \\\n",
    "    --output-dir tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "T1jbqaW1hrmQ"
   },
   "source": [
    "We can create a visualization for the tree using the `empress` QIIME 2 plugin. See [documentation](https://github.com/biocore/empress) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "2Ttf0RTyhaPH",
    "outputId": "093a0e79-c6e1-4bfb-9804-3676d0a87cce"
   },
   "outputs": [],
   "source": [
    "!qiime empress tree-plot \\\n",
    "    --i-tree tree/rooted_tree.qza \\\n",
    "    --o-visualization tree/empress.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Diversity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H4fTVpeToufK",
    "tags": []
   },
   "source": [
    "## Part 1: The Basics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dyp1iAoW0fAV"
   },
   "source": [
    "The QIIME2 `diversity` plugin has an action called `core-metrics-phylogenetic`, which calculates _several_ phylogenetic and non-phylogenetic diversity metrics.\n",
    "\n",
    "We can use our table and tree to calculate several diversity metrics. To account for variations in **sampling depth**, we'll provide QIIME2 with a cutoff at which we will **rarefy** all our samples. Rarefaction  **randomly** selects sequences, so your results might look a little different from mine. We'll also pass in our metadata file, so we can keep track of which samples come from each group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "1YqaTqS5huAz",
    "outputId": "2c082204-f04f-406d-b6a2-f4f401fef260"
   },
   "outputs": [],
   "source": [
    "!qiime diversity core-metrics-phylogenetic \\\n",
    "    --i-table dada/table.qza \\\n",
    "    --i-phylogeny tree/rooted_tree.qza \\\n",
    "    --p-sampling-depth 5000 \\\n",
    "    --m-metadata-file data/metadata.tsv \\\n",
    "    --output-dir diversity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7kiYQC_30fAW"
   },
   "source": [
    "If you open the `diversity` folder, you'll see that we calculated all of the diversity metrics I just mentioned. For more information on how these diversity metrics are calculated and interpreted, check out this QIIME [forum post](https://forum.qiime2.org/t/alpha-and-beta-diversity-explanations-and-commands/2282).\n",
    "\n",
    "But what we're really interested in is whether diversity differs *significantly* between disease and healthy samples, or other covariates. Therefore, we will need to do statistical tests. \n",
    "\n",
    "Head back to the [slides](https://gibbons-lab.github.io/isb_course_2024/16S), and we'll go over some types of statistical tests, and when to use them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XslTNgFBg5qp",
    "tags": []
   },
   "source": [
    "### Alpha diversity: richness and evenness _within_ a sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AVT3NfDT0fAW"
   },
   "source": [
    "We get a bunch of outputs from the previous command - measures of both alpha and beta diversity. To start, let's use the Shannon vector in the output directory to create a visualization of alpha diversity across samples. Generally, healthy, long-living individuals have balanced diverse microbiomes. However, this isn't necessarily a direct indicator of health or disease. Let's see how it looks in our samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "ne2WybKgh6fv",
    "outputId": "c4f2a54d-83af-4c37-c51e-e78c652dbaf2"
   },
   "outputs": [],
   "source": [
    "!qiime diversity alpha-group-significance \\\n",
    "    --i-alpha-diversity diversity/shannon_vector.qza \\\n",
    "    --m-metadata-file data/metadata.tsv \\\n",
    "    --o-visualization diversity/alpha_groups.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HzF0iOjD0fAW"
   },
   "source": [
    "Like before, we can download the visualization and open it with the [QIIME2 viewer](http://view.qiime2.org)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "V5eJbSNo0fAX",
    "tags": []
   },
   "source": [
    "### Beta diversity: distance between samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VNmxg2lmjyHE"
   },
   "source": [
    "Let's visualize the beta diversity and see how they separate. For this we'll look at weighted UniFrac. This time, we'll have to download the file ⬅️\n",
    "\n",
    "We can check for 'significant' separation between samples using PERMANOVA. We can do this with the diversity plugin in QIIME2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "WMASNLaWjxwD",
    "outputId": "969fdd0d-b65c-4351-8982-86c2972d37b7"
   },
   "outputs": [],
   "source": [
    "!qiime diversity adonis \\\n",
    "    --i-distance-matrix diversity/weighted_unifrac_distance_matrix.qza \\\n",
    "    --m-metadata-file data/metadata.tsv \\\n",
    "    --p-formula \"parkinson_disease\" \\\n",
    "    --p-n-jobs 2 \\\n",
    "    --o-visualization diversity/permanova.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YtogNk5h0fAY"
   },
   "source": [
    "We already ran a PCoAs for each distance metric, and we can look at them if we open `weighted_unifrac_emperor.qzv`, `unweighted_unifrac_emperor.qzv`, `bray_curtis_emperor.qzv`, or `jaccard_emperor.qzv` in [QIIME2 viewer](http://view.qiime2.org)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Are these results informative? A discussion of sample size, effect size, and confounding"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results aren't very informative, and in some cases they are counterintuitive. There are several reasons this may be:\n",
    "\n",
    "1. **Small sample size:** We are using a very small sample size (n = 10). Even if we could observe trends, we wouldn't expect statistical significance.\n",
    "\n",
    "2. **Alpha diversity is a very coarse-grain metric:** Shannon Index attempts to summarize microbiome diversity (very complicated) into a single number.\n",
    "\n",
    "3. **Strong confounders:** This data is from an observational study, were confounders were not controlled. Some confounders (like diet and drug exposure) are known to have strong effects on diversity, which may conceal smaller relationships between the microme and a variable of interest (like disease status).\n",
    "\n",
    "In observational studies, it is common to **stratify** across confounders, or include them as **covariates** in regressions. But with a sample size n = 10, we can really only test 1 variable without breaking assumptions. \n",
    "\n",
    "**But what if we had a bigger sample size?** This would open the door to more types of analyses. We couldn't process the entire dataset on Colab due to limited computational power, but luckily, **we've pre-processed the entire dataset for you**!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Dealing with confounders"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Alpha Diversity - complete dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full Hill-Burns 2017 dataset files are in the MATERIALS FOLDER. Let's load them into QIIME2 now.\n",
    "\n",
    "We'll keep this analysis in the `complete_study` folder so we don't confuse our files from earlier.\n",
    "\n",
    "If you get lost at any point, the outputs generated from the complete dataset can be found in `treasure_chest/complete_study_outputs`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy full dataset directory into your local directory\n",
    "![ -d complete_study ] || cp -r treasure_chest/complete_study ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a LOT more metadata than before. Let's take a peek."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "complete_metadata = pd.read_csv('complete_study/complete_metadata.tsv', sep='\\t')\n",
    "print(complete_metadata.shape)\n",
    "complete_metadata.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the original paper, the authors found that anticholinergic medications, catechol-O-methyltransferase (comt) inhibitors, location, sex, age, and PD had significant, independent signals on microbiome composition. How does alpha diversity vary across these variables?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!qiime diversity alpha-group-significance \\\n",
    "    --i-alpha-diversity complete_study/diversity/shannon_vector.qza \\\n",
    "    --m-metadata-file complete_study/complete_metadata.tsv \\\n",
    "    --o-visualization complete_study/diversity/alpha_groups.qzv"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, view in [QIIME2 viewer](http://view.qiime2.org).\n",
    "\n",
    "Notice that `alpha-group-significance` performs a _separate_ Kruskal-Wallis test for each variable, and outputs an _un-adjusted p-value_. If we want to analyze _multiple_ alpha diversity-covariate relationships (i.e. perform multiple tests), then we need to _control the false discovery rate (FDR) by correcting our p-values_. We will go over this later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### PCoA + PERMANOVA for Beta diversity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we now have n = 300+, we can include more variables in our PERMANOVA."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's run PERMANOVA for just disease status so we can really see what adding covariates does to our results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!qiime diversity adonis \\\n",
    "    --i-distance-matrix complete_study/diversity/weighted_unifrac_distance_matrix.qza \\\n",
    "    --m-metadata-file complete_study/complete_metadata.tsv \\\n",
    "    --p-formula \"parkinson_disease\" \\\n",
    "    --p-n-jobs 2 \\\n",
    "    --o-visualization complete_study/diversity/permanova_PD.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you look at the table in [QIIME2 viewer](http://view.qiime2.org), you'll notice that since our sample size increased, the effect of PD on UniFrac beta diversity (though very small at R-squared = 0.01) is now significant (at least before FDR corrections)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use PERMANOVA to identify confounders. PERMANOVA tells us how much variance in the community composition is explained by each variable. Common confounders include sex, age, BMI, diet, and antibiotic use. In the original study, authors identified that Parkinson's Disease medications were associated with different microbiome compositions. Let's take a look at these variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But wait - `adonis` requires any variables in the formula to have no NaN/missing values! If we tried to run `adonis` with all of these variables, it would throw an error. So we will have to **filter** out samples that have NA values for any variables we want to test. \n",
    "\n",
    "Let's look at the metadata file to see which variables have NaN values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata_perm = complete_metadata[['parkinson_disease', 'sex', 'age', 'bmi', 'comt_inhibitor', 'carbidopa_levodopa', 'anticholinergic', 'location']].copy()\n",
    "metadata_perm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what variables have N/A values for the PERMANOVA we want to run?\n",
    "col_na_num = metadata_perm.isna().sum()\n",
    "col_na = col_na_num[col_na_num > 0]\n",
    "col_na"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like a subset of people did not fill out these questions. How many samples will we have if we filter them out?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata_perm.dropna(axis=0, subset=['bmi', 'comt_inhibitor', 'carbidopa_levodopa', 'anticholinergic'], how='any')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the QIIME action `filter-distance-matrix` this time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter QIIME table\n",
    "!qiime diversity filter-distance-matrix \\\n",
    "    --i-distance-matrix complete_study/diversity/weighted_unifrac_distance_matrix.qza \\\n",
    "    --m-metadata-file complete_study/complete_metadata.tsv \\\n",
    "    --p-where \"[bmi] IS NOT NULL AND [comt_inhibitor] IS NOT NULL AND [carbidopa_levodopa] IS NOT NULL AND [anticholinergic] IS NOT NULL\" \\\n",
    "    --o-filtered-distance-matrix complete_study/diversity/weighted_unifrac_distance_matrix_no_NA.qza"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!qiime diversity adonis \\\n",
    "    --i-distance-matrix complete_study/diversity/weighted_unifrac_distance_matrix_no_NA.qza \\\n",
    "    --m-metadata-file complete_study/complete_metadata.tsv \\\n",
    "    --p-formula \"parkinson_disease + sex + age + bmi + location + comt_inhibitor + carbidopa_levodopa + anticholinergic\" \\\n",
    "    --p-n-jobs 2 \\\n",
    "    --o-visualization complete_study/diversity/permanova_with_confounders.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before, we did not see a significant p-value for the effect of Parkinson's disease on beta diversity. However, when we add certain covariates, we might find that they were confounding a relationship. \n",
    "\n",
    "However, most of our variance remains unexplained. Microbiome composition is affected by many things, and this is an uncontrolled cohort study, so we would not expect any single variable to explain most of the variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dimensionality reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to __visually__ show this separation between samples, we can't just plot the entire UniFrac distance matrix, because it has 100+ dimensions! Instead, we can use **dimensionality reduction** to \"compress\" our data into a few dimensions that explain most of the variance. There are several types of dimensionality reduction (like UMAP and tSNE), but the preferred method of dimensionality reduction for microbiome communities is Principal Coordinate Analysis (PCoA). This is because PCoA is linear, and thus preserves the global structure of the data and is reproducible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We already ran a PCoAs for each distance metric, and we can look at them if we download `weighted_unifrac_emperor.qzv`, `unweighted_unifrac_emperor.qzv`, `bray_curtis_emperor.qzv`, or `jaccard_emperor.qzv`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XTXoJhjs0fAY"
   },
   "source": [
    "---\n",
    "# Taxonomic Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### From ASV to taxonomic count table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "CR2f0q4trGfA"
   },
   "source": [
    "We can learn a lot from diversity metrics, alpha and beta. But to really dig into the data, we need to know what microbes are in each sample 🦠. To do this, we'll classify the reads in QIIME2 using a pre-trained Bayesian classifier. Several such classifiers are available at https://docs.qiime2.org/2024.5/data-resources/\n",
    "\n",
    "From the QIIME2 plugin `feature-classifier`, we will use the action the action `classifiy-sklearn` which takes your ASVs and a pre-trained scikit-learn classifier as inputs. We are using a Naive Bayesian classifier specific to the 16S V4 region (515f-806r) because that is the region we sequence, but also it will make it easier to run MICOM tomorrow!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "firXAVIHrEzT",
    "outputId": "44fd0067-ad26-4083-ea0b-2945c9efc345"
   },
   "outputs": [],
   "source": [
    "# run a classifier on our 10-sample mini dataset\n",
    "!qiime feature-classifier classify-sklearn \\\n",
    "    --i-reads dada/representative_sequences.qza \\\n",
    "    --i-classifier ncbi-refseq-genus-515f-806r.qza \\\n",
    "    --p-n-jobs 2 \\\n",
    "    --o-classification taxa.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ftEKkIx-uJtA"
   },
   "source": [
    "Now we've classified the reads, we can visualize the taxonomic breakdown of our samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "RSRrxhX7x4dE",
    "outputId": "09dbb1eb-8b03-453c-d3ed-1aa0642a7574"
   },
   "outputs": [],
   "source": [
    "!qiime taxa barplot \\\n",
    "    --i-table dada/table.qza \\\n",
    "    --i-taxonomy taxa.qza \\\n",
    "    --m-metadata-file data/metadata.tsv \\\n",
    "    --o-visualization taxa_barplot.qzv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LbZ2f-kOx4jY"
   },
   "source": [
    "Now, we can use ```table.qza```, which contains our reads, and ```taxa.qza```, which contains taxonomic classifications for reads, and collapse the data onto the genus level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "XTClvRnXx4ob",
    "outputId": "ae841b68-ecc2-4e0b-cd23-7a2e415bd206"
   },
   "outputs": [],
   "source": [
    "!qiime taxa collapse \\\n",
    "    --i-table dada/table.qza \\\n",
    "    --i-taxonomy taxa.qza \\\n",
    "    --p-level 6 \\\n",
    "    --o-collapsed-table genus.qza"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exporting your taxa table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zZPj_eiux9DP"
   },
   "source": [
    "Taxonomic feature tables in QIIME2 can be exported as .biom files, which can then be converted to .tsv files if desired. This is important, because some analysis methods not available on QIIME2 may require .biom or .tsv formats.\n",
    "\n",
    "We'll export this as a .tsv, which will be more usable for the next portion of the course that you'll see tomorrow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "deJ-FLcTx9Jq",
    "outputId": "fd93fc05-677b-4ff0-da72-eeeeaf5c906e"
   },
   "outputs": [],
   "source": [
    "!qiime tools export \\\n",
    "    --input-path genus.qza \\\n",
    "    --output-path exported\n",
    "    \n",
    "!biom convert -i exported/feature-table.biom -o genus.tsv --to-tsv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_1jjxJv5Pzwr"
   },
   "source": [
    "Let's peek at the results 🔭"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 474
    },
    "id": "nrq3gtn42kWv",
    "outputId": "5e8e3328-325f-4ceb-9330-d29beeab7112"
   },
   "outputs": [],
   "source": [
    "abundances = pd.read_table(\"genus.tsv\", skiprows=1, index_col=0)\n",
    "abundances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ROKU7_tcyFKD"
   },
   "source": [
    "This is easier to interpret by visualizing the results. We can use the file we just exported from QIIME2 to build a visualization using any tool we like, such as seaborn or plotnine.\n",
    "\n",
    "Here is an example of building a visualization (a heatmap) in seaborn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "vC5qKf680fAa"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 548
    },
    "id": "YOHG-0nXquu5",
    "outputId": "878d6d19-62e1-43ef-ff02-feca30ee0c09"
   },
   "outputs": [],
   "source": [
    "abund_to_plot = abundances.copy()\n",
    "abund_to_plot.index = abund_to_plot.index.str.split(\";\").str[5]          # Use only the genus name\n",
    "abund_to_plot = abund_to_plot[~abund_to_plot.index.isin([\"g__\", \"__\"])]  # remove unclassified genera\n",
    "abund_to_plot = abund_to_plot.sample(50, axis=0)                         # use 50 random genera (rows)\n",
    "\n",
    "# Let's do a centered log-ratio transform: log x_i - log mean(x)\n",
    "transformed = abund_to_plot.apply(\n",
    "    lambda xs: np.log(xs + 0.5) - np.log(xs.mean() + 0.5),\n",
    "    axis=0)\n",
    "\n",
    "sns.clustermap(transformed.T, cmap=\"magma\", xticklabels=True, figsize=(16, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GfzsEjtS0fAb"
   },
   "source": [
    "Now, our data is starting to be interpretable. Each row is a sample, and each column is a bacterial genus. The table values are \"counts\", or the number of times a genus was detected in a certain sample. We can use relative abundance data to test hypotheses, but it requires special statistical methods because it is ___compositional___.\n",
    "\n",
    "Let's discuss what we can do about that in the [presentation](https://gibbons-lab.github.io/isb_course_2024/16S)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-N5E8lW40fAb"
   },
   "source": [
    "---\n",
    "# Differential Abundance Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at abundances in the _complete_ dataset, so you can see the abundance distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to be switching to python to show you these transformations, so we'll need a few more packages 📦"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Counts to relative abundance:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we do any testing, we need to divide genus counts by total sample counts. This gives us _relative abundance_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "complete_abundances = pd.read_table('complete_study/genus.tsv', skiprows=1, index_col=0)\n",
    "complete_abundances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rel_abund = complete_abundances.divide(complete_abundances.sum(axis=0), axis=1)     # divide all values in a row by the sum of the row\n",
    "rel_abund"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Preprocessing:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do any rows have entirely zero? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rows_with_zeros = rel_abund.eq(0).all(axis=1)\n",
    "count_row_zeros = rows_with_zeros.sum()\n",
    "count_row_zeros"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many rows (genera) are not present in most samples (zero in >80% of samples?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# identify sparse rows\n",
    "num_zeros_row = rel_abund.eq(0).sum(axis=1)\n",
    "num_cols = rel_abund.shape[1]\n",
    "sparse_rows = num_zeros_row > (0.8 * num_cols)\n",
    "print(sparse_rows.sum())\n",
    "\n",
    "# remove sparse rows \n",
    "abund_filt = rel_abund[~sparse_rows]\n",
    "\n",
    "# transpose to make samples rows rather than genera (match formatting of our metadata)\n",
    "abund_filt = abund_filt.T\n",
    "abund_filt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the column names of genera in a list, then merge the abundance and metadata tables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "genus_columns = abund_filt.columns                                                            # get a list of abundance column names\n",
    "complete_metadata.set_index('id', inplace=True)                                               # change the index (row names) of metadata to the column 'id'\n",
    "merged_rel_abund = complete_metadata.merge(abund_filt, left_index=True, right_index=True)     # merge the two dataframes by their index (row names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a histogram of Faecalibacterium split by PD status\n",
    "\n",
    "sns.displot(data=merged_rel_abund,\n",
    "            x='k__Bacteria;p__Firmicutes;c__Clostridia;o__Eubacteriales;f__Oscillospiraceae;g__Faecalibacterium',\n",
    "            hue='parkinson_disease')\n",
    "\n",
    "plt.xlabel('Faecalibacterium relative abundance')\n",
    "plt.ylabel('Number of samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1. Nonparametric tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wanted to test whether Faecalibacterium specifically was differentially abundant in PD, we could do this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# split the variable of interest into 2 lists by PD status\n",
    "faecali_pd = merged_rel_abund[merged_rel_abund['parkinson_disease'] == 'Yes']['k__Bacteria;p__Firmicutes;c__Clostridia;o__Eubacteriales;f__Oscillospiraceae;g__Faecalibacterium']\n",
    "faecali_hc = merged_rel_abund[merged_rel_abund['parkinson_disease'] == 'No']['k__Bacteria;p__Firmicutes;c__Clostridia;o__Eubacteriales;f__Oscillospiraceae;g__Faecalibacterium']\n",
    "\n",
    "u_stat, p_value = stats.mannwhitneyu(faecali_pd, faecali_hc, alternative='two-sided')\n",
    "\n",
    "print(f\"U-statistic: {u_stat}\")\n",
    "print(f\"p-value: {p_value}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would see that yes, if we only perform one test, it appears that Faecalibacterium mean relative abundance is significantly different between PD and healthy controls (with p < 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Normalization + parametric tests:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can do a CLR transform, we have to add a pseudocount, because log(0) does not exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# keep track of which values are nonzero\n",
    "abund_nonzero = abund_filt > 0\n",
    "\n",
    "# add a pseudocount\n",
    "abund_filt_pseudo = abund_filt + 0.0000001\n",
    "abund_filt_pseudo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have pseudocounts, we can do the CLR transform. Here, we use the scipy.stats.mstats function gmean and compute the CLR manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a new col for gmean\n",
    "abund_clr = abund_filt_pseudo.copy()\n",
    "abund_clr['gmean'] = abund_clr.apply(stats.mstats.gmean, axis=1, nan_policy='omit')\n",
    "\n",
    "# actual CLR\n",
    "abund_clr = np.log(abund_clr.divide(abund_clr['gmean'], axis=0))\n",
    "\n",
    "# check gmean to make sure CLR went correctly\n",
    "# gmean should be all zeros (because log(x)/log(x) = log(x/x) = log(1) = 0)\n",
    "abund_clr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop gmean\n",
    "abund_clr.drop(columns='gmean', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's see what the transformation did to our data. We can make the same plot as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.displot(data=abund_clr,\n",
    "            x='k__Bacteria;p__Firmicutes;c__Clostridia;o__Eubacteriales;f__Oscillospiraceae;g__Faecalibacterium')\n",
    "plt.xlabel('Faecalibacterium CLR')\n",
    "plt.ylabel('Number of samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we have a bimodal distribution due to the presence of zeros. At this point, the most conservative option is to remove these values, because we cannot be sure if they are true zeros. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop the samples that were zero (make NA so they will not be plotted)\n",
    "abund_clr_nonzero = abund_clr[abund_nonzero]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now try plotting again\n",
    "sns.displot(data=abund_clr_nonzero,\n",
    "            x='k__Bacteria;p__Firmicutes;c__Clostridia;o__Eubacteriales;f__Oscillospiraceae;g__Faecalibacterium')\n",
    "plt.xlabel('Faecalibacterium CLR')\n",
    "plt.ylabel('Number of samples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our data is _normal-ish_ (though still a bit skewed) but probably safe to use with parametric tests! This means we are ready to merge it back to the metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_clr = complete_metadata.merge(abund_clr_nonzero, left_index=True, right_index=True)     # merge the two dataframes by their index (row names)\n",
    "merged_clr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code to perform a t-test is very similar to that of a Mann-Whitney U-test, but we need to also remove NAs.\n",
    "\n",
    "We can accomplish this by adding `dropna()` to the end of our group assignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# how we would do one t-test\n",
    "faecali_clr_pd = merged_clr[merged_clr['parkinson_disease'] == 'Yes']['k__Bacteria;p__Firmicutes;c__Clostridia;o__Eubacteriales;f__Oscillospiraceae;g__Faecalibacterium'].dropna()\n",
    "faecali_clr_hc = merged_clr[merged_clr['parkinson_disease'] == 'No']['k__Bacteria;p__Firmicutes;c__Clostridia;o__Eubacteriales;f__Oscillospiraceae;g__Faecalibacterium'].dropna()\n",
    "\n",
    "t_stat, p_val = stats.ttest_ind(faecali_clr_pd, faecali_clr_hc, equal_var=False)\n",
    "\n",
    "print(f\"T-stat: {t_stat}\")\n",
    "print(f\"p-value: {p_val}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In people with _measurable_ amounts of Faecalibacterium, people with PD have significantly lower Faecalibacterium than people without. \n",
    "\n",
    "Notice how in trimming our data, we are now testing a different question."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiple testing and FDR correction:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this much interesting data, it would be pretty tedious to test variables one-by-one. This is why we made a list `genus_column` earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty list to store test results\n",
    "t_test_results = []\n",
    "\n",
    "# Loop through each genus column and perform the t-test\n",
    "for genus in genus_columns:\n",
    "    # Separate the data for Parkinson's disease 'Yes' and 'No' groups\n",
    "    pd_group = merged_clr[merged_clr['parkinson_disease'] == 'Yes'][genus].dropna()\n",
    "    hc_group = merged_clr[merged_clr['parkinson_disease'] == 'No'][genus].dropna()\n",
    "    \n",
    "    # Perform the t-test\n",
    "    if len(pd_group) > 0 and len(hc_group) > 0:  # Only perform test if there are values in both groups\n",
    "        t_stat, p_val = stats.ttest_ind(pd_group, hc_group, equal_var=False)\n",
    "        \n",
    "        # Store the results\n",
    "        t_test_results.append({\n",
    "            'genus': genus,\n",
    "            't_stat': t_stat,\n",
    "            'p_val': p_val,\n",
    "            'n_pd': len(pd_group),\n",
    "            'n_hc': len(hc_group)\n",
    "        })\n",
    "\n",
    "# Convert the results into a DataFrame\n",
    "ttest_results_df = pd.DataFrame(t_test_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a look at our results dataframe\n",
    "ttest_results_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does our p-value distribution look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(data=ttest_results_df, x='p_val', bins=13)\n",
    "plt.xlabel('p-value')\n",
    "plt.ylabel('number of t-tests')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks ok? We expect a left skew given both true and false positives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we have a list of p-values, we can easily do FDR correction. We will use a Benjamini-Hochberg (BH) correction, which is less conservative than Bonferroni. The BH corrected p-value, q, is calculated as a function of the entire list of q-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.multitest import multipletests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# do the FDR correction (this returns q-values in the order they were given)\n",
    "qvals = multipletests(ttest_results_df['p_val'], method='fdr_bh')[1]\n",
    "\n",
    "# append q-values to dataframe\n",
    "ttest_results_df['q_val'] = qvals\n",
    "\n",
    "# show the new dataframe, but sort by q-value\n",
    "ttest_results_df.sort_values(by='q_val')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's pretty cool! Let's separate just the significant hits, and make a column with some easier-to-plot taxonomic names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ttest_signif = ttest_results_df[ttest_results_df['q_val'] < 0.1]\n",
    "ttest_signif.loc[:, 'genus_plot'] = ttest_signif['genus'].str.split(';').apply(lambda x: ';'.join(x[-2:])) # take only the last 2 taxonomic levels\n",
    "ttest_signif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How can we visualize this information? Here, I've plotted the T-statistic for each significant hit, colored by q-value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort by absolute t-statistic for better visualization\n",
    "ttest_signif = ttest_signif.sort_values(by='t_stat', key=abs, ascending=False)\n",
    "\n",
    "# make a custom color palette that goes from light to dark based on the q-values\n",
    "norm = plt.Normalize(0.0, 0.1)\n",
    "cmap = plt.cm.Blues_r  # color map\n",
    "bar_colors = cmap(norm(ttest_signif['q_val'])).tolist()\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "sns.barplot(data=ttest_signif, x='t_stat', y='genus_plot', palette=bar_colors, hue='genus_plot', ax=ax, legend=False)\n",
    "\n",
    "# Make labels\n",
    "ax.set_title('Differentially Abundant Genera in PD patients')\n",
    "ax.set_xlabel('T-Statistic')\n",
    "ax.set_ylabel('Genus')\n",
    "\n",
    "# make a key for q-values\n",
    "sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)  # create ScalarMappable object\n",
    "sm.set_array([])  # Empty array for ScalarMappable\n",
    "cbar = plt.colorbar(sm, ax=ax)\n",
    "cbar.set_label('q-value')\n",
    "\n",
    "# Show the plot\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WPMlG8w9aCsW",
    "tags": []
   },
   "source": [
    "---\n",
    "# Exercise: Plant a Tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hZrVCQE80fAc"
   },
   "source": [
    "One visualization that we did not spend a lot of time on was the phylogentic tree of our ASVs. Let's change that! We have seen that there are genera that appear in multiple populations in the previous step. But are the organisms in that genus actually the same?\n",
    "\n",
    "Let's annotate the tree with our taxonomic classifications and abundances. We will use the empress plugin again but this time with the `community-plot` option. I filled in a template of the command for you. Can you figure out what has to go in the empty spaces?\n",
    "\n",
    "**QUESTIONS:**\n",
    "\n",
    "1) Are some of the branch lengths on the tree longer than you would expect? Do you notice anything interesting or suspicious about the taxonomic identities of these branches?\n",
    "\n",
    "2) Can you find examples of phyla that are polyphyletic (i.e. where clusters of ASVs from the same phylum are found in different locations on the tree, showing different commmon ancestors)? What about polyphyletic taxa at lower taxonomic levels, like at the family or genus levels? Why do you think these patterns exist?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qf-bY5aci4bv"
   },
   "outputs": [],
   "source": [
    "# This won't run until you fill in the [EMPTY] spots with the right files ;)\n",
    "\n",
    "!qiime empress community-plot \\\n",
    "    --i-tree [EMPTY] \\\n",
    "    --i-feature-table dada/table.qza \\\n",
    "    --m-sample-metadata-file [EMPTY] \\\n",
    "    --m-feature-metadata-file taxa.qza \\\n",
    "    --o-visualization community-tree-viz.qzv"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}