Region_Classification_Example.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:efd768856ce394ae0d7e074d792c56bd296887e1763bc21003f5dd941e6ce9b7"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from collections import OrderedDict, Counter\n",
      "from clipper.src import CLIP_analysis\n",
      "import pybedtools\n",
      "from gscripts.general import region_helpers"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "regions = OrderedDict()\n",
      "regions[\"cds\"] = \"CDS\"\n",
      "regions[\"three_prime_utrs\"] = \"3' UTR\"\n",
      "regions[\"five_prime_utrs\"] = \"5' UTR\"\n",
      "regions[\"proxintron500\"] = \"Proximal\\nIntron\"\n",
      "regions[\"distintron500\"] = \"Distal\\nIntron\"\n",
      "\n",
      "bedtool = pybedtools.BedTool(\"/Users/gpratt/bioinformatics/yeolab/projects/public_clip/analysis/v6/293T_hnRNPA1_WT_rep1.peaks.bed\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Classify each peak to a different region"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cluster_regions = CLIP_analysis.assign_to_regions(bedtool, \"overlapping\", \n",
      "                                    regions, \"assigned\", nrand=0, species=\"mm9\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "There are a total 88 clusters I'll examine\n",
        "For region: cds found 6 that overlap and 82 that don't"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "ignoring three_prime_utrs "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "ignoring five_prime_utrs "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "For region: proxintron500 found 14 that overlap and 68 that don't"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "For region: distintron500 found 34 that overlap and 34 that don't"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "After assigning 34 un-categorized regions\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print cluster_regions"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "{'proxintron500': {'real': <BedTool(assigned/overlapping.proxintron500.real.BED)>, 'rand': {}}, 'all': {'real': <BedTool(assigned/overlapping.all.real.BED)>, 'rand': {}}, 'uncatagorized': {'real': <BedTool(assigned/overlapping.uncatagorized.real.BED)>}, 'distintron500': {'real': <BedTool(assigned/overlapping.distintron500.real.BED)>, 'rand': {}}, 'cds': {'real': <BedTool(assigned/overlapping.cds.real.BED)>, 'rand': {}}}\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Count Number of Reads in each Region"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for region in cluster_regions.keys():\n",
      "    print region, len(cluster_regions[region]['real'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "proxintron500 14\n",
        "all 54\n",
        "uncatagorized 34\n",
        "distintron500 34\n",
        "cds 6\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Count Transcript Types"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gene_id_to_type = region_helpers.gene_id_to_type(\"/Users/gpratt/bioinformatics/genomes/hg19/gencode.v17.annotation.gtf.db\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bedtool.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "chr1\t11968283\t11968303\tENSG00000199347.1_1_10\t5.00579720837e-05\t+\t11968287\t11968291\n",
        " chr1\t36850227\t36850254\tENSG00000196182.6_1_8\t0.0226923149445\t-\t36850232\t36850236\n",
        " chr1\t71534206\t71534232\tENSG00000132485.8_1_16\t8.04450376755e-11\t-\t71534213\t71534217\n",
        " chr1\t173834525\t173834542\tENSG00000200729.1_1_10\t5.00579720837e-05\t-\t173834532\t173834536\n",
        " chr1\t173834527\t173834544\tENSG00000234741.2_1_13\t1.0111184346e-07\t-\t173834530\t173834534\n",
        " chr1\t173836041\t173836063\tENSG00000234741.2_3_9\t0.00206019486484\t-\t173836048\t173836052\n",
        " chr1\t235273608\t235273638\tENSG00000173726.6_1_8\t0.0259646963199\t-\t235273621\t235273625\n",
        " chr1\t245013700\t245013724\tENSG00000188206.5_4_8\t0.0197801382447\t-\t245013707\t245013711\n",
        " chr1\t245015966\t245015997\tENSG00000188206.5_7_10\t0.000357718415964\t-\t245015978\t245015982\n",
        " chr1\t245015966\t245015997\tENSG00000153187.12_3_10\t0.000357718415964\t-\t245015978\t245015982\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The Name before the underscore is the gene id, we get it by futzing with the string"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "[interval.name.split(\"_\")[0] for interval in bedtool][:10]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "['ENSG00000199347.1',\n",
        " 'ENSG00000196182.6',\n",
        " 'ENSG00000132485.8',\n",
        " 'ENSG00000200729.1',\n",
        " 'ENSG00000234741.2',\n",
        " 'ENSG00000234741.2',\n",
        " 'ENSG00000173726.6',\n",
        " 'ENSG00000188206.5',\n",
        " 'ENSG00000188206.5',\n",
        " 'ENSG00000153187.12']"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now Lets get the types, I have a dict / helper function that maps gene_ids to gene type"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "[gene_id_to_type[interval.name.split(\"_\")[0]] for interval in bedtool][:10]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "[u'snRNA',\n",
        " u'protein_coding',\n",
        " u'protein_coding',\n",
        " u'snoRNA',\n",
        " u'processed_transcript',\n",
        " u'processed_transcript',\n",
        " u'protein_coding',\n",
        " u'processed_transcript',\n",
        " u'processed_transcript',\n",
        " u'protein_coding']"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And now we can count the total number of these things, there are more obvious ways to do this, but python provides a counter class that does this for us"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Counter([gene_id_to_type[interval.name.split(\"_\")[0]] for interval in bedtool])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "Counter({u'lincRNA': 71, u'protein_coding': 15, u'Mt_rRNA': 14, u'processed_transcript': 5, u'rRNA': 1, u'snRNA': 1, u'sense_overlapping': 1, u'miRNA': 1, u'misc_RNA': 1, u'snoRNA': 1})"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}