16S_2023.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_2023/blob/main/16S_2023.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2inlS7SNn0-o"
      },
      "source": [
        "# 🦠 Amplicon Sequencing Data Analysis with QIIME 2\n",
        "\n",
        "This notebook will accompany the first session of the 2023 ISB Virtual Microbiome Series. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2023/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/2023.7/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": "9Cyg41Tcn1Ua"
      },
      "outputs": [],
      "source": [
        "!git clone https://github.com/gibbons-lab/isb_course_2023 materials"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MQ9XVvOToEBR"
      },
      "source": [
        "This repository, called __materials__, contains all the relevant data and other resources we'll need for this course. Let's navigate to that directory now:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GvNHs-EroFNM"
      },
      "outputs": [],
      "source": [
        "%cd materials"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Notice here we use ```%``` instead of ```!``` to run out 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"
      ],
      "metadata": {
        "id": "IHvKnTbdypyp"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rIOk9J0AoYXX"
      },
      "source": [
        "## Install QIIME2\n",
        "Now that we have all our materials, we're _almost_ ready to get started - but not quite. Remember QIIME2? 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": {
        "id": "Inkvkt_DoSQT"
      },
      "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_2023/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"
      },
      "source": [
        "# Let's Get Started!\n",
        "\n",
        "Now we're on to the fun part. Let's begin by taking a look at our data. In the _data_ folder, you'll find eight FASTQ files, a file manifest, and a metadata file. Let's take a look at the manifest, first. This is a file that contains the name, and filepath of all our samples, and we'll need it later on when we use QIIME2 📝."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_pFaCnO6iqqW"
      },
      "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 the metadata file, which will give us more context on our samples 🔬"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "sAfnoRnxgZOO"
      },
      "outputs": [],
      "source": [
        "metadata = pd.read_csv('data/metadata.tsv', sep = '\\t')\n",
        "metadata"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qc10jMUYsXyh"
      },
      "source": [
        "Looks good, all eight FASTQ files are accounted for, four healthy and four with recurrent CDI. We can use the manifest file to import our files into QIIME2.\n",
        "\n",
        "## QIIME2 Pipeline\n",
        "\n",
        "Let's remind ourselves what the QIIME2 pipeline will do:\n",
        "![our workflow](https://github.com/Gibbons-Lab/isb_course_2023/raw/main/docs/16S/assets/steps.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HLF-zdrJXq7q"
      },
      "source": [
        "To use sequencing data in QIIME2, 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": {
        "id": "QmLebzqFX_0N"
      },
      "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 QIIME2 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\n",
        "\n",
        "---\n",
        "\n",
        "## Visualizing our Data 🔎\n",
        "\n",
        "Before we move on, let's use QIIME2 to visualize our sequencing data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nVPsDHTNvNUw"
      },
      "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.\n",
        "\n",
        "---\n",
        "\n",
        "## Quality Filtering\n",
        "\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 three 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": {
        "id": "HclmzHqyacLU"
      },
      "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": {
        "id": "EhU2aNz4_MDO"
      },
      "source": [
        "Let's check to see how that went. 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_2023/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_2023/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",
        "\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": "z-iW0UBd_Iiu"
      },
      "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.\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uDreq_YROCHQ"
      },
      "source": [
        "## Diversity and Phylogenetics\n",
        "An important metric to consider when studying microbial ecology is __diversity__. Diversity comes in two flavors: ⍺ (alpha) and β (beta).\n",
        "\n",
        "Alpha diversity is pretty simple - how diverse is a single sample? You might consider measures like richness and evenness.\n",
        "\n",
        "![alpha diversity](https://gibbons-lab.github.io/isb_course_2023/16S/assets/alpha_diversity.png)\n",
        "\n",
        "Beta diversity instead looks at how different two samples are from each other - what taxa are shared, and how their abundances differ.\n",
        "\n",
        "![beta diversity](https://gibbons-lab.github.io/isb_course_2023/16S/assets/beta_diversity.png)\n",
        "\n",
        "\n",
        "### Starting our Tree\n",
        "\n",
        "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": {
        "id": "Ym7blrPpKQuJ"
      },
      "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](https://github.com/biocore/empress) QIIME 2 plugin."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2Ttf0RTyhaPH"
      },
      "outputs": [],
      "source": [
        "!qiime empress tree-plot \\\n",
        "    --i-tree tree/rooted_tree.qza \\\n",
        "    --o-visualization tree/empress.qzv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H4fTVpeToufK"
      },
      "source": [
        "\n",
        "## Calculating Diversity\n",
        "Using the Diversity plugin, 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 rarefy all our samples. Since this randomly selects sequences, your results might look a little different. We'll also pass in our metadata file, so we can keep track how which samples come from each group."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1YqaTqS5huAz"
      },
      "outputs": [],
      "source": [
        "!qiime diversity core-metrics-phylogenetic \\\n",
        "    --i-table dada/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": "markdown",
      "metadata": {
        "id": "XslTNgFBg5qp"
      },
      "source": [
        "## Alpha Diversity\n",
        "\n",
        "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": {
        "id": "ne2WybKgh6fv"
      },
      "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": "VNmxg2lmjyHE"
      },
      "source": [
        "Like before, we can download the visualization and open it with the QIIME2 viewer.\n",
        "\n",
        "## Beta Diversity\n",
        "\n",
        "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",
        "<br>\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": {
        "id": "WMASNLaWjxwD"
      },
      "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 \"disease_state\" \\\n",
        "    --p-n-jobs 2 \\\n",
        "    --o-visualization diversity/permanova.qzv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CR2f0q4trGfA"
      },
      "source": [
        "---\n",
        "\n",
        "## Taxonomic Classification\n",
        "\n",
        "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 Bayesian classifier. Several such classifiers are available at https://docs.qiime2.org/2023.7/data-resources"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "firXAVIHrEzT"
      },
      "outputs": [],
      "source": [
        "!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": {
        "id": "RSRrxhX7x4dE"
      },
      "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": {
        "id": "XTClvRnXx4ob"
      },
      "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": {
        "id": "zZPj_eiux9DP"
      },
      "source": [
        "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": {
        "id": "deJ-FLcTx9Jq"
      },
      "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": "_1jjxJv5Pzwr"
      },
      "source": [
        "Let's peek at the results 🔭"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nrq3gtn42kWv"
      },
      "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. Here is an example of building a visualization (a heatmap) in seaborn:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YOHG-0nXquu5"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import seaborn as sns\n",
        "\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",
      "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?"
      ],
      "metadata": {
        "id": "WPMlG8w9aCsW"
      }
    },
    {
      "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": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}