16S_2022.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Gibbons-Lab/isb_course_2022/blob/main/16S_2022.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nQRlv0ypvthh"
      },
      "source": [
        "# 🦠 Amplicon Sequencing Data Analysis with Qiime 2\n",
        "\n",
        "This notebook will accompany the first session of the 2022 ISB Virtual Microbiome Series. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2022/16S). \n",
        "\n",
        "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",
        "**Disclaimer:**\n",
        "\n",
        "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. \n",
        "\n",
        "## Setup\n",
        "\n",
        "QIIME 2 is usually installed by following the [official installation instructions](https://docs.qiime2.org/2022.8/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."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QFsZdvOuvZmq",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!git clone https://github.com/gibbons-lab/isb_course_2022 materials"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2eenIEY76mcv"
      },
      "source": [
        "We will switch to working within the `materials` directory for the rest of the notebook."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "h3NALJ7u6mBP",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "%cd materials"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fQVnDjFauZ-n"
      },
      "source": [
        "Now we are ready to set up our environment. This will take about 10-15 minutes. \n",
        "\n",
        "**Note**: This setup is only relevant for Google Colaboratory and will not work on your local machine. Please follow the [official installation instructions](https://docs.qiime2.org/2021.8/install/) for that."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "BAyTazMCvrpf",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "%run setup_qiime2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BGWsqUfq5oSi"
      },
      "source": [
        "---\n",
        "\n",
        "While this is running let's get back to the presentation. 📚\n",
        "\n",
        "---\n",
        "\n",
        "<br>\n",
        "\n",
        "## Our first QIIME 2 command\n",
        "\n",
        "Let's remember our workflow for today.\n",
        "\n",
        "![our workflow](https://github.com/Gibbons-Lab/isb_course_2022/raw/main/docs/16S/assets/steps.png)\n",
        "\n",
        "The first thing we have to do is to get the data into an 'artifact'.\n",
        "We can import the data with the `import` action from the tools. For that we have to give\n",
        "QIIME 2 a *manifest* (list of raw files) and tell it what *type of data* we\n",
        "are importing and what *type of artifact* we want. \n",
        "\n",
        "**QoL Tip:** QIIME 2 commands can get very long. To split them up over several lines we can use `\\` which means \"continue on the next line\"."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6mGgUrRf5pyb",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime tools import \\\n",
        "  --type 'SampleData[PairedEndSequencesWithQuality]' \\\n",
        "  --input-path data/manifest.tsv \\\n",
        "  --output-path sequences.qza \\\n",
        "  --input-format PairedEndFastqManifestPhred33V2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lw7TmzbiAoLj"
      },
      "source": [
        "Since we have quality information for the sequencing reads, let's also generate\n",
        "our first visualization to inspect sequence quality. \n",
        "\n",
        "---\n",
        "\n",
        "Qiime 2 commands can become pretty long. Here are some pointers to remember the\n",
        "structure of a command:\n",
        "\n",
        "```\n",
        "qiime plugin action --i-argument1 ... --o-argument2 ...\n",
        "```\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\n",
        "\n",
        "---\n",
        "\n",
        "In this case we will use the `summarize` action from the `demux` plugin with the previously generated artifact as input and output the resulting visualization to the `qualities.qzv` file."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "feGa41cEApGs",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime demux summarize --i-data sequences.qza --o-visualization qualities.qzv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O2GenFy3C8gu"
      },
      "source": [
        "You can view the plot by downloading the .qzv 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.\n",
        "\n",
        "🤔 What do you observe across the read? Where would you truncate the reads?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "os5MeRvLIq0v"
      },
      "source": [
        "# Denoising amplicon sequence variants\n",
        "\n",
        "We will now run the DADA2 plugin which will do 3 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",
        "Since this step takes a bit, let's start the process and use the time to\n",
        "understand what is happening:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RdLSMp9-Dl-d",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime dada2 denoise-paired \\\n",
        "    --i-demultiplexed-seqs sequences.qza \\\n",
        "    --p-trunc-len-f 150 \\\n",
        "    --p-trunc-len-r 150 \\\n",
        "    --p-n-threads 2 \\\n",
        "    --output-dir dada2 --verbose"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aFpQQ5GHdbzX"
      },
      "source": [
        "---\n",
        "\n",
        "While this is running let's get back to the presentation. 📚\n",
        "\n",
        "---\n",
        "\n",
        "<br>\n",
        "\n",
        "\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": "PFrfo-uCdoRz",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "# obscure magic that will only copy if the previous command failed\n",
        "![ -d dada2 ] || cp -r treasure_chest/dada2 ."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4YzrovxVj7U_"
      },
      "source": [
        "Ok, this step ran, but we should also make sure it kind of worked. One good way to tell if the identified ASVs are representative of the sample is to see how many reads were maintained throughout the pipeline. Here, the most common issues and solutions are:\n",
        "\n",
        "**Large fraction of reads is lost during merging (only paired-end)**\n",
        "\n",
        "![read overlap](https://gibbons-lab.github.io/isb_course_2022/16S/assets/read_overlap.png)\n",
        "\n",
        "In order to merge ASVs DADA2 uses an overlap of 12 bases between forward and reverse reads by default. Thus, your reads must allow for sufficient overlap *after* trimming. So if your amplified region is 450bp long and you have 2x250bp reads and you trim the last 30 bases of each read, truncating the length to 220bp, the total length of covered sequence is 2x220 = 440 which is shorter than 450bp so there will be no overlap. To solve this issue trim less of the reads or adjust the `--p-min-overlap` parameters to something lower (but not too low).\n",
        "\n",
        "<br>\n",
        "\n",
        "**Most of the reads are lost as chimeric**\n",
        "\n",
        "![read overlap](https://gibbons-lab.github.io/isb_course_2022/16S/assets/chimera.png)\n",
        "\n",
        "This is usually an experimental issue as chimeras are introduced during amplification. If you can adjust your PCR, try to run fewer cycles. Chimeras can also be introduced by incorrect merging. If your minimum overlap is too small ASVs may be merged randomly. Possible fixes are to increase the `--p-min-overlap` parameter or run the analysis on the forward reads only (in our empirical observations, chimeras are more likely to be introduced in the joined reads). *However, losing between 5-25% of your reads to chimeras is normal and does not require any adjustments.*\n",
        "\n",
        "Our denoising stats are contained in an artifact. To convert it to a visualization we can use `qiime metadata tabulate`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "n_Tr3Ujcj61w",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime metadata tabulate \\\n",
        "    --m-input-file dada2/denoising_stats.qza \\\n",
        "    --o-visualization dada2/denoising-stats.qzv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o75g7DSO2tu5"
      },
      "source": [
        "What proprotion of reads was retained throughout the entire pipeline? Look at the final number of used reads (non-chimeric). What do you observe when comparing those values between samples and how might that affect diversity metrics?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "c0ZV3_sGQ20u"
      },
      "source": [
        "# Phylogeny and ecological diversity metrics"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aEit8p7V9Ht9"
      },
      "source": [
        "## Building a tree\n",
        "\n",
        "We can build a phylogenetic tree for our sequences using the following command:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2PVNlD0g9MsX",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime phylogeny align-to-tree-mafft-fasttree \\\n",
        "    --i-sequences dada2/representative_sequences.qza \\\n",
        "    --output-dir tree"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4lB8NI0IiJgk"
      },
      "source": [
        "We can create a visualization for the tree using the [empress](https://github.com/biocore/empress) Qiime 2 plugin."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cSS_y3Kligrc",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime empress tree-plot \\\n",
        "    --i-tree tree/rooted_tree.qza \\\n",
        "    --o-visualization tree/empress.qzv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UmxTRou99a_n"
      },
      "source": [
        "This looks nice but is not particularly informative as is. The main utility of our tree will be in complementing our diversity analyses. It will tell us which ASVs are more or less related to one another, which will allow us to calculate different kinds of ecological diversity metrics. But what are those?\n",
        "\n",
        "---\n",
        "\n",
        "Let's get back to the presentation. 📚\n",
        "\n",
        "---\n",
        "\n",
        "<br>\n",
        "\n",
        "## Alpha and Beta Diversity\n",
        "\n",
        "![sample sources](https://github.com/Gibbons-Lab/isb_course_2022/raw/main/docs/16S/assets/sample_sources.png)\n",
        "\n",
        "One of our main goals will be to compare the microbial composition across different populations. QIIME 2 has \"run-all\" command for diversity analyses. This will\n",
        "\n",
        "1. Subsample our samples to the same total number of reads (Why?) \n",
        "2. Calculate alpha and beta diversity measures\n",
        "3. Visualize PCoA projections"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p-uZBNVX-y2L"
      },
      "source": [
        "## Statistical analyses\n",
        "\n",
        "Let's first have a look at alpha diversity. Can we see a difference in the per-sample diversity across populations? "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0js6xHuw-NqF",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime diversity core-metrics-phylogenetic \\\n",
        "    --i-table dada2/table.qza \\\n",
        "    --i-phylogeny tree/rooted_tree.qza \\\n",
        "    --p-sampling-depth 8000 \\\n",
        "    --m-metadata-file data/metadata.tsv \\\n",
        "    --output-dir diversity"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GnkeX6iY--g-",
        "vscode": {
          "languageId": "python"
        }
      },
      "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": "ngH5fUeDQMO0"
      },
      "source": [
        "Now, let's use beta diversity to see how different the samples are from one another. First download `diversity/weighted_unifrac_emperor.qzv` and take a look. Do samples separate based on the population?\n",
        "\n",
        "We can check whether that separation is 'significant' by using a PERMANOVA test."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KzNW92riQguX",
        "vscode": {
          "languageId": "python"
        }
      },
      "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 \"ethnic_group\" \\\n",
        "    --p-n-jobs 2 \\\n",
        "    --o-visualization diversity/permanova.qzv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5OgGG9yabL4d"
      },
      "source": [
        "---\n",
        "\n",
        "Let's get back to the presentation. 📚\n",
        "\n",
        "---\n",
        "\n",
        "<br>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FB7WYPZ1DEoR"
      },
      "source": [
        "# Taxonomy\n",
        "\n",
        "We will use a Bayes classifier trained on the NCBI Refseq database version 212. However, there are many more databases which can be downloaded from https://docs.qiime2.org/2022.8/data-resources/. For instance, the SILVA database is very popular and probably the largest 16S database available right now."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "EEBbgvnhDGTs",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime feature-classifier classify-sklearn \\\n",
        "    --i-reads dada2/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": "xISxzBo6Mpka"
      },
      "source": [
        "Now let's have a look at the relative abundances of the different bacterial taxa we have in each sample:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4VTL2yEBMqfO",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime taxa barplot \\\n",
        "    --i-table dada2/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": "4_FVM8VfX_yQ"
      },
      "source": [
        "We can also collapse data on a particular taxonomic rank using the QIIME 2 [taxa plugin](https://docs.qiime2.org/2021.4/plugins/available/taxa/). Why might we want to look at different taxonomic ranks, rather than just looking at ASVs?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gl72G8KaYmBm",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime taxa collapse \\\n",
        "    --i-table treasure_chest/dada2/table.qza \\\n",
        "    --i-taxonomy treasure_chest/taxa.qza \\\n",
        "    --p-level 6 \\\n",
        "    --o-collapsed-table genus.qza"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dlVdzgfbZcjt"
      },
      "source": [
        "We can export the table and convert it to a .csv file so that we can analyze these data using tools outside of the QIIME 2 environment."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZsxKf3yAZiHk",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "!qiime tools export \\\n",
        "    --input-path genus.qza \\\n",
        "    --output-path exported\n",
        "!biom convert -i exported/feature-table.biom -o genus.tsv --to-tsv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TezGLSEGZrcn"
      },
      "source": [
        "Now the data are in a common format and we can use them, for instance, to draw a heatmap using Pandas and Seaborn. No need to understand every bit of code here. This just serves to illustrate that you can get data out of QIIME 2 for custom visualizations."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "p5zUjU8JZqL4",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "import seaborn as sns\n",
        "\n",
        "abundances = pd.read_table(\"genus.tsv\", skiprows=1, index_col=0)\n",
        "abundances.index = abundances.index.str.split(\";\").str[5]       # Use only the genus name\n",
        "abundances = abundances[~abundances.index.isin([\"g__\", \"__\"])]  # remove unclassified genera\n",
        "abundances = abundances.sample(50)                              # use 50 random genera\n",
        "\n",
        "# Let's do a centered log-ratio transform: log x_i - log mean(x)\n",
        "transformed = abundances.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=(18, 6))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QbRoCh1BR0GZ"
      },
      "source": [
        "# Exercises\n",
        "\n",
        "Okay, that's enough time in the back seat. \n",
        "\n",
        "It's time to take the wheel 🚗 \n",
        "\n",
        "Now you can dive into the data 🏊"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hy5HN4ZoS2-9"
      },
      "source": [
        "## Exercise - Plant a Tree\n",
        "\n",
        "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?\n",
        "\n",
        "3) Use the \"Shear Tree\" functionality to only visualize the phylogenetic tree of the genus *Prevotella*. What do you observe? Can you figure out which of the unique V4 sequences are found only in one population? What impact might this have on studies trying to link *Prevotella* to human health?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tq9TWRMjkER5",
        "vscode": {
          "languageId": "python"
        }
      },
      "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 dada2/table.qza \\\n",
        "    --m-sample-metadata-file d[EMPTY] \\\n",
        "    --m-feature-metadata-file taxa.qza \\\n",
        "    --o-visualization community-tree-viz.qzv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tFuoGkWtptmR"
      },
      "source": [
        "<br><br><br>\n",
        "\n",
        "---\n",
        "\n",
        "# 🦠 Space for your own exploration\n",
        "\n",
        "Feel free to use this space to run some additional analyses. If you rerun all the previous cells in the notebook (`Runtime > Run All`) you should be all set to start working on the previously generated artifacts. Remember that you can get help in the Slack channel and also in the [Qiime 2 documentation](https://docs.qiime2.org)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pq919yDYqYoe",
        "vscode": {
          "languageId": "python"
        }
      },
      "outputs": [],
      "source": [
        "# You can add more code cells with the \"+ Code\" button on the top right\n"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}