16S.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "16S.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "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_2020/blob/master/16S.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 2020 ISB Microbiome Course. The presentation slides can be [found here](https://gibbons-lab.github.io/isb_course_2020/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 promted to certify the notebook is safe. We promise that it is 🤞\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 leading `!`. 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/2021.4/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",
      "metadata": {
        "id": "QFsZdvOuvZmq"
      },
      "source": [
        "!git clone https://github.com/gibbons-lab/isb_course_2020 materials"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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.4/install/) for that."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BAyTazMCvrpf"
      },
      "source": [
        "%run materials/setup_qiime2.py"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "h3NALJ7u6mBP"
      },
      "source": [
        "%cd materials"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BGWsqUfq5oSi"
      },
      "source": [
        "## 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_2020/raw/master/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",
      "metadata": {
        "id": "6mGgUrRf5pyb"
      },
      "source": [
        "!qiime tools import \\\n",
        "  --type 'SampleData[SequencesWithQuality]' \\\n",
        "  --input-path manifest.tsv \\\n",
        "  --output-path cdiff.qza \\\n",
        "  --input-format SingleEndFastqManifestPhred33V2"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "feGa41cEApGs"
      },
      "source": [
        "!qiime demux summarize --i-data cdiff.qza --o-visualization qualities.qzv"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "RdLSMp9-Dl-d"
      },
      "source": [
        "!qiime dada2 denoise-single \\\n",
        "    --i-demultiplexed-seqs cdiff.qza \\\n",
        "    --p-trunc-len 150 \\\n",
        "    --output-dir dada2 --verbose"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aFpQQ5GHdbzX"
      },
      "source": [
        "If this step takes too long or fails, you can also copy the results from the treasure chest. *However, don't run the next cell if this cell completes successfully.*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PFrfo-uCdoRz"
      },
      "source": [
        "# only run if the previous cell did not work\n",
        "!cp -r treasure_chest/dada2 ."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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_2020/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_2020/16S/assets/chimera.png)\n",
        "\n",
        "This is usualluy 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",
      "metadata": {
        "id": "n_Tr3Ujcj61w"
      },
      "source": [
        "!qiime metadata tabulate \\\n",
        "    --m-input-file dada2/denoising_stats.qza \\\n",
        "    --o-visualization dada2/denoising-stats.qzv"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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 diversity"
      ]
    },
    {
      "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",
      "metadata": {
        "id": "2PVNlD0g9MsX"
      },
      "source": [
        "!qiime phylogeny align-to-tree-mafft-fasttree \\\n",
        "    --i-sequences dada2/representative_sequences.qza \\\n",
        "    --output-dir tree"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UmxTRou99a_n"
      },
      "source": [
        "You can visualize your tree using iTOL (https://itol.embl.de/). Open iTol and upload the artifact from `materials/tree/tree.qza`.\n",
        "\n",
        "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.\n",
        "\n",
        "## Alpha and Beta Diversity\n",
        "\n",
        "Qiime 2 has \"run-all\" command for diversity analyses. This will\n",
        "\n",
        "1. Subsampe 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": "code",
      "metadata": {
        "id": "0js6xHuw-NqF"
      },
      "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 metadata.tsv \\\n",
        "    --output-dir diversity"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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 between healthy and sick individuals?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GnkeX6iY--g-"
      },
      "source": [
        "!qiime diversity alpha-group-significance \\\n",
        "    --i-alpha-diversity diversity/shannon_vector.qza \\\n",
        "    --m-metadata-file metadata.tsv \\\n",
        "    --o-visualization diversity/alpha_groups.qzv"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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 disease state? \n",
        "\n",
        "We can check whether that separation is 'significant' by using a PERMANOVA test."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KzNW92riQguX"
      },
      "source": [
        "!qiime diversity adonis \\\n",
        "    --i-distance-matrix diversity/weighted_unifrac_distance_matrix.qza \\\n",
        "    --m-metadata-file metadata.tsv \\\n",
        "    --p-formula \"disease_stat\" \\\n",
        "    --p-n-jobs 2 \\\n",
        "    --o-visualization permanova.qzv"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FB7WYPZ1DEoR"
      },
      "source": [
        "# Taxonomy\n",
        "\n",
        "We will use a Bayes classifier trained on the GreenGenes database which can be downloaded from https://docs.qiime2.org/2021.4/data-resources/."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "drgYtfLbDosv"
      },
      "source": [
        "!wget https://data.qiime2.org/2021.4/common/gg-13-8-99-515-806-nb-classifier.qza"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EEBbgvnhDGTs"
      },
      "source": [
        "!qiime feature-classifier classify-sklearn \\\n",
        "    --i-reads dada2/representative_sequences.qza \\\n",
        "    --i-classifier gg-13-8-99-515-806-nb-classifier.qza \\\n",
        "    --o-classification taxa.qza"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "4VTL2yEBMqfO"
      },
      "source": [
        "!qiime taxa barplot \\\n",
        "    --i-table dada2/table.qza \\\n",
        "    --i-taxonomy taxa.qza \\\n",
        "    --m-metadata-file metadata.tsv \\\n",
        "    --o-visualization taxa_barplot.qzv"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "gl72G8KaYmBm"
      },
      "source": [
        "!qiime taxa collapse \\\n",
        "    --i-table dada2/table.qza \\\n",
        "    --i-taxonomy taxa.qza \\\n",
        "    --p-level 6 \\\n",
        "    --o-collapsed-table genus.qza"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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",
      "metadata": {
        "id": "ZsxKf3yAZiHk"
      },
      "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"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p5zUjU8JZqL4"
      },
      "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",
        "\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=1)\n",
        "\n",
        "sns.clustermap(transformed.T, cmap=\"magma\", xticklabels=True, figsize=(16, 5))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "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": "Moq2YOvDSDap"
      },
      "source": [
        "## Exercise 1 - Sample Classification\n",
        "\n",
        "One pretty basic question we can ask is whether the microbial composition in the gut is predictive of recurrent *C. difficile* infection. \n",
        "Could you predict if someone had a recurrent infection from 16S data alone? \n",
        "\n",
        "Metadata about the samples is again contained in `metadata.tsv` which looks like this:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T16jC6U1ONow"
      },
      "source": [
        "import pandas as pd\n",
        "\n",
        "pd.read_csv(\"metadata.tsv\", sep=\"\\t\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j1Q2ZpBEON42"
      },
      "source": [
        "Start with the `classify-samples` action and follow it up by finding and looking at the `heatmap` visualization afterwards that shows important genera. You can also find accuracy results that will tell you how well the phenotype can be predicted from the microbial composition."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Nhf1g2RpSCak"
      },
      "source": [
        "# Your code here"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hy5HN4ZoS2-9"
      },
      "source": [
        "## Exercise 2 - 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! \n",
        "\n",
        "First, annotate the tree with taxonomy by going back to the tree (uploading `tree.qza`) and visualizing it using iTOL (https://itol.embl.de/). You can add the predicted taxa (`taxa.qza`) using the `+` button on the lower right of the visualization. \n",
        "\n",
        "Also add your feature table, but don't use the raw one (why?). Rather, use `diversity/rarefied_table.qza` where each sample has been subsampled to 10000 reads (why?). \n",
        "\n",
        "What other normalization methods could you use?\n",
        "\n",
        "What do you observe?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bx6qkBT-85VE"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}