ngs-31-python-and-short-reads.ipynb

{
 "metadata": {
  "name": "ngs-31-python-and-short-reads"
 },
 "nbformat": 3,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "input": [
      "!df",
      "!mount /dev/xvdf /data"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "cd /mnt"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!head /data/drosophila/RAL357_1.fastq"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "source": [
      "# Introducing screed",
      "",
      "screed is a Python library built in the Brown Lab to work with sequences.  It includes both functionality for reading in sequences in FASTA and FASTQ, and creating a random-access sequence database for efficiently retrieving specific named sequences.",
      "",
      "One of the most convenient things about screed is that it will automatically figure out whether or not it's a FASTA or FASTQ file, and automatically read gz and bz2 files.",
      "",
      "Read more about it here:",
      "",
      "http://readthedocs.org/docs/screed/en/latest/"
     ]
    },
    {
     "cell_type": "code",
     "input": [
      "import screed",
      "",
      "for read in screed.open('/data/drosophila/RAL357_1.fastq'):",
      "    print read",
      "    break",
      ""
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "for read in screed.open('/data/drosophila/RAL357_1.fastq'):",
      "    print read.name, read.sequence",
      "    break"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "for n, read in enumerate(screed.open('/data/drosophila/RAL357_1.fastq')):",
      "    if n > 5:",
      "        break",
      "    print read.name, read.sequence"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "lengths = []",
      "for n, read in enumerate(screed.open('/data/drosophila/RAL357_1.fastq')):",
      "    if n > 10000: break",
      "        ",
      "    seq = read.sequence",
      "    first_n = seq.find('N')  # or you could do seq = seq.rstrip('N')",
      "    if first_n >= 0:",
      "        seq = seq[0:first_n]",
      "        ",
      "    read_len = len(seq)",
      "    lengths.append(read_len)",
      ""
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "print seq",
      "print seq[:5]"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "hist(lengths, cumulative=True)",
      "axis(ymax=100, xmin=0)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "nlocs = []",
      "for n, read in enumerate(screed.open('/data/drosophila/RAL357_1.fastq')):",
      "    if n > 10000: break",
      "        ",
      "    seq = read.sequence",
      "    first_n = seq.find('N')  # or you could do seq = seq.rstrip('N')",
      "    if first_n >= 0:",
      "        nlocs.append(first_n)",
      ""
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "hist(nlocs, cumulative=True, bins=100)",
      "axis(ymax=100, xmin=0)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "# pick off the first 100,000 sequences from the FASTQ file and align them.",
      "!curl -O ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_r5.37_FB2011_05/fasta/dmel-all-chromosome-r5.37.fasta.gz",
      "!gunzip -f dmel-all-chromosome-r5.37.fasta.gz",
      "!bowtie-index dmel-all-chromosome-r5.37.fasta drosophila_bowtie",
      "",
      "!head -400000 /data/drosophila/RAL357_1.fastq > RAL357_1.100k.fastq",
      "!bowtie -p 2 -q drosophila_bowtie RAL357_1.100k.fastq > ral357-random.map"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "!head -1 ral357-random.map"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "# define a function that we can use to read in the mapping output file",
      "def read_mapping_file(filename):",
      "    for line in open(filename):",
      "        _, orient, refname, pos, read, _, _, mismatches = line.split('\\t')[:8]",
      "        pos = int(pos)",
      "        mismatches = mismatches.split(',')",
      "        x = []",
      "        for posn in mismatches:",
      "            posn = posn.split(':')[0]",
      "            posn = int(posn)",
      "            x.append(posn)",
      "        yield orient, refname, pos, read, x"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "d = {}",
      "",
      "for (orient, refname, pos, read, mismatches) in read_mapping_file('ral357-random.map'):",
      "    for posn in mismatches:",
      "        d[posn] = d.get(posn, 0) + 1"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "items = d.items()",
      "items.sort()",
      "",
      "x = []",
      "y = []",
      "for (a, b) in items:",
      "    x.append(a)",
      "    y.append(b)",
      ""
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      "plot(x, y)",
      "axis(ymax=4000)"
     ],
     "language": "python",
     "outputs": []
    },
    {
     "cell_type": "code",
     "input": [
      ""
     ],
     "language": "python",
     "outputs": []
    }
   ]
  }
 ]
}