05.Module_5_taxonomy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module 5: Taxonomy Sequence Classification\n",
    "\n",
    "## Overview\n",
    "\n",
    "Kraken is a taxonomic sequence classifier that assigns taxonomic labels to short DNA reads. It does this by examining the k-mers within a read and querying a database with those k-mers. This database contains a mapping of every k-mer in Kraken's genomic library to the lowest common ancestor (LCA) in a taxonomic tree of all genomes that contain that k-mer. The set of LCA taxa that correspond to the k-mers in a read are then analyzed to create a single taxonomic label for the read; this label can be any of the nodes in the taxonomic tree. Kraken is designed to be rapid, sensitive, and highly precise. Our tests on various real and simulated data have shown Kraken to have sensitivity slightly lower than Megablast with precision being slightly higher. On a set of simulated 100 bp reads, Kraken processed over 1.3 million reads per minute on a single core in normal operation, and over 4.1 million reads per minute in quick operation.\n",
    "\n",
    "![kraken](images/kraken.png)\n",
    "\n",
    "*Taken from: https://link.springer.com/article/10.1186/gb-2014-15-3-r46*\n",
    "\n",
    "Further reading: https://github.com/DerrickWood/kraken2/wiki/Manual\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install condacolab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -q condacolab\n",
    "import condacolab\n",
    "condacolab.install()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Install software"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install kraken2\n",
    "!conda install bioconda::kraken2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check if kraken2 is installed\n",
    "!kraken2 -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download SRA tools from bioconda\n",
    "!conda install -c bioconda sra-tools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://zenodo.org/records/14231070/files/Module_5.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract the .tar.gz file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar xvf Module_5.tar.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will navigate to the folder containing the paired FASTQ files ERR2667698_1.fastq.gz and ERR2667698_2.fastq.gz, which correspond to the run accession ERR2667694 from the project [PRJEB3084](https://www.ebi.ac.uk/ena/browser/view/PRJEB3084).\n",
    "\n",
    "Some important data about the sample:\n",
    "\n",
    "- Country of origin: Peru\n",
    "- Organism: *Streptococcus pneumoniae*\n",
    "- Instrument Platform: ILLUMINA\n",
    "- Instrument Model: Illumina MiSeq\n",
    "- Read Count: 2445004\n",
    "- Base Count: 369195604\n",
    "- Center Name: Wellcome Sanger Institute; SC\n",
    "- Library Layaout: PAIRED\n",
    "- Library strategy: WGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd Module_5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification\n",
    "\n",
    "To classify a set of sequences (reads), we will use the command: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run kraken2\n",
    "!kraken2 --use-names --db Kraken/k2_standard_08gb_20221209  --paired ERR2667698_1.fastq.gz ERR2667698_2.fastq.gz --report kraken_report.txt --output -"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An explanation of this command is as follows:\n",
    "\n",
    "`kraken2`: is the tool\n",
    "\n",
    "`--use-names`: Print scientific names instead of just taxids\n",
    "\n",
    "`--db` k2_standard_08gb_20221209: Name for Kraken 2 database\n",
    "\n",
    "`--paired` ERR2667698_1.fastq.gz ERR2667698_2.fastq.gz: Provide paired-end reads\n",
    "\n",
    "`--report` kraken_report.txt: Generate a report with aggregrate counts/clade to file\n",
    "\n",
    "`--output -`: Suppress output for standard output\n",
    "\n",
    "You can check the kraken output using `!head kraken_report.txt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should get:\n",
    "\n",
    "- **1st column**: Percentage of fragments covered by the clade rooted at this taxon\n",
    "\n",
    "- **2nd column**: Number of fragments covered by the clade rooted at this taxon\n",
    "\n",
    "- **3rd column**: Number of fragments assigned directly to this taxon\n",
    "\n",
    "- **4th column**: A rank code, indicating (U)nclassified, (R)oot, (D)omain, (K)ingdom, (P)hylum, (C)lass, (O)rder, (F)amily, (G)enus, or (S)pecies. Taxa that are not at any of these 10 ranks have a rank code that is formed by using the rank code of the closest ancestor rank with a number indicating the distance from that rank. E.g., \"G2\" is a rank code indicating a taxon is between genus and species and the grandparent taxon is at the genus rank\n",
    "\n",
    "- **5th column**: NCBI taxonomic ID number\n",
    "\n",
    "- **6th column**: Indented scientific name\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the percentage of fragments covered by the *Streptoccocus pneumoniae*, run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!grep \"Streptococcus pneumoniae\" kraken_report.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should get:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QC of Taxonomy\n",
    "\n",
    "Now we will calculate the percentage of *Streptococcus pneumoniae* reads and based on that we will assess if the sample is **Pass** or **Fail**\n",
    "\n",
    "The cut-off limit for *Streptococcus pneumoniae* is 60%\n",
    "\n",
    "≤60% = FAIL\n",
    "\n",
    "*>60% = PASS\n",
    "\n",
    "The sample should have more than 60% of the reads from *Streptococcus pneumoniae*.\n",
    "\n",
    "We will run a custom bash script to do the assessment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!bash percentage_calculator.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BONUS!\n",
    "\n",
    "If you are working with BASH in your computer or in a HPC and you have too many files you can optimize commands, loops are very useful for large datasets.\n",
    "\n",
    "Here's a way to do it. \n",
    "\n",
    "Create a new bash script using nano named `taxonomy.sh`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dont run this cell\n",
    "# create a script to run the trimming\n",
    "!nano taxonomy.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then copy and paste the following scrip in your new file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/bin/bash\n",
    "#Author: Nathalia Portilla\n",
    "\n",
    "for i in $(ls *_1.trimmed.fastq.gz); do\n",
    "\n",
    "NAME=$(basename $i _1.trimmed.fastq.gz)\n",
    "echo \"$NAME\"\n",
    "j=\"${NAME}_1.trimmed.fastq.gz\"\n",
    "echo \"$j\"\n",
    "k=\"${NAME}_2.trimmed.fastq.gz\"\n",
    "echo \"$k\"\n",
    "\n",
    "kraken2 --use-names --db k2_standard_08gb_20221209 --paired $j $k --report ${NAME}.kraken.txt --output -\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the file.\n",
    "\n",
    "Finally, you can execute it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dont run this cell\n",
    "!bash taxonomy.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Adapted from:*\n",
    "\n",
    "- Advanced Bioinformatics Course developed for the GPS and JUNO projects - Wellcome Sanger Insitute\n",
    "\n",
    "*Modified by Luisa Sacristán (Universidad de los Andes-CABANA)*"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}