{
 "metadata": {
  "name": "(Galaxy) Example 1. finding coding exon which has the highest number of single nucleotide polymorphisms on chromosome 22"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "This is an example from the galaxy project in IPython notebook. https://main.g2.bx.psu.edu/u/aun1/p/galaxy101"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "1. Getting data from UCSC\n==========================\n\n 1.0. Getting coding exons\n--------------------------\n\nFirst thing we will do is to obtain data from UCSC"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from IPython.core.display import HTML\nHTML(\"<iframe src='http://genome.ucsc.edu/cgi-bin/hgTables?GALAXY_URL=http%3A//galaxyproject.cloudapp.net%3A8080/tool_runner&tool_id=ucsc_table_direct1&hgta_compressType=none&sendToGalaxy=1&hgta_outputType=bed' width=880 height=600>\")",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": "<iframe src='http://genome.ucsc.edu/cgi-bin/hgTables?GALAXY_URL=http%3A//galaxyproject.cloudapp.net%3A8080/tool_runner&tool_id=ucsc_table_direct1&hgta_compressType=none&sendToGalaxy=1&hgta_outputType=bed' width=880 height=600>",
       "output_type": "pyout",
       "prompt_number": 1,
       "text": "<IPython.core.display.HTML at 0x18701d0>"
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Make sure that your settings are exactly the same as shown on the screen (in particular, position should be set to \"chr22\", output format should be set to \"BED - browser extensible data\", and \"Galaxy\" should be checked by Send output to option). Click get output and you will see the next screen:"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "here make sure Create one BED record per is set to \"Coding Exons\" and click Send Query to Galaxy. After this you will see your first History Item in Galaxy's right pane. It will go through gray (preparing) and yellow (running) states to become green:\n"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "The first step will download Exons on a local directory. The filepath will be used as an input to a next pipeline.\nIn this example, first_query is used to indicate Exons location."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "first_query = \"/home/azureuser/galaxy-dist/database/files/000/dataset_37.dat\"",
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "1.1. Getting SNPs\n------------------\n\nNow is the time to obtain SNP data. This is done almost exactly the same way. First thing we will do is to again click on \"Get Data -> UCSC Main\":"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "from IPython.core.display import HTML\nHTML(\"<iframe src='http://genome.ucsc.edu/cgi-bin/hgTables?GALAXY_URL=http%3A//galaxyproject.cloudapp.net%3A8080/tool_runner&tool_id=ucsc_table_direct1&hgta_compressType=none&sendToGalaxy=1&hgta_outputType=bed' width=880 height=600>\")",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": "<iframe src='http://genome.ucsc.edu/cgi-bin/hgTables?GALAXY_URL=http%3A//galaxyproject.cloudapp.net%3A8080/tool_runner&tool_id=ucsc_table_direct1&hgta_compressType=none&sendToGalaxy=1&hgta_outputType=bed' width=880 height=600>",
       "output_type": "pyout",
       "prompt_number": 2,
       "text": "<IPython.core.display.HTML at 0x1931950>"
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "where you need to make sure that Whole Gene is selected (\"Whole Gene\" here really means \"Whole Feature\") and click Send Query to Galaxy. You will get your second item in the history:\n"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Now we will rename the two history items to \"Exons\" and \"SNPs\" by clicking on the Pencil icon adjacent to each item. Also we will rename history to \"Galaxy 101\" (or whatever you want) by clicking on \"Unnamed history\" so everything looks like this:\n"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "The second step downloads SNPs to a local directory as well. second_query is used to indicate the path of the data."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "second_query = \"/home/azureuser/galaxy-dist/database/files/000/dataset_38.dat\"",
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "2. Finding Exons with the highest number of SNPs\n================================================\n\n2.0. Joining exons with SNPs\n------------------------------\n\nLet's remind ourselves that our objective was to find which exon contains the most SNPs. This first step in answering this question will be joining exons with SNPs (a fancy word for printing exons and SNPs that overlap side by side). This is done using \"Operate on Genomics Intervals -> Join\" tool:"
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Instead of selecting input values on the web, python variables used to indicate the downloaded Exons and SNPs."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "first_query = \"/home/azureuser/galaxy-dist/database/files/000/dataset_37.dat\"\nsecond_query = \"/home/azureuser/galaxy-dist/database/files/000/dataset_38.dat\"\noutput = \"/home/azureuser/galaxy-dist/database/files/000/dataset_39.dat\"",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Galaxy libraries will be loaded in main python path."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import os, sys, logging\ngalaxy_dir = \"/home/azureuser/galaxy-dist/\"\nnew_path = [ os.path.join( galaxy_dir, \"lib\" ) ]\nnew_path.extend( sys.path[1:] ) # remove scripts/ from the path\nsys.path = new_path",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Next part is a python code for joing exons with SNPs. Original script is gops_join.py under galaxy-dist/tools/new_operations"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#!/usr/bin/env python\n\"\"\"\nJoin two sets of intervals using their overlap as the key.\n\nusage: %prog bed_file_1 bed_file_2 out_file\n    -1, --cols1=N,N,N,N: Columns for start, end, strand in first file\n    -2, --cols2=N,N,N,N: Columns for start, end, strand in second file\n    -m, --mincols=N: Require this much overlap (default 1bp)\n    -f, --fill=N: none, right, left, both\n\"\"\"\nfrom galaxy import eggs\nimport pkg_resources\npkg_resources.require( \"bx-python\" )\nimport sys, traceback, fileinput\nfrom warnings import warn\nfrom bx.intervals import *\nfrom bx.intervals.io import *\nfrom bx.intervals.operations.join import *\nfrom bx.cookbook import doc_optparse\nfrom galaxy.tools.util.galaxyops import *\n\nassert sys.version_info[:2] >= ( 2, 4 )\n\n\nmincols = 1\nupstream_pad = 0\ndownstream_pad = 0\nleftfill = False\nrightfill = False\n\n#options, args = doc_optparse.parse( __doc__ )\n\n\ntry:\n    #chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )\n    #chr_col_2, start_col_2, end_col_2, strand_col_2 = parse_cols_arg( options.cols2 )\n    #if options.mincols: mincols = int( options.mincols )\n    #if options.fill:\n    #    if options.fill == \"both\":\n    #        rightfill = leftfill = True\n    #    else:\n    #        rightfill = options.fill == \"right\"\n    #        leftfill = options.fill == \"left\"\n    chr_col_1, start_col_1, end_col_1, strand_col_1 = 0,1,2,5\n    chr_col_2, start_col_2, end_col_2, strand_col_2 = 0,1,2,5\n    mincols = 1\n    args = [ first_query, second_query, output ]\n    in_fname, in2_fname, out_fname = args\nexcept:\n    doc_optparse.exception()\n\ng1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),\n                        chrom_col=chr_col_1,\n                        start_col=start_col_1,\n                        end_col=end_col_1,\n                        strand_col=strand_col_1,\n                        fix_strand=True )\ng2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),\n                        chrom_col=chr_col_2,\n                        start_col=start_col_2,\n                        end_col=end_col_2,\n                        strand_col=strand_col_2,\n                        fix_strand=True )\n\nout_file = open( out_fname, \"w\" )\n\ntry:\n    for outfields in join(g1, g2, mincols=mincols, rightfill=rightfill, leftfill=leftfill):\n        if type( outfields ) is list:\n            out_file.write( \"%s\\n\" % \"\\t\".join( outfields ) )\n        else:\n            out_file.write( \"%s\\n\" % outfields )\nexcept ParseError, exc:\n    out_file.close()\n    fail( \"Invalid file format: %s\" % str( exc ) )\nexcept MemoryError:\n    out_file.close()\n    fail( \"Input datasets were too large to complete the join operation.\" )\n\nout_file.close()\n\nif g1.skipped > 0:\n    print skipped( g1, filedesc=\" of 1st dataset\" )\nif g2.skipped > 0:\n    print skipped( g2, filedesc=\" of 2nd dataset\" )\n",
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "!head $output",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "chr22\t16258185\t16258303\tuc002zlh.1_cds_1_0_chr22_16258186_r\t0\t-\tchr22\t16258278\t16258279\trs2845178\t0\t+\r\nchr22\t16266928\t16267095\tuc002zlh.1_cds_2_0_chr22_16266929_r\t0\t-\tchr22\t16267031\t16267032\trs7292200\t0\t+\r\nchr22\t16266928\t16267095\tuc002zlh.1_cds_2_0_chr22_16266929_r\t0\t-\tchr22\t16267011\t16267012\trs7290262\t0\t+\r\nchr22\t16266928\t16267095\tuc002zlh.1_cds_2_0_chr22_16266929_r\t0\t-\tchr22\t16266963\t16266964\trs10154680\t0\t+\r\nchr22\t16266928\t16267095\tuc002zlh.1_cds_2_0_chr22_16266929_r\t0\t-\tchr22\t16267037\t16267038\trs2818572\t0\t+\r\nchr22\t16269872\t16269943\tuc002zlh.1_cds_4_0_chr22_16269873_r\t0\t-\tchr22\t16269933\t16269934\trs2845206\t0\t+\r\nchr22\t16275206\t16275277\tuc002zlh.1_cds_5_0_chr22_16275207_r\t0\t-\tchr22\t16275252\t16275253\trs8142076\t0\t+\r\nchr22\t16275206\t16275277\tuc002zlh.1_cds_5_0_chr22_16275207_r\t0\t-\tchr22\t16275237\t16275238\trs2845214\t0\t+\r\nchr22\t16277747\t16277885\tuc002zlh.1_cds_6_0_chr22_16277748_r\t0\t-\tchr22\t16277818\t16277819\trs2073406\t0\t+\r\nchr22\t16277747\t16277885\tuc002zlh.1_cds_6_0_chr22_16277748_r\t0\t-\tchr22\t16277756\t16277757\trs79385954\t0\t+\r\n"
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Let's take a look at this dataset. The first six columns correspond to exons. The last six correspond to SNPs. You can see that exon with ID uc002zlh.1_cds_2_0_chr22_16266929_r contains four SNPs with IDs rs7290262, rs10154680, rs2818572, and rs7292200. "
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "2.1. Counting the number of SNPs per exon\n------------------------------------------\n\nAbove we've seen that exon uc002zlh.1_cds_2_0_chr22_16266929_r is repeated four times in the above dataset. Thus we can easily compute the number of SNPs per exon by simply counting the number of repetitions of name for each exon. This can be easily done with the \"Join, Subtract, and Group -> Group\" tool:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# python /home/azureuser/galaxy-dist/tools/stats/grouping.py        /home/azureuser/galaxy-dist/database/files/000/dataset_40.dat       /home/azureuser/galaxy-dist/database/files/000/dataset_39.dat       4       0        'length         4         no'\n\nparam1 = \"/home/azureuser/galaxy-dist/database/files/000/dataset_40.dat\"\nparam2 = \"/home/azureuser/galaxy-dist/database/files/000/dataset_39.dat\"\nparam3 = \"4\"\nparam4 = \"0\"\nparam5 = 'length         4         no'",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "if you look at the above image you will see that the result of grouping (dataset #4) contains two columns. This first contains the exon name while the second shows the number of times this name has been repeated in dataset #3. \n"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#!/usr/bin/env python\n# Guruprasad Ananda\n# Refactored 2011 to use numpy instead of rpy, Kanwei Li\n\"\"\"\nThis tool provides the SQL \"group by\" functionality.\n\"\"\"\nimport sys, commands, tempfile, random\ntry:\n    import numpy\nexcept:\n    from galaxy import eggs\n    eggs.require( \"numpy\" )\n    import numpy\n\nfrom itertools import groupby\n\ndef stop_err(msg):\n    sys.stderr.write(msg)\n    sys.exit()\n\ndef mode(data):\n    counts = {}\n    for x in data:\n        counts[x] = counts.get(x,0) + 1\n    maxcount = max(counts.values())\n    modelist = []\n    for x in counts:\n        if counts[x] == maxcount:\n            modelist.append( str(x) )\n    return ','.join(modelist)\n\ndef main():\n    #inputfile = sys.argv[2]\n    #ignorecase = int(sys.argv[4])\n    global param1, param2, param3, param4, param5\n    inputfile = param2\n    ignorecase = int(param4)\n    ops = []\n    cols = []\n    round_val = []\n    data_ary = []\n\n    #for var in sys.argv[5:]:\n    var = param5\n    op, col, do_round = var.split()\n    ops.append(op)\n    cols.append(col)\n    round_val.append(do_round)\n\n    \"\"\"\n    At this point, ops, cols and rounds will look something like this:\n    ops:  ['mean', 'min', 'c']\n    cols: ['1', '3', '4']\n    round_val: ['no', 'yes' 'no']\n    \"\"\"\n   \n    try:\n        #group_col = int( sys.argv[3] )-1\n        group_col = int( param3 )-1\n    except:\n        stop_err( \"Group column not specified.\" )\n\n    str_ops = ['c', 'length', 'unique', 'random', 'cuniq', 'Mode'] #ops that can handle string/non-numeric inputs\n\n    tmpfile = tempfile.NamedTemporaryFile()\n\n    try:\n        \"\"\"\n        The -k option for the Posix sort command is as follows:\n        -k, --key=POS1[,POS2]\n        start a key at POS1, end it at POS2 (origin 1)\n        In other words, column positions start at 1 rather than 0, so\n        we need to add 1 to group_col.\n        if POS2 is not specified, the newer versions of sort will consider the entire line for sorting. To prevent this, we set POS2=POS1.\n        \"\"\"\n        case = ''\n        if ignorecase == 1:\n            case = '-f'\n        command_line = \"sort -t '       ' %s -k%s,%s -o %s %s\" % (case, group_col+1, group_col+1, tmpfile.name, inputfile)\n    except Exception, exc:\n        stop_err( 'Initialization error -> %s' %str(exc) )\n\n    error_code, stdout = commands.getstatusoutput(command_line)\n\n    if error_code != 0:\n        stop_err( \"Sorting input dataset resulted in error: %s: %s\" %( error_code, stdout ))\n\n    #fout = open(sys.argv[1], \"w\")\n    fout = open(param1, \"w\")\n\n    def is_new_item(line):\n        try:\n            item = line.strip().split(\"\\t\")[group_col]\n        except IndexError:\n            stop_err( \"The following line didn't have %s columns: %s\" % (group_col+1, line) )\n\n        if ignorecase == 1:\n            return item.lower()\n        return item\n\n    for key, line_list in groupby(tmpfile, key=is_new_item):\n        op_vals = [ [] for op in ops ]\n        out_str = key\n        multiple_modes = False\n        mode_index = None\n\n        for line in line_list:\n            fields = line.strip().split(\"\\t\")\n            for i, col in enumerate(cols):\n                col = int(col)-1 # cXX from galaxy is 1-based\n                try:\n                    val = fields[col].strip()\n                    op_vals[i].append(val)\n                except IndexError:\n                    sys.stderr.write( 'Could not access the value for column %s on line: \"%s\". Make sure file is tab-delimited.\\n' % (col+1, line) )\n                    sys.exit( 1 )\n\n        # Generate string for each op for this group\n        for i, op in enumerate( ops ):\n            data = op_vals[i]\n            rval = \"\"\n            if op == \"mode\":\n                rval = mode( data )\n            elif op == \"length\":\n                rval = len( data )\n            elif op == \"random\":\n                rval = random.choice(data)\n            elif op in ['cat', 'cat_uniq']:\n                if op == 'cat_uniq':\n                    data = numpy.unique(data)\n                rval = ','.join(data)\n                try:\n                    val = fields[col].strip()\n                    op_vals[i].append(val)\n                except IndexError:\n                    sys.stderr.write( 'Could not access the value for column %s on line: \"%s\". Make sure file is tab-delimited.\\n' % (col+1, line) )\n                    sys.exit( 1 )\n\n        # Generate string for each op for this group\n        for i, op in enumerate( ops ):\n            data = op_vals[i]\n            rval = \"\"\n            if op == \"mode\":\n                rval = mode( data )\n            elif op == \"length\":\n                rval = len( data )\n            elif op == \"random\":\n                rval = random.choice(data)\n            elif op in ['cat', 'cat_uniq']:\n                if op == 'cat_uniq':\n                    data = numpy.unique(data)\n                rval = ','.join(data)\n            elif op == \"unique\":\n                rval = len( numpy.unique(data) )\n            else:\n                # some kind of numpy fn\n                try:\n                    data = map(float, data)\n                except ValueError:\n                    sys.stderr.write( \"Operation %s expected number values but got %s instead.\\n\" % (op, data) )\n                    sys.exit( 1 )\n                rval = getattr(numpy, op)( data )\n                if round_val[i] == 'yes':\n                    rval = round(rval)\n                else:\n                    rval = '%g' % rval\n\n            out_str += \"\\t%s\" % rval\n\n        fout.write(out_str + \"\\n\")\n\n    # Generate a useful info message.\n    msg = \"--Group by c%d: \" %(group_col+1)\n    for i, op in enumerate(ops):\n        if op == 'cat':\n            op = 'concat'\n        elif op == 'cat_uniq':\n            op = 'concat_distinct'\n        elif op == 'length':\n            op = 'count'\n        elif op == 'unique':\n            op = 'count_distinct'\n        elif op == 'random':\n            op = 'randomly_pick'\n\n        msg += op + \"[c\" + cols[i] + \"] \"\n\n    print msg\n    fout.close()\n    tmpfile.close()\n\nif __name__ == \"__main__\":\n    main()",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "SystemExit",
       "evalue": "",
       "output_type": "pyerr",
       "traceback": [
        "An exception has occurred, use %tb to see the full traceback.\n",
        "\u001b[0;31mSystemExit\u001b[0m\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": "Sorting input dataset resulted in error: 512: sort: multi-character tab `       'To exit: use 'exit', 'quit', or Ctrl-D."
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "!head $param1",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "uc002zlh.1_cds_1_0_chr22_16258186_r\t1\r\nuc002zlh.1_cds_2_0_chr22_16266929_r\t4\r\nuc002zlh.1_cds_4_0_chr22_16269873_r\t1\r\nuc002zlh.1_cds_5_0_chr22_16275207_r\t2\r\nuc002zlh.1_cds_6_0_chr22_16277748_r\t5\r\nuc002zlh.1_cds_7_0_chr22_16279195_r\t2\r\nuc002zlj.1_cds_0_0_chr22_16266931_r\t4\r\nuc002zlj.1_cds_2_0_chr22_16269873_r\t1\r\nuc002zlj.1_cds_3_0_chr22_16275207_r\t2\r\nuc002zlj.1_cds_4_0_chr22_16277748_r\t5\r\n"
      }
     ],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "2.2. Sorting exons by SNP count\n--------------------------------\n\nTo see which exon has the highest number of SNPs we can simply sort the dataset #4 on the second column in descending order. This is done with \"Filter and Sort -> Sort\":"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#python /home/azureuser/galaxy-dist/tools/filters/sorter.py         --input=/home/azureuser/galaxy-dist/database/files/000/dataset_40.dat         --output=/home/azureuser/galaxy-dist/database/files/000/dataset_41.dat                  --key=2,2nr\ninputfile = \"/home/azureuser/galaxy-dist/database/files/000/dataset_40.dat\"\noutputfile = \"/home/azureuser/galaxy-dist/database/files/000/dataset_41.dat\"\nkeyinput = \"2,2nr\"",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "\"\"\"\n    Sorts tabular data on one or more columns. All comments of the file are collected\n    and placed at the beginning of the sorted output file.\n\n    usage: sorter.py [options]\n    -i, --input: Tabular file to be sorted\n    -o, --output: Sorted output file\n    -k, --key: Key (see manual for bash/sort)\n\n    usage: sorter.py input output [key ...]\n\"\"\"\n# 03/05/2013 guerler\n\n# imports\nimport os, re, string, sys\nfrom optparse import OptionParser\n\n# error\ndef stop_err( msg ):\n    sys.stderr.write( \"%s\\n\" % msg )\n    sys.exit()\n\n# main\ndef main():\n    # define options\n    parser = OptionParser()\n    parser.add_option(\"-i\", \"--input\")\n    parser.add_option(\"-o\", \"--output\")\n    parser.add_option(\"-k\", \"--key\", action=\"append\")\n\n    # parse\n    #options, args = parser.parse_args()\n    \n    global inputfile, outputfile, keyinput\n\n    try:\n        # retrieve options\n        #input   = options.input\n        #output  = options.output\n        #key     = [\" -k\" + k for k in options.key]\n        input = inputfile\n        output = outputfile\n        key = [\" -k\" + k for k in [keyinput]]\n\n        # grep comments\n        grep_comments = \"(grep '^#' %s) > %s\" % (input, output)\n        #print grep_comments\n\n        # grep and sort columns\n        sort_columns  = \"(grep '^[^#]' %s | sort -f -t '\\t' %s) >> %s\" % (input, ' '.join(key), output)\n        #print sort_columns\n\n        # execute\n        os.system(grep_comments)\n        os.system(sort_columns)\n\n    except Exception, ex:\n        stop_err('Error running sorter.py\\n' + str(ex))\n\n    # exit\n    sys.exit(0)\n\nif __name__ == \"__main__\":\n    main()",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "SystemExit",
       "evalue": "0",
       "output_type": "pyerr",
       "traceback": [
        "An exception has occurred, use %tb to see the full traceback.\n",
        "\u001b[0;31mSystemExit\u001b[0m\u001b[0;31m:\u001b[0m 0\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": "To exit: use 'exit', 'quit', or Ctrl-D."
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "!head $outputfile",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "uc010gsw.2_cds_1_0_chr22_21480537_r\t67\r\nuc021wmb.1_cds_0_0_chr22_21480537_r\t67\r\nuc002zoc.3_cds_0_0_chr22_18834445_f\t58\r\nuc021wnd.1_cds_0_0_chr22_24647973_f\t50\r\nuc021wmc.1_cds_0_0_chr22_21637809_f\t47\r\nuc003bhh.3_cds_0_0_chr22_46652458_r\t46\r\nuc002zsd.4_cds_0_0_chr22_20456382_r\t41\r\nuc002zuq.4_cds_0_0_chr22_21738148_f\t41\r\nuc002zuy.4_cds_0_0_chr22_21900346_r\t41\r\nuc003atq.1_cds_12_0_chr22_38119192_f\t37\r\n"
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "2.4. Selecting top five\n-----------------------\n\nNow let's select top five exons with the highest number of SNPs. For this we will use \"Text Manipulation -> Select First\" tool:"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "inputfilename = \"/home/azureuser/galaxy-dist/database/files/000/dataset_41.dat\"\nlinenum = 5\noutputfilename = \"/home/azureuser/galaxy-dist/database/files/000/dataset_42.dat\"",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 37
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#! /usr/bin/perl -w\n\nuse strict;\nuse warnings;\n\n# a wrapper for head for use in galaxy\n# headWrapper.pl [filename] [# lines to show] [output]\n\ndie \"Check arguments\" unless @ARGV == 3;\ndie \"Line number must be an integer\\n\" unless $ARGV[1]=~ m/^\\d+$/;\n\nopen (OUT, \">$ARGV[2]\") or die \"Cannot create $ARGV[2]:$!\\n\";\nopen (HEAD, \"head -n $ARGV[1] $ARGV[0]|\") or die \"Cannot run head:$!\\n\";\nwhile (<HEAD>) {\n    print OUT;\n}\nclose OUT;\nclose HEAD;\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "SyntaxError",
       "evalue": "invalid syntax (<ipython-input-36-1b4ab6e90b60>, line 3)",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-36-1b4ab6e90b60>\"\u001b[0;36m, line \u001b[0;32m3\u001b[0m\n\u001b[0;31m    use strict;\u001b[0m\n\u001b[0m             ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
       ]
      }
     ],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "!head $outputfilename",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "uc010gsw.2_cds_1_0_chr22_21480537_r\t67\r\nuc021wmb.1_cds_0_0_chr22_21480537_r\t67\r\nuc002zoc.3_cds_0_0_chr22_18834445_f\t58\r\nuc021wnd.1_cds_0_0_chr22_24647973_f\t50\r\nuc021wmc.1_cds_0_0_chr22_21637809_f\t47\r\n"
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "2.5. Recovering exon info and displaying data in genome browsers\n----------------------------------------------------------------\n\nNow we know that in this dataset the five top exons contain between 41 and 67 SNPs. But what else can we learn about these? To know more we need to get back the positional information (coordinates) of these exons. This information was lost at the grouping step and now all we have is just two columns. To get coordinates back we will match the names of exons in dataset #6 (column 1) against names of the exons in the original dataset #1 (column 4). This can be done with \"Join, Subtract and Group -> Compare two Queries\" tool (note the settings of the tool in the middle pane)"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#python /home/azureuser/galaxy-dist/tools/filters/joinWrapper.py /home/azureuser/galaxy-dist/database/files/000/dataset_37.dat /home/azureuser/galaxy-dist/database/files/000/dataset_42.dat 4 1 N /home/azureuser/galaxy-dist/database/files/000/dataset_43.dat\ninfile1 = \"/home/azureuser/galaxy-dist/database/files/000/dataset_37.dat\"\ninfile2 = \"/home/azureuser/galaxy-dist/database/files/000/dataset_42.dat\"\nfield1 = 4\nfield2 = 1\nmode = \"N\"\noutfile = \"/home/azureuser/galaxy-dist/database/files/000/dataset_43.dat\"",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "#!/usr/bin/env python\n#Guruprasad Ananda\n\"\"\"\nThis tool provides the UNIX \"join\" functionality.\n\"\"\"\nimport sys, os, tempfile, subprocess\n\ndef stop_err(msg):\n    sys.stderr.write(msg)\n    sys.exit()\n\ndef main():\n    global infile1, infile2, field1, field2, mode, outfile\n    \n    #infile1 = sys.argv[1]\n    #infile2 = sys.argv[2]\n    #field1 = int(sys.argv[3])\n    #field2 = int(sys.argv[4])\n    #mode =sys.argv[5]\n    #outfile = sys.argv[6]\n\n    tmpfile1 = tempfile.NamedTemporaryFile()\n    tmpfile2 = tempfile.NamedTemporaryFile()\n\n    try:\n        #Sort the two files based on specified fields\n        os.system(\"sort -t '    ' -k %d,%d -o %s %s\" %(field1, field1, tmpfile1.name, infile1))\n        os.system(\"sort -t '    ' -k %d,%d -o %s %s\" %(field2, field2, tmpfile2.name, infile2))\n    except Exception, exc:\n        stop_err( 'Initialization error -> %s' %str(exc) )\n\n    option = \"\"\n    for line in file(tmpfile1.name):\n        line = line.strip()\n        if line:\n            elems = line.split('\\t')\n            for j in range(1,len(elems)+1):\n                if j == 1:\n                    option = \"1.1\"\n                else:\n                    option = option + \",1.\" + str(j)\n            break\n\n    #check if join has --version option. BSD join doens't have this option, while GNU join does.\n    #The return value in the latter case will be 0, and non-zero in the latter case.\n    ret = subprocess.call('join --version 2>/dev/null', shell=True)\n    # check if we are a version later than 7 of join. If so, we want to skip\n    # checking the order since join will raise an error with duplicated items in\n    # the two files being joined.\n    if ret == 0:\n        cl = subprocess.Popen([\"join\", \"--version\"], stdout=subprocess.PIPE)\n        (stdout, _) = cl.communicate()\n        version_line = stdout.split(\"\\n\")[0]\n        (version, _) = version_line.split()[-1].split(\".\")\n        if int(version) >= 7:\n            flags = \"--nocheck-order\"\n        else:\n            flags = \"\"\n    else:\n        flags = \"\"\n\n    if mode == \"V\":\n        cmdline = \"join %s -t ' ' -v 1 -o %s -1 %d -2 %d %s %s > %s\" %(flags, option, field1, field2, tmpfile1.name, tmpfile2.name, outfile)\n    else:\n        cmdline = \"join %s -t ' ' -o %s -1 %d -2 %d %s %s > %s\" %(flags, option, field1, field2, tmpfile1.name, tmpfile2.name, outfile)\n\n    try:\n        os.system(cmdline)\n    except Exception, exj:\n        stop_err('Error joining the two datasets -> %s' %str(exj))\n\nif __name__ == \"__main__\":\n    main()\n",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "join --nocheck-order -t ' ' -o  -1 4 -2 1 /tmp/tmpwCn7om /tmp/tmpwZTKoj > /home/azureuser/galaxy-dist/database/files/000/dataset_43.dat"
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "\n"
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "!head $outfile",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "chr22\t18834444\t18835833\tuc002zoc.3_cds_0_0_chr22_18834445_f\t0\t+\r\nchr22\t21480536\t21481925\tuc010gsw.2_cds_1_0_chr22_21480537_r\t0\t-\r\nchr22\t21480536\t21481925\tuc021wmb.1_cds_0_0_chr22_21480537_r\t0\t-\r\nchr22\t21637808\t21638558\tuc021wmc.1_cds_0_0_chr22_21637809_f\t0\t+\r\nchr22\t24647972\t24649256\tuc021wnd.1_cds_0_0_chr22_24647973_f\t0\t+\r\n"
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "The best way to learn about these exons is to look at their genomic surrounding. There is really no better way to do this than using genome browsers. Because this analysis was performed on \"standard\" human genome, you have two choices - UCSC Genome Browser and Ensembl:\n\n\n\nFor example, clicking on \"display at UCSC main\" will show this (to see your regions look at \"User Track\" on top of browser image):"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "",
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}