_3_Aligning-ONT-Reads.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Alignment:\n",
    "The good news: Oxford nanopore reads are extrememly long. Repeat regions are no match for long reads, the variable flanking regions around the repeat can determine where in the genome this read belongs. It also means that any ambiguity in alignment is cleared up through other regions in the read.\n",
    "The bad news: Oxford nanopore reads lack accuracy compared to other second generation sequencing technologies.\n",
    "Therefore, seed-and-extend algorithms have had to change tack a bit to allow for more mismatches in the seed before dismissing a read entirely.\n",
    "We will use a bwa mem variant that has the parameters optimised for Oxford Nanopore reads.\n",
    "This tutorial will take you through the basics of converting a set of fastq files to a sam file to a sorted bam file.\n",
    "It is important for bam files to be sorted for many downstream analysis.\n",
    "It's \"sorted\" in the computer way. Bam files are sam files in binary format so it isn't easy to show how the end result.\n",
    "We can use the fastq files that we extracted from our fast5 files and compare the accuracies over the different chemistries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the modules we need.\n",
    "import os\n",
    "import subprocess\n",
    "from Bio import Entrez, SeqIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Set the directories:\n",
    "HOME_DIR = os.environ.get('PORECAMPAU_ANALYSIS_PATH') + \"/\"\n",
    "PORECAMPAU_DATA_DIR = os.environ.get('PORECAMPAU_DATA_PATH') + \"/\"\n",
    "FASTQ_DIRECTORY = PORECAMPAU_DATA_DIR + \"data/fastq/\"\n",
    "ALIGNMENT_DIRECTORY = HOME_DIR + \"alignment/\"\n",
    "\n",
    "# Create the alignment directory if it doesn't already exist.\n",
    "if not os.path.isdir(ALIGNMENT_DIRECTORY):\n",
    "    os.mkdir(ALIGNMENT_DIRECTORY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Download the reference \n",
    "from Bio import Entrez, SeqIO\n",
    "\n",
    "# Use your own email here\n",
    "Entrez.email = \"<insert_email_address_here>\"\n",
    "\n",
    "# Create reference directory and file name.\n",
    "reference_directory = \"/home/researcher/references/\"\n",
    "if not os.path.isdir(reference_directory):\n",
    "    os.mkdir(reference_directory)\n",
    "reference_name = \"Escherichia_coli_k12_MG1655\"\n",
    "reference_file = reference_directory + reference_name + \".fa\"\n",
    "uid = \"U00096.3\"  # This is the uid for E coli genome strain K-12 MG1655\n",
    "handle = Entrez.efetch(db=\"nucleotide\", id=uid, rettype=\"fasta\")\n",
    "fasta_handler = SeqIO.read(handle, \"fasta\")\n",
    "reference_handler = open(reference_file, \"w\")\n",
    "SeqIO.write(fasta_handler, reference_handler, \"fasta\")\n",
    "reference_handler.close()  # Always close the door behind you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Use samtools to generate a fasta index for the reference file.\n",
    "command = \"samtools faidx %s\" % reference_file\n",
    "stderr = subprocess.check_output(command, shell=True, stderr=subprocess.STDOUT)\n",
    "\n",
    "if not stderr == \"\":\n",
    "    print(stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Before we run bwa mem, we will also need to run bwa index on the reference file\n",
    "bwa_index_command = \"bwa index %s\" % reference_file\n",
    "stderr = subprocess.check_output(bwa_index_command, shell=True, stderr=subprocess.STDOUT)\n",
    "\n",
    "if not stderr == \"\":\n",
    "    print(\"Stderr:\\n%s\" % stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Now we can align the genome using bwa-mem\n",
    "# The output is in sam format and printed to the command line.\n",
    "# We can redirect the output using the > sign.\n",
    "SAMPLE_NAME = \"e_coli_R9\"\n",
    "\n",
    "fastq_file = FASTQ_DIRECTORY + \"e_coli_R9/metrichor/pass/2d/2016-11-07_E_COLI_R9_pass.2d.fastq\"\n",
    "sam_file = ALIGNMENT_DIRECTORY + SAMPLE_NAME + \".sam\"\n",
    "bam_file = ALIGNMENT_DIRECTORY + SAMPLE_NAME + \".bam\"\n",
    "sorted_bam_file = ALIGNMENT_DIRECTORY + SAMPLE_NAME + \".sorted.bam\"\n",
    "sorted_bam_file_index = ALIGNMENT_DIRECTORY + SAMPLE_NAME + \".sorted.bai\"\n",
    "\n",
    "bwa_command = \"bwa mem -x ont2d %s %s > %s\" % (reference_file, fastq_file, sam_file)\n",
    "stderr = subprocess.check_output(bwa_command, shell=True, stderr=subprocess.STDOUT)\n",
    "\n",
    "if not stderr == \"\":\n",
    "    print(\"Stderr:\\n%s\" % stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now let's turn that sam file into a bam file.\n",
    "sam_to_bam_command = \"samtools view -b %s -o %s\" % (sam_file, bam_file)\n",
    "stderr = subprocess.check_output(sam_to_bam_command, shell=True, stderr=subprocess.STDOUT)\n",
    "\n",
    "if not stderr == \"\":\n",
    "    print(\"Error: \", stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now sort the sam file\n",
    "sort_bam_command = \"samtools sort -o %s %s\" % (sorted_bam_file, bam_file)\n",
    "stderr = subprocess.check_output(sort_bam_command, shell=True, stderr=subprocess.STDOUT)\n",
    "\n",
    "if not stderr == \"\":\n",
    "    print(\"Error: \", stderr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now index the bam file\n",
    "index_bam_command = \"samtools index %s %s\" % (sorted_bam_file, sorted_bam_file_index)\n",
    "stderr = subprocess.check_output(index_bam_command, shell=True, stderr=subprocess.STDOUT)\n",
    "\n",
    "if not stderr == \"\":\n",
    "    print(\"Error: \", stderr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Yay! We've managed to generate a sorted bam file with an index.\n",
    "Unfortunately I've fooled you into generating one for a subset of reads.\n",
    "For the IGV Jupyter session we will be using the pre-intalled bam files."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}