ngs-41-velvet-genome-assembly.ipynb

{
 "metadata": {
  "name": "ngs-41-velvet-genome-assembly"
 },
 "nbformat": 3,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "source": [
      "# Using Velvet for genome assembly.",
      "",
      "Velvet is a commonly-used genome assembler; below, we'll use it to assemble some E. coli sequences.",
      "",
      "See: http://genome.cshlp.org/content/18/5/821.short"
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "## Requirements",
      "",
      "You'll need to have screed and khmer installed (see the ngs-02-install-screed and ngs-03-install-khmer notebooks).",
      "",
      "Also, mount a copy of the snapshot snap-f69fda89 as /assembly-data: briefly, make a volume in the same zone as your instance; attach it to your instance as (for example) /dev/sdf; then do"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!mkdir /assembly-data",
      "!mount /dev/xvdf /assembly-data"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "First, download and install Velvet."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "cd /mnt"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!curl -O http://www.ebi.ac.uk/~zerbino/velvet/velvet_1.2.06.tgz",
      "!tar xzf velvet_1.2.06.tgz",
      "!cd velvet_1.2.06 && make MAXKMERLENGTH=51 OPENMP=2",
      "!cp velvet_1.2.06/velvet? /usr/local/bin"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Now, copy the raw reads into the ecoli directory."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!mkdir ecoli",
      "!cp /assembly-data/ecoli*.fq.gz ecoli",
      ""
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "cd /mnt/ecoli"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "ls"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "The first command you run, velveth, creates the de Bruijn graph from the reads.  Here is where you have to specify the kinds of reads you're loading in, from where, and the k size for the graph."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!velveth ecoli.31 31 -fastq.gz -shortPaired ecoli_ref-1m.fq.gz"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "The second command, velvetg, explores the graph and builds contigs and scaffolds from it."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!velvetg ecoli.31 -exp_cov auto -cov_cutoff auto"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "And voila, you're done!  The final output file containing your contigs or scaffolds is in the ecoli.31 directory as 'contigs.fa'.  You can use a program included with khmer, 'assemstats3.py', to look at the number of contigs, the total sum of assembly, and the max sequence length in your assembly:"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/sandbox/assemstats3.py 1000 ecoli.31/contigs.fa"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "Let's quickly calculate the contig length distribution, too."
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "import screed",
      "",
      "lengths = []",
      "for read in screed.open('ecoli.31/contigs.fa'):",
      "    lengths.append(len(read.sequence))",
      ""
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "hist(lengths, bins=20)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "The main thing to note here is that there's a strong tendency towards really short contigs.  This is true of ALL short-read assemblers, and it's one reason why picking a lower cutoff is important...",
      "",
      "One neat trick -- use khmer's \"extract-long-sequences\" to pick out only those sequences you think are interesting:"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/sandbox/extract-long-sequences.py 1000 ecoli.31/contigs.fa > ecoli-31-1k.fa"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "These can now be BLASTed or used for further analysis of some sort."
     ]
    },
    {
     "cell_type": "markdown",
     "source": [
      "Another thing that is commonly done is to explore what different parameters -- in particular, what k values -- do to the assembly.  For example, you can run multiple assemblies with multiple different k values:"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "for k in [33, 35]:",
      "    !velveth ecoli.$k $k -fastq.gz -shortPaired ecoli_ref-1m.fq.gz",
      "    !velvetg ecoli.$k -exp_cov auto -cov_cutoff auto"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!python /usr/local/src/khmer/sandbox/assemstats3.py 1000 ecoli.*/contigs.fa"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "So which one of these is \"best\"?  Well, I would argue that the 35 is looking pretty good, based on small number of contigs, biggest sum of bases, and longest contig.  But what happens if you move up to, say, 37?"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      ""
     ],
     "language": "python",
     "outputs": []
    }
   ]
  }
 ]
}