adar1_QC.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:602be4138db2bd02415d71c8193da50b9a2ee68f58ec274c9ead9386c483734e"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from gscripts.general import parsers\n",
      "from IPython.core.display import HTML\n",
      "import pandas as pd\n",
      "\n",
      "from gscripts.general import dataviz\n",
      "from gscripts.clipseq import clipseq_barcode_metrics"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "iclip_metrics = parsers.clipseq_metrics(\"/nas3/gpratt/projects/adar1/analysis/v1/\")\n",
      "iclip_metrics = iclip_metrics.ix[['index12', 'index6']]\n",
      "iclip_metrics = iclip_metrics.rename({'index12': 'Wild_type_ADR_1_CLIP', 'index6': 'ds1_2_mutant_ADR_1'})"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Adar1 FastQC Metrics\n",
      "\n",
      "http://sauron.ucsd.edu/4HH/index6.GCCAAT_fastqc/fastqc_report.html -- ds1_2 mutant\n",
      " \n",
      "http://sauron.ucsd.edu/4HH/index12.CTTGTA_fastqc/fastqc_report.html --Wild Type adar1 clip\n",
      "\n",
      "Main Takeaways:\n",
      "----------------\n",
      "The majority of Reads are Primer Dimers in both datasets, and this CLIP is unlikely to work\n",
      "\n",
      "For due dillagance I've mapped and called peaks on the Wild Type and ds1 mutants, mapping and peak calling QC is below"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Adar1 Mapping Metrics"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "filtered_df = iclip_metrics[[\"Input Reads\", \"Reads Passing Quality Filter\", \"Uniquely Mapped Reads\", \"Usable Reads\", 'Num Peaks']]                            \n",
      "HTML(filtered_df.to_html(formatters={\"Input Reads\" : parsers.commas, \"Reads After Triming\" : parsers.commas,\n",
      "                                     \"Reads Passing Quality Filter\" : parsers.commas, \"Uniquely Mapped Reads\" : parsers.commas,\n",
      "                                     \"Usable Reads\" : parsers.commas, \"Num Peaks\": parsers.commas,} ))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>Input Reads</th>\n",
        "      <th>Reads Passing Quality Filter</th>\n",
        "      <th>Uniquely Mapped Reads</th>\n",
        "      <th>Usable Reads</th>\n",
        "      <th>Num Peaks</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>Wild_type_ADR_1_CLIP</th>\n",
        "      <td>82,621,817</td>\n",
        "      <td> 81,779,038</td>\n",
        "      <td>   331,276</td>\n",
        "      <td>5,273</td>\n",
        "      <td>16</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>ds1_2_mutant_ADR_1</th>\n",
        "      <td>57,144,044</td>\n",
        "      <td> 47,768,030</td>\n",
        "      <td> 2,379,127</td>\n",
        "      <td>5,173</td>\n",
        "      <td> 9</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<IPython.core.display.HTML at 0x2aaaac2b7910>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Main Takeaways:\n",
      "----------------\n",
      "\n",
      "The Majority of Adar1 CLIP-seq reads do not map, due to the large number of primer dimers.  Those that do map are mostly PCR duplicates and removed during barcode collapsing. Leaving peak calling to be performed on 5000 reads.  When peak calling is performed very few peaks are identified (16 and 9).  I would not recommend this data for followup analysis"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#MSI CLIP Results"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!scp -q tscc-login.sdsc.edu:/home/gpratt/projects/adar1/analysis/v1/*{metrics,peaks.bed,final.out} /nas3/gpratt/projects/adar1/analysis/v1/"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "discarding /projects/ps-yeolab/software/anaconda-2.0.1_2014-08-12/bin from PATH\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "prepending /projects/ps-yeolab/software/anaconda-2.0.1_2014-08-12/envs/gpratt/bin to PATH\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "discarding /projects/ps-yeolab/software/anaconda-2.0.1_2014-08-12/bin from PATH\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "prepending /projects/ps-yeolab/software/anaconda-2.0.1_2014-08-12/envs/gpratt/bin to PATH\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "discarding /projects/ps-yeolab/software/anaconda-2.0.1_2014-08-12/bin from PATH\r\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "prepending /projects/ps-yeolab/software/anaconda-2.0.1_2014-08-12/envs/gpratt/bin to PATH\r\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "iclip_metrics = parsers.clipseq_metrics(\"/nas3/gpratt/projects/adar1/analysis/v1/\")\n",
      "iclip_metrics = iclip_metrics.drop([\"Wild_type_ADR_1_CLIP\", \"ds1_2_mutant_ADR_1\"])\n",
      "iclip_metrics = iclip_metrics.rename({'index12': 'Wild_type_ADR_1_CLIP', 'index6': 'ds1_2_mutant_ADR_1'})\n",
      "filtered_df = iclip_metrics[[\"Input Reads\", \"Reads Passing Quality Filter\", \"Uniquely Mapped Reads\", \"Usable Reads\", 'Num Peaks', u'Percent Usable / Input',\n",
      "]]                            \n",
      "HTML(filtered_df.to_html(formatters={\"Input Reads\" : parsers.commas, \"Reads After Triming\" : parsers.commas,\n",
      "                                     \"Reads Passing Quality Filter\" : parsers.commas, \"Uniquely Mapped Reads\" : parsers.commas,\n",
      "                                     \"Usable Reads\" : parsers.commas, \"Num Peaks\": parsers.commas,} ))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>Input Reads</th>\n",
        "      <th>Reads Passing Quality Filter</th>\n",
        "      <th>Uniquely Mapped Reads</th>\n",
        "      <th>Usable Reads</th>\n",
        "      <th>Num Peaks</th>\n",
        "      <th>Percent Usable / Input</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>MSI2_ACAGTG_ACAGTG_L008_R1</th>\n",
        "      <td>60,059,965</td>\n",
        "      <td> 52,081,382</td>\n",
        "      <td> 34,034,711</td>\n",
        "      <td>406,600</td>\n",
        "      <td>23,781</td>\n",
        "      <td>  0.01194663</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>MSI2_CAGATC_CAGATC_L008_R1</th>\n",
        "      <td>59,001,980</td>\n",
        "      <td> 53,837,810</td>\n",
        "      <td>  5,024,516</td>\n",
        "      <td> 31,995</td>\n",
        "      <td> 1,010</td>\n",
        "      <td> 0.006367778</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>MSI2_GCCAAT_GCCAAT_L008_R1</th>\n",
        "      <td>   152,525</td>\n",
        "      <td>    123,931</td>\n",
        "      <td>     26,252</td>\n",
        "      <td> 14,738</td>\n",
        "      <td>   279</td>\n",
        "      <td>   0.5614048</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>MSI2_TAGCTT_TAGCTT_L008_R1</th>\n",
        "      <td>    77,513</td>\n",
        "      <td>     60,517</td>\n",
        "      <td>     25,779</td>\n",
        "      <td> 18,625</td>\n",
        "      <td>   396</td>\n",
        "      <td>   0.7224873</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>MSI2_TGACCA_TGACCA_L008_R1</th>\n",
        "      <td>   101,367</td>\n",
        "      <td>     78,734</td>\n",
        "      <td>     37,698</td>\n",
        "      <td> 28,202</td>\n",
        "      <td>   693</td>\n",
        "      <td>   0.7481033</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>MSI2_TTAGGC_TTAGGC_L008_R1</th>\n",
        "      <td>86,861,473</td>\n",
        "      <td> 67,159,228</td>\n",
        "      <td> 37,499,177</td>\n",
        "      <td>288,619</td>\n",
        "      <td>21,809</td>\n",
        "      <td> 0.007696676</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>Wild_type_ADR_1_CLIP</th>\n",
        "      <td>82,621,817</td>\n",
        "      <td> 81,779,038</td>\n",
        "      <td>    331,276</td>\n",
        "      <td>  5,273</td>\n",
        "      <td>    16</td>\n",
        "      <td>  0.01591724</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>ds1_2_mutant_ADR_1</th>\n",
        "      <td>57,144,044</td>\n",
        "      <td> 47,768,030</td>\n",
        "      <td>  2,379,127</td>\n",
        "      <td>  5,173</td>\n",
        "      <td>     9</td>\n",
        "      <td> 0.002174327</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "<IPython.core.display.HTML at 0x2aaaac3410d0>"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Main Takeaways:\n",
      "----------------\n",
      "\n",
      "1. Samples GCCAAT, TAGCTT and TGACCA failed due to a lack of sequencing depth\n",
      "2. Sample CAGATC failed due to a low percent usable usable reads and a low number of usable reads overall\n",
      "3. The other two samples are borderline.  If it were my data I wouldn't trust them due to a low number of usable reads, but they might be able to be salvaged if there is enrichment in specific known regions or results are consistant with previously published data\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}