notebooks/read-spectrum-lib.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:d32b83417b6968946bb27448a5eb50cc410af80910308cbb56f6b638db279856"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "How to read a PRIDE Cluster spectrum library into a Spark RDD"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this notebook we will explain how to read a [PRIDE Cluster](http://wwwdev.ebi.ac.uk/pride/cluster/#/libraries) spectrum library into a Spark RDD. This will allow us to perform further operations with spectral data into a Spark cluster.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Splitting a library into individual MGF files"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In order to better read our RDD, we need individual MGF files containing a single spectra per file. Right now we have a library with multiple spectra per file. We will split it using normal Python file I/O.   "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we define our split function, that basically reads the source file secuantially. It transfers every line to a destination file. Every time an empty line appears, the destination file is closed and a new one is created with a different name.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import gzip \n",
      "\n",
      "def split_library(lib_path, destination_folder):\n",
      "    # create destination folder if needed\n",
      "    if not os.path.exists(destination_folder):\n",
      "        os.mkdir(destination_folder)\n",
      "        \n",
      "    # open library file for reading\n",
      "    input_f = gzip.open(lib_path, 'rb')\n",
      "    \n",
      "    # open destination file for writing\n",
      "    i = 0\n",
      "    target_f = open(os.path.join(destination_folder,str(i)+\".mgf\"), 'w')\n",
      "    \n",
      "    # process input library file\n",
      "    for line in input_f.readlines():\n",
      "        if line==\"\\n\": # end of spectrum, move to next file\n",
      "            target_f.close()\n",
      "            i+=1\n",
      "            target_f = open(os.path.join(destination_folder,str(i)+\".mgf\"), 'w')\n",
      "        else: # write line\n",
      "            target_f.write(line)\n",
      "\n",
      "    target_f.close()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can use the split function to split our human spectrum library.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import urllib\n",
      "\n",
      "human_ftp = \"ftp://ftp.pride.ebi.ac.uk/pride/data/cluster/spectrum-libraries/1.0.1/Human.msp.gz\"\n",
      "f = urllib.urlretrieve (human_ftp, \"../Human.msp.gz\")\n",
      "\n",
      "split_library(\"../Human.msp.gz\", \"../human\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Read multiple MGF files into an RDD"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We want to read multiple MGF files into a single RDD of key-value pairs where the key is the name of the peptide, and the value is a list of peaks (pairs of numbers). Having spectra in this shape will allow us later on to perform distance calculation in parallel using our Spark cluster.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will start by using `SparkContext.wholeTextFiles`, that lets you read a directory containing multiple small text files, and returns each of them as (filename, content) pairs. The loading process takes a while, depending on the size of the library itself, but it is a long process anyway. The current human library contains a lot of spectra.   "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"There are {} spectrum files in the folder\".format(len([name for name in os.listdir('../human') if os.path.isfile(os.path.join('../human', name))]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "There are 80593 spectrum files in the folder\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's read them into our RDD.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from time import time\n",
      "\n",
      "raw_file = sc.wholeTextFiles (\"../human\")\n",
      "t0 = time()\n",
      "raw_file.take(2)\n",
      "tt = time() - t0\n",
      "\n",
      "print \"Raw data read in {} seconds\".format(round(tt,3))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Raw data read in 441.324 seconds\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We need to to parse those entries using `map` into the right format. This is the parse function we are going to apply to each RDD element.    "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def parse_file_entry(file_entry):\n",
      "    # each file entry is a tuple with the first element being the file path and the second its contents\n",
      "    file_split = file_entry[1].split(\"\\n\")\n",
      "    # name is the first element in the split, we need to remove the 'name: ' prefix\n",
      "    name = file_split[0].split(\" \")[1]\n",
      "    # number of peaks is the third element, removed 'Num peaks: ' prefix\n",
      "    num_peaks = int(file_split[2].split(\" \")[2])\n",
      "    # read the peaks themselves, that are in the next num_peaks lines, space separated (there is an extra space after each peak)\n",
      "    peak_list = [tuple(map(float, peak.split(\" \")[0:2])) for peak in file_split[3:3+num_peaks]]\n",
      "    \n",
      "    return (name, peak_list)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So we proceed to apply it using `map`.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "human_spectrum_library = raw_file.map(parse_file_entry)\n",
      "t0 = time()\n",
      "human_spectrum_library.take(2)\n",
      "tt = time() - t0\n",
      "\n",
      "print \"Raw data parsed in {} seconds\".format(round(tt,3))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Raw data parsed in 4.59 seconds\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, we have our RDD containing the whole spectrum library for Human, ready to use.  "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Saving the library"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We want to reuse the library later on, withouth having to read and parse the MGF files again. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "human_spectrum_library.saveAsPickleFile(\"../human/lib.file\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    }
   ],
   "metadata": {}
  }
 ]
}