Outline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nanopore DNA sequencing is portable and relatively cheap, allowing real time sequencing in the field.  We see the potential to use nanopore sequencing as an accessible educational experience. With a clear pipeline that Just Works(TM), a citizen scientist could do WIMP (What's in my pot?) analysis on their own samples without the need for any external tools.  Undergrad or high school students could follow the steps of the pipeline to learn about the basics of genome assembly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Why Nanopore is amazing\n",
    "* Long reads, cheap, portable\n",
    "* Comparison to the current standard (Illumina)\n",
    "* Used for detecting ebola/zika?  Sent to space\n",
    "* Sequencing in the jungle (tweet below!)\n",
    "* Idea we care about: you should be able to take a random sample of stuff (ocean water?  dirt?) and sequence it cheaply and easily to find out what's there\n",
    "\n",
    "To see a video of nanopore sequencing in the jungle, click on the below block of code and click the \"Run\" button at the top of this page.  You can see on the left side of the code block that it's been run x number of times (\"In [x]\").  You'll see a star while the code is running (\"In [*]\") and once it's done, you'll be able to see the video."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<blockquote class=\"twitter-tweet\" data-lang=\"en\"><p lang=\"en\" dir=\"ltr\">Welcome to my laboratory :)<br><br>Sequencing long ribosomal cluster from plants, insects &amp; fungi in real-time in the Amazon rainforest. Within a few mins of <a href=\"https://twitter.com/nanopore?ref_src=twsrc%5Etfw\">@nanopore</a> data generated, performed BLAST &amp; got correct hits! Dual indexing looks great for pooling many samples<a href=\"https://twitter.com/hashtag/junglegenomics?src=hash&amp;ref_src=twsrc%5Etfw\">#junglegenomics</a> <a href=\"https://t.co/UQVjYfmU8U\">pic.twitter.com/UQVjYfmU8U</a></p>&mdash; Aaron Pomerantz (@AaronPomerantz) <a href=\"https://twitter.com/AaronPomerantz/status/980873273348038656?ref_src=twsrc%5Etfw\">April 2, 2018</a></blockquote>\n",
       "<script async src=\"https://platform.twitter.com/widgets.js\" charset=\"utf-8\"></script>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%html\n",
    "<blockquote class=\"twitter-tweet\" data-lang=\"en\"><p lang=\"en\" dir=\"ltr\">Welcome to my laboratory :)<br><br>Sequencing long ribosomal cluster from plants, insects &amp; fungi in real-time in the Amazon rainforest. Within a few mins of <a href=\"https://twitter.com/nanopore?ref_src=twsrc%5Etfw\">@nanopore</a> data generated, performed BLAST &amp; got correct hits! Dual indexing looks great for pooling many samples<a href=\"https://twitter.com/hashtag/junglegenomics?src=hash&amp;ref_src=twsrc%5Etfw\">#junglegenomics</a> <a href=\"https://t.co/UQVjYfmU8U\">pic.twitter.com/UQVjYfmU8U</a></p>&mdash; Aaron Pomerantz (@AaronPomerantz) <a href=\"https://twitter.com/AaronPomerantz/status/980873273348038656?ref_src=twsrc%5Etfw\">April 2, 2018</a></blockquote>\n",
    "<script async src=\"https://platform.twitter.com/widgets.js\" charset=\"utf-8\"></script>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How Nanopore actually works\n",
    "Insert a good paragraph description here:\n",
    "\n",
    "My current understanding is just like, DNA strand pulled through a pore, current across pore changes depending on the nucleotide, you end up with an electrical signal in fast5 format that can be converted to a fasta file (could be fastq but no one trusts quality scores?)\n",
    "\n",
    "## Poretools\n",
    "The raw data we get from nanopore sequencing is in the fast5 format.  This is just a series of current values that were read across the pore as the DNA strand passed through it.\n",
    "\n",
    "We are going to start by looking at this fast5 data, containing current values, and converting it to a fasta file that contains nucleotides. \n",
    "\n",
    "This poretools tutorial is adapted from here: http://nbviewer.jupyter.org/github/arq5x/poretools/blob/master/poretools/ipynb/test_run_report.ipynb\n",
    "\n",
    "First we're going to find our fast5 files.  Our sample fast5 file is in the \"data\" folder, so we set the variable dataDirectory to \"data/\".  If you are using your own data, change dataDirectory to the path to your .fast5 files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\r\n"
     ]
    }
   ],
   "source": [
    "# dataDirectory is the path to our fast5 file.\n",
    "# If you are using your own data, change dataDirectory to the path to your .fast5 files.\n",
    "dataDirectory = 'data/'\n",
    "\n",
    "# Print the number of fast5 files in the dataDirectory.\n",
    "# Click the \"Run\" button at the top of this page to run this code.\n",
    "!find $dataDirectory -maxdepth 1 -name \"*.fast5\" | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poretools has a number of command line options.  Let's start with the stats command, which will give us some basic statistics about our reads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "total reads\t6\n",
      "total base pairs\t25217\n",
      "mean\t4202.83\n",
      "median\t4205\n",
      "min\t2940\n",
      "max\t5826\n",
      "N25\t5079\n",
      "N50\t5011\n",
      "N75\t3399\n"
     ]
    }
   ],
   "source": [
    "# The -q option stops poretools outputting any warning messages.\n",
    "!poretools stats -q $dataDirectory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our sample data has 6 reads and 25,217 base pairs. (Anything else of interest to say about this info?)\n",
    "\n",
    "We have 3 reads per fast5 because forward, reverse, and two-directional reads are all counted separately. (Is this correct?) We can see the stats about each of these types of reads using the below code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "total reads\t2\n",
      "total base pairs\t8019\n",
      "mean\t4009.50\n",
      "median\t4009\n",
      "min\t2940\n",
      "max\t5079\n",
      "N25\t5079\n",
      "N50\t5079\n",
      "N75\t2940\n"
     ]
    }
   ],
   "source": [
    "# Look at stats for forward strands\n",
    "!poretools stats -q --type fwd $dataDirectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "total reads\t2\n",
      "total base pairs\t7973\n",
      "mean\t3986.50\n",
      "median\t3986\n",
      "min\t2962\n",
      "max\t5011\n",
      "N25\t5011\n",
      "N50\t5011\n",
      "N75\t2962\n"
     ]
    }
   ],
   "source": [
    "# Look at stats for reverse strands\n",
    "!poretools stats -q --type rev $dataDirectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "total reads\t2\n",
      "total base pairs\t9225\n",
      "mean\t4612.50\n",
      "median\t4612\n",
      "min\t3399\n",
      "max\t5826\n",
      "N25\t5826\n",
      "N50\t5826\n",
      "N75\t3399\n"
     ]
    }
   ],
   "source": [
    "# Look at two-directional reads\n",
    "!poretools stats -q --type 2D $dataDirectory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully we are going to add how to make a squiggle plot here at some point, that shows the current changing and gives a good idea of what signal actually looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add squiggle plot here!!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to convert our fast5 file to fasta.  Fasta is a common format for storing DNA sequences.  The below code will take each of the fast5 files in dataDirectory, create a fasta file of that sequence, and store it in a folder called fastaOutput."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\r\n",
      "  from ._conv import register_converters as _register_converters\r\n"
     ]
    }
   ],
   "source": [
    "# Make a folder to store our fasta files in.\n",
    "!mkdir fastaOutput\n",
    "\n",
    "# Convert our fast5 files to fasta.\n",
    "!poretools fasta $dataDirectory > fastaOutput/outputPoretoolData.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look at the first few lines of this fasta file to see what's in it.  Each of the sequences has a line containing \">\" and then a unique identfier, followed by a line containing the nucleotide sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">b233f432-7786-4b0b-8b2d-03c2e168a45b_Basecall_2D_2d CPHG_CNU4299G4G_louse_library_2016_3_4_3507_1_ch120_read240_strand data/2016_3_4_3507_1_ch120_read240_strand.fast5\r\n",
      "GAAATTGCTCCCGCTCTCAGTTCTGCTTTAACAGATAAATTAATAACATATCAATAAAGCATCAAAATCACGTGATTGGAACGCCGTACTTCGAAGAGGAGGATGGAGACGAGGATGGGAGCAGAGGGGAGGATGTGCACTTCTCCCCACGTCAGTTGGGATTCGAAGGAAGTTTGCGGCTTGTTTTAGAGTGGAGGACA\r\n"
     ]
    }
   ],
   "source": [
    "# This will show us the first 200 characters of the first two lines of our file.\n",
    "# We don't want to look at the whole sequence because it's going to be really long!\n",
    "!cut -c -200 fastaOutput/outputPoretoolData.fasta | head -2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This fasta file containing our raw reads will be the input to the next steps in our pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\r\n",
      "  from ._conv import register_converters as _register_converters\r\n"
     ]
    }
   ],
   "source": [
    "!poretools winner $dataDirectory > winner.fasta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "ImportError",
     "evalue": "No module named 'PIL'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-6c308242257f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mPIL\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mImageFile\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msystem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"head -1 winner.fasta | sed 's/>.* //' | xargs poretools squiggle --saveas png --num-facets 12\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0msquiggle\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetoutput\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"head -1 winner.fasta | sed 's/>.* //'\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0msquiggle\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msquiggle\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'.png'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mImageFile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msquiggle\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mImportError\u001b[0m: No module named 'PIL'"
     ]
    }
   ],
   "source": [
    "# Alignment and assembly\n",
    "\n",
    "(I don't really know how to explain this so I'll put Brandon's explanation in here)\n",
    "\n",
    "Minimap2, we take our fasta file of raw reads and we try to align all of the reads to each other (all vs. all).  In order to do this we tell minimap2 to compare our fasta file to itself.  We write this to another file that contains the (alignment?) called overlap.paf.\n",
    "\n",
    "Because real fast5 files are humongous (very high frequency resolution), we didn't use a real fast5 sample for poretools.  For this step, we want real sample data, so we're actually going to use a different fasta file.  We're going to wget this and then run minimap2 on it.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-04-04 21:36:25--  https://nanopore.s3.climb.ac.uk/MAP006-1_2D_pass.fasta\n",
      "Resolving nanopore.s3.climb.ac.uk (nanopore.s3.climb.ac.uk)... 137.205.52.6\n",
      "Connecting to nanopore.s3.climb.ac.uk (nanopore.s3.climb.ac.uk)|137.205.52.6|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 255168226 (243M) [text/plain]\n",
      "Saving to: ‘MAP006-1_2D_pass.fasta.2’\n",
      "\n",
      "MAP006-1_2D_pass.fa 100%[===================>] 243.35M  5.64MB/s    in 55s     \n",
      "\n",
      "2018-04-04 21:37:21 (4.43 MB/s) - ‘MAP006-1_2D_pass.fasta.2’ saved [255168226/255168226]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://nanopore.s3.climb.ac.uk/MAP006-1_2D_pass.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're going to tell minimap2 to try to align these reads to each other and write the overlaps to a file called overlap.paf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[M::mm_idx_gen::14.047*0.86] collected minimizers\n",
      "Killed\n"
     ]
    }
   ],
   "source": [
    "!minimap2 -x ava-ont -t 1 MAP006-1_2D_pass.fasta MAP006-1_2D_pass.fasta.2 > overlap.paf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Miniasm, we take these overlapping reads and assemble them into contigs.  Need the raw reads from poretools and the overlap file from minimap2 (I don't really understand what this means so once again gonna use a paragraph from brandon).  (because this outputs a .gfa file, we're piping the output into awk and using that to convert it to fasta format, then writing it to asm.fa)\n",
    "\n",
    "A contig is a combination of several shorter reads that overlap to make one continuous sequence.  This \"contiguous\" segment of DNA is thought to be a true representation of a piece of the organism's genome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[M::main] ===> Step 1: reading read mappings <===\r\n",
      "[M::ma_hit_read::0.003*0.00] read 3 hits; stored 6 hits and 5 sequences (22277 bp)\r\n",
      "[M::main] ===> Step 2: 1-pass (crude) read selection <===\r\n",
      "[M::ma_hit_sub::0.003*0.00] 0 query sequences remain after sub\r\n",
      "[M::ma_hit_cut::0.003*0.00] 0 hits remain after cut\r\n",
      "[M::ma_hit_flt::0.003*0.00] 0 hits remain after filtering; crude coverage after filtering: -nan\r\n",
      "[M::main] ===> Step 3: 2-pass (fine) read selection <===\r\n",
      "[M::ma_hit_sub::0.003*0.00] 0 query sequences remain after sub\r\n",
      "[M::ma_hit_cut::0.003*0.00] 0 hits remain after cut\r\n",
      "[M::ma_hit_contained::0.004*0.00] 0 sequences and 0 hits remain after containment removal\r\n",
      "[M::main] ===> Step 4: graph cleaning <===\r\n",
      "[M::ma_sg_gen] read 0 arcs\r\n",
      "[M::main] ===> Step 4.1: transitive reduction <===\r\n",
      "[M::asg_arc_del_trans] transitively reduced 0 arcs\r\n",
      "[M::main] ===> Step 4.2: initial tip cutting and bubble popping <===\r\n",
      "[M::asg_cut_tip] cut 0 tips\r\n",
      "[M::asg_arc_del_multi] removed 0 multi-arcs\r\n",
      "[M::asg_arc_del_asymm] removed 0 asymmetric arcs\r\n",
      "[M::asg_pop_bubble] popped 0 bubbles and trimmed 0 tips\r\n",
      "[M::main] ===> Step 4.3: cutting short overlaps (3 rounds in total) <===\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::main] ===> Step 4.4: removing short internal sequences and bi-loops <===\r\n",
      "[M::asg_cut_internal] cut 0 internal sequences\r\n",
      "[M::asg_cut_biloop] cut 0 small bi-loops\r\n",
      "[M::asg_cut_tip] cut 0 tips\r\n",
      "[M::asg_pop_bubble] popped 0 bubbles and trimmed 0 tips\r\n",
      "[M::main] ===> Step 4.5: aggressively cutting short overlaps <===\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::main] ===> Step 5: generating unitigs <===\r\n",
      "[M::main] Version: 0.2-r168-dirty\r\n",
      "[M::main] CMD: miniasm -f fastaOutput/outputPoretoolData.fasta overlap.paf\r\n",
      "[M::main] Real time: 0.008 sec; CPU: 0.000 sec\r\n"
     ]
    }
   ],
   "source": [
    "!miniasm -f fastaOutput/outputPoretoolData.fasta overlap.paf | awk '/^S/{print \">\"$2\"\\n\"$3}' | fold > asm.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alignment and assembly\n",
    "\n",
    "(I don't really know how to explain this so I'll put Brandon's explanation in here)\n",
    "\n",
    "Minimap2, we take our fasta file of raw reads and we try to align all of the reads to each other (all vs. all).  In order to do this we tell minimap2 to compare our fasta file to itself.  We write this to another file that contains the (alignment?) called overlap.paf.\n",
    "\n",
    "Because real fast5 files are humongous (very high frequency resolution), we didn't use a real fast5 sample for poretools.  For this step, we want real sample data, so we're actually going to use a different fasta file.  We're going to wget this and then run minimap2 on it.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2018-04-04 21:36:25--  https://nanopore.s3.climb.ac.uk/MAP006-1_2D_pass.fasta\n",
      "Resolving nanopore.s3.climb.ac.uk (nanopore.s3.climb.ac.uk)... 137.205.52.6\n",
      "Connecting to nanopore.s3.climb.ac.uk (nanopore.s3.climb.ac.uk)|137.205.52.6|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 255168226 (243M) [text/plain]\n",
      "Saving to: ‘MAP006-1_2D_pass.fasta.2’\n",
      "\n",
      "MAP006-1_2D_pass.fa 100%[===================>] 243.35M  5.64MB/s    in 55s     \n",
      "\n",
      "2018-04-04 21:37:21 (4.43 MB/s) - ‘MAP006-1_2D_pass.fasta.2’ saved [255168226/255168226]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://nanopore.s3.climb.ac.uk/MAP006-1_2D_pass.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're going to tell minimap2 to try to align these reads to each other and write the overlaps to a file called overlap.paf."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[M::mm_idx_gen::14.047*0.86] collected minimizers\n",
      "Killed\n"
     ]
    }
   ],
   "source": [
    "!minimap2 -x ava-ont -t 1 MAP006-1_2D_pass.fasta MAP006-1_2D_pass.fasta.2 > overlap.paf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Miniasm, we take these overlapping reads and assemble them into contigs.  Need the raw reads from poretools and the overlap file from minimap2 (I don't really understand what this means so once again gonna use a paragraph from brandon).  (because this outputs a .gfa file, we're piping the output into awk and using that to convert it to fasta format, then writing it to asm.fa)\n",
    "\n",
    "A contig is a combination of several shorter reads that overlap to make one continuous sequence.  This \"contiguous\" segment of DNA is thought to be a true representation of a piece of the organism's genome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[M::main] ===> Step 1: reading read mappings <===\r\n",
      "[M::ma_hit_read::0.003*0.00] read 3 hits; stored 6 hits and 5 sequences (22277 bp)\r\n",
      "[M::main] ===> Step 2: 1-pass (crude) read selection <===\r\n",
      "[M::ma_hit_sub::0.003*0.00] 0 query sequences remain after sub\r\n",
      "[M::ma_hit_cut::0.003*0.00] 0 hits remain after cut\r\n",
      "[M::ma_hit_flt::0.003*0.00] 0 hits remain after filtering; crude coverage after filtering: -nan\r\n",
      "[M::main] ===> Step 3: 2-pass (fine) read selection <===\r\n",
      "[M::ma_hit_sub::0.003*0.00] 0 query sequences remain after sub\r\n",
      "[M::ma_hit_cut::0.003*0.00] 0 hits remain after cut\r\n",
      "[M::ma_hit_contained::0.004*0.00] 0 sequences and 0 hits remain after containment removal\r\n",
      "[M::main] ===> Step 4: graph cleaning <===\r\n",
      "[M::ma_sg_gen] read 0 arcs\r\n",
      "[M::main] ===> Step 4.1: transitive reduction <===\r\n",
      "[M::asg_arc_del_trans] transitively reduced 0 arcs\r\n",
      "[M::main] ===> Step 4.2: initial tip cutting and bubble popping <===\r\n",
      "[M::asg_cut_tip] cut 0 tips\r\n",
      "[M::asg_arc_del_multi] removed 0 multi-arcs\r\n",
      "[M::asg_arc_del_asymm] removed 0 asymmetric arcs\r\n",
      "[M::asg_pop_bubble] popped 0 bubbles and trimmed 0 tips\r\n",
      "[M::main] ===> Step 4.3: cutting short overlaps (3 rounds in total) <===\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::main] ===> Step 4.4: removing short internal sequences and bi-loops <===\r\n",
      "[M::asg_cut_internal] cut 0 internal sequences\r\n",
      "[M::asg_cut_biloop] cut 0 small bi-loops\r\n",
      "[M::asg_cut_tip] cut 0 tips\r\n",
      "[M::asg_pop_bubble] popped 0 bubbles and trimmed 0 tips\r\n",
      "[M::main] ===> Step 4.5: aggressively cutting short overlaps <===\r\n",
      "[M::asg_arc_del_short] removed 0 short overlaps\r\n",
      "[M::main] ===> Step 5: generating unitigs <===\r\n",
      "[M::main] Version: 0.2-r168-dirty\r\n",
      "[M::main] CMD: miniasm -f fastaOutput/outputPoretoolData.fasta overlap.paf\r\n",
      "[M::main] Real time: 0.008 sec; CPU: 0.000 sec\r\n"
     ]
    }
   ],
   "source": [
    "# Retrieve genome sequences for each species\n",
    "!wget ftp://ftp.ensemblgenomes.org/pub/bacteria/release-33/fasta/bacteria_9_collection/escherichia_coli_0_1304/dna/Escherichia_coli_0_1304.GCA_000303235_1.dna.nonchromosomal.fa.gz\n",
    "!wget ftp://ftp.ensemblgenomes.org/pub/release-38/bacteria//fasta/bacteria_15_collection/_clostridium_asparagiforme_dsm_15981/dna/_clostridium_asparagiforme_dsm_15981.ASM15807v1.dna.nonchromosomal.fa.gz\n",
    "!wget ftp://ftp.ensemblgenomes.org/pub/release-38/bacteria//fasta/bacteria_176_collection/_bacillus_aminovorans/dna/_bacillus_aminovorans.ASM164324v1.dna.nonchromosomal.fa.gz\n",
    "!wget https://downloads.yeastgenome.org/sequence/S288C_reference/orf_dna/orf_genomic_1000_all.fasta.gz \n",
    "# Give yeast file a better name\n",
    "!mv  orf_genomic_1000_all.fasta.gz yeast_genomic_1000_all.fa.gz "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should be able to see these four files in your directory now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_bacillus_aminovorans.ASM164324v1.dna.nonchromosomal.fa.gz\r\n",
      "_clostridium_asparagiforme_dsm_15981.ASM15807v1.dna.nonchromosomal.fa.gz\r\n",
      "Escherichia_coli_0_1304.GCA_000303235_1.dna.nonchromosomal.fa.gz\r\n",
      "yeast_genomic_1000_all.fa.gz\r\n"
     ]
    }
   ],
   "source": [
    "# See our four FASTA genome files\n",
    "ls *.fa.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many techniques to identify which species you have sequenced. One approach would be to try to align each of our contigs to each reference in our database. However, that would be extremely computationally intensive. Instead here we will use a computationally light application called MASH, which calculates an estimate for the distance between our query sequence and each reference.  In this context, distance refers to how similar or different two genomes are.\n",
    "\n",
    "First we must form a sketch from each of our reference genomes. A sketch is a minified representation of a reference. This contains less information than a whole genome sequence, but still allows us to calculate the distance between each of the references and our query.\n",
    "\n",
    "We have already downloaded the sequences for each sample species described above. Now we must build the sketch of each reference using the command sketch from the mash tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sketching /work/data/Escherichia_coli_0_1304.GCA_000303235_1.dna.nonchromosomal.fa.gz...\n",
      "Sketching /work/data/_bacillus_aminovorans.ASM164324v1.dna.nonchromosomal.fa.gz...\n",
      "Sketching /work/data/_clostridium_asparagiforme_dsm_15981.ASM15807v1.dna.nonchromosomal.fa.gz...\n",
      "Sketching /work/data/yeast_genomic_1000_all.fa.gz...\n",
      "Writing to reference.msh...\n"
     ]
    }
   ],
   "source": [
    "# The -o option allows us to name the sketch output file.\n",
    "!mash sketch /work/data/*.fa.gz -o reference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have our database of references in one file, reference.msh.  We want to take our assembly from before and compare it to the reference database.\n",
    "\n",
    "When we assembled our contigs using miniasm, it tried to piece together all of the individual strands of DNA sequenced by the Nanopore.  If more than one species was present or the assembler couldn't figure out where to place small pieces of the organism's genome, then the assembly will be made of multiple contigs.  If we sequenced human DNA, we would naturally get many contigs, one for each chromosome.  In practice, it's hard to piece together the whole genome so even within one chromosome some gaps may remain unfilled.\n",
    "\n",
    "To build a sketch for each contig we must us the -i option of the sketch command. This option tells mash to treat each contig as its own seperate entity. Here we again use the -o option as described above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sketching asm.fa...\r\n",
      "Writing to contigs.msh...\r\n"
     ]
    }
   ],
   "source": [
    "# Note: if there's no star on the left side of the code block (\"In [*]:\") then you know it's done running.\n",
    "!mash sketch -i asm.fa -o contigs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are prepared to compute the distances between each contig and each reference. This is accomplished using the dist command. The '>' signifies to take the output from the following program and create a file called distance.tab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mash dist reference.msh contigs.msh > distance.tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the head command to look at the first five lines of this distance file. This file is organized into the following columns: [reference-ID, query-ID, distance, p-value, shared-hashes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head distance.tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}