{
 "metadata": {
  "name": "",
  "signature": "sha256:852f4ae9fec9eb5882c43884225447b904940bba1079514a1ee20027bbc6c732"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Identification of Differentially Methylated loci in Oyster larvae (mix v control)"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "methratio output"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#mix\n",
      "#details @ http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/fish546/BSseq_workflow.ipynb\n",
      "!head /Volumes/web/cnidarian/YE_mix_22sm_methratio_out.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "chr\tpos\tstrand\tcontext\tratio\teff_CT_count\tC_count\tCT_count\trev_G_count\trev_GA_count\tCI_lower\tCI_upper\r\n",
        "C12706\t112\t+\tTTCTA\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t135\t+\tGACTC\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t137\t+\tCTCTT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t148\t+\tATCAA\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t151\t+\tAACTT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t159\t+\tAGCCT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t160\t+\tGCCTA\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t163\t+\tTACAA\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12706\t170\t-\tAAGGT\t0.000\t1.00\t0\t1\t2\t2\t0.000\t0.793\r\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#control\n",
      "#details @ http://nbviewer.ipython.org/github/sr320/ipython_nb/blob/master/fish546/BSseq_YE_control.ipynb\n",
      "!head /Volumes/web/cnidarian/YE_control_22sm_methratio_out.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "chr\tpos\tstrand\tcontext\tratio\teff_CT_count\tC_count\tCT_count\trev_G_count\trev_GA_count\tCI_lower\tCI_upper\r\n",
        "C12798\t41\t+\tAACAA\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t45\t+\tAACTT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t52\t+\tATCCC\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t53\t+\tTCCCC\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t54\t+\tCCCCT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t55\t+\tCCCTT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t58\t+\tTTCAC\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t60\t+\tCACCT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12798\t61\t+\tACCTA\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "CG only filter"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!grep \"[A-Z][A-Z]CG[A-Z]\" </Volumes/web/cnidarian/YE_mix_22sm_methratio_out.txt> /Volumes/web/cnidarian/YE_mix_22sm_methratio_CG.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!head /Volumes/web/cnidarian/YE_mix_22sm_methratio_CG.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C12722\t104\t+\tTACGT\t0.500\t2.00\t1\t2\t2\t2\t0.095\t0.905\r\n",
        "C12728\t40\t+\tGCCGC\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12728\t62\t+\tGCCGT\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12734\t61\t+\tACCGT\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12734\t136\t+\tACCGC\t0.000\t2.00\t0\t2\t2\t2\t0.000\t0.658\r\n",
        "C12734\t178\t+\tGACGG\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12734\t183\t+\tGGCGG\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12748\t71\t+\tAACGA\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12758\t82\t+\tTCCGC\t1.000\t1.00\t1\t1\t1\t1\t0.207\t1.000\r\n",
        "C12758\t125\t+\tTTCGC\t1.000\t1.00\t1\t1\t1\t1\t0.207\t1.000\r\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!grep \"[A-Z][A-Z]CG[A-Z]\" </Volumes/web/cnidarian/YE_control_22sm_methratio_out.txt> /Volumes/web/cnidarian/YE_control_22sm_methratio_CG.txt\n",
      "!head /Volumes/web/cnidarian/YE_control_22sm_methratio_CG.txt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C12802\t82\t+\tACCGG\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12802\t90\t+\tTGCGA\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12824\t150\t+\tGGCGT\t0.000\t1.00\t0\t1\t0\t0\t0.000\t0.793\r\n",
        "C12824\t160\t+\tATCGC\t0.000\t2.00\t0\t2\t2\t2\t0.000\t0.658\r\n",
        "C12824\t183\t+\tTTCGA\t0.000\t2.00\t0\t2\t2\t2\t0.000\t0.658\r\n",
        "C12832\t9\t+\tTGCGT\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12832\t25\t+\tAGCGC\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12832\t129\t+\tGTCGA\t0.000\t2.00\t0\t2\t2\t2\t0.000\t0.658\r\n",
        "C12842\t123\t+\tCACGA\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n",
        "C12860\t54\t+\tAACGT\t0.000\t1.00\t0\t1\t1\t1\t0.000\t0.793\r\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "convert methratio to methylkit format"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!awk '{print $1,$2,$2+1,$3,$8,($7/$8),(1-($7/$8))}' </Volumes/web/cnidarian/YE_mix_22sm_methratio_CG.txt> /Volumes/web/cnidarian/YE_mix_22sm_methkit_outCG.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!head /Volumes/web/cnidarian/YE_mix_22sm_methkit_outCG.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C12722 104 105 + 2 0.5 0.5\r\n",
        "C12728 40 41 + 1 0 1\r\n",
        "C12728 62 63 + 1 0 1\r\n",
        "C12734 61 62 + 1 0 1\r\n",
        "C12734 136 137 + 2 0 1\r\n",
        "C12734 178 179 + 1 0 1\r\n",
        "C12734 183 184 + 1 0 1\r\n",
        "C12748 71 72 + 1 0 1\r\n",
        "C12758 82 83 + 1 1 0\r\n",
        "C12758 125 126 + 1 1 0\r\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#convert space to tab\n",
      "!tr ' ' \"\\t\" </Volumes/web/cnidarian/YE_mix_22sm_methkit_outCG.txt> /Volumes/web/cnidarian/YE_mix_22sm_methkit.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!head /Volumes/web/cnidarian/YE_mix_22sm_methkit.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C12722\t104\t105\t+\t2\t0.5\t0.5\r\n",
        "C12728\t40\t41\t+\t1\t0\t1\r\n",
        "C12728\t62\t63\t+\t1\t0\t1\r\n",
        "C12734\t61\t62\t+\t1\t0\t1\r\n",
        "C12734\t136\t137\t+\t2\t0\t1\r\n",
        "C12734\t178\t179\t+\t1\t0\t1\r\n",
        "C12734\t183\t184\t+\t1\t0\t1\r\n",
        "C12748\t71\t72\t+\t1\t0\t1\r\n",
        "C12758\t82\t83\t+\t1\t1\t0\r\n",
        "C12758\t125\t126\t+\t1\t1\t0\r\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#processing control sample\n",
      "!awk '{print $1,$2,$2+1,$3,$8,($7/$8),(1-($7/$8))}' </Volumes/web/cnidarian/YE_control_22sm_methratio_CG.txt> /Volumes/web/cnidarian/YE_control_22sm_methkit_outCG.txt\n",
      "!tr ' ' \"\\t\" </Volumes/web/cnidarian/YE_control_22sm_methkit_outCG.txt> /Volumes/web/cnidarian/YE_control_22sm_methkit.txt \n",
      "!head /Volumes/web/cnidarian/YE_control_22sm_methkit.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C12802\t82\t83\t+\t1\t0\t1\r\n",
        "C12802\t90\t91\t+\t1\t0\t1\r\n",
        "C12824\t150\t151\t+\t1\t0\t1\r\n",
        "C12824\t160\t161\t+\t2\t0\t1\r\n",
        "C12824\t183\t184\t+\t2\t0\t1\r\n",
        "C12832\t9\t10\t+\t1\t0\t1\r\n",
        "C12832\t25\t26\t+\t1\t0\t1\r\n",
        "C12832\t129\t130\t+\t2\t0\t1\r\n",
        "C12842\t123\t124\t+\t1\t0\t1\r\n",
        "C12860\t54\t55\t+\t1\t0\t1\r\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#filter mk files for 3x coverage\n",
      "!awk '{if ($5 >= 3) print $1,$2,$3,$4,$5,$6,$7}' </Volumes/web/cnidarian/YE_mix_22sm_methkit.txt> /Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt \n",
      "!head /Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C12764 119 120 + 3 0 1\r\n",
        "C12764 129 130 + 4 0 1\r\n",
        "C12764 132 133 + 4 0 1\r\n",
        "C12778 114 115 + 5 0 1\r\n",
        "C12802 119 120 + 4 0 1\r\n",
        "C12860 105 106 + 3 0 1\r\n",
        "C12860 113 114 + 4 0 1\r\n",
        "C12860 116 117 + 4 0 1\r\n",
        "C12860 121 122 + 3 0 1\r\n",
        "C12860 159 160 + 3 0 1\r\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#filter mk files for 3x coverage\n",
      "!awk '{if ($5 >= 3) print $1,$2,$3,$4,$5,$6,$7}' </Volumes/web/cnidarian/YE_control_22sm_methkit.txt> /Volumes/web/cnidarian/YE_control_22sm_methkit3x.txt \n",
      "!head /Volumes/web/cnidarian/YE_control_22sm_methkit3x.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C13288 101 102 + 4 0 1\r\n",
        "C13288 108 109 + 4 0 1\r\n",
        "C13448 37 38 + 3 0 1\r\n",
        "C13800 98 99 + 3 0 1\r\n",
        "C14308 114 115 + 3 0 1\r\n",
        "C14308 121 122 + 3 0 1\r\n",
        "C14966 94 95 + 3 0 1\r\n",
        "C15066 82 83 + 3 0.333333 0.666667\r\n",
        "C15066 89 90 + 3 0 1\r\n",
        "C15066 261 262 + 4 0 1\r\n"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Methylkit (R package)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext rpy2.ipython"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%R library(methylKit)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "<StrVector - Python:0x10a62c560 / R:0x10f303400>\n",
        "[str, str, str, ..., str, str, str]"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%R library(data.table)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "data.table 1.9.2  For help type: help(\"data.table\")\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%R library(GenomicRanges)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "Loading required package: BiocGenerics\n",
        "Loading required package: parallel\n",
        "\n",
        "Attaching package: \u2018BiocGenerics\u2019\n",
        "\n",
        "The following objects are masked from \u2018package:parallel\u2019:\n",
        "\n",
        "    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,\n",
        "    clusterExport, clusterMap, parApply, parCapply, parLapply,\n",
        "    parLapplyLB, parRapply, parSapply, parSapplyLB\n",
        "\n",
        "The following object is masked from \u2018package:stats\u2019:\n",
        "\n",
        "    xtabs\n",
        "\n",
        "The following objects are masked from \u2018package:base\u2019:\n",
        "\n",
        "    anyDuplicated, append, as.data.frame, as.vector, cbind, colnames,\n",
        "    duplicated, eval, evalq, Filter, Find, get, intersect, is.unsorted,\n",
        "    lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int,\n",
        "    pmin, pmin.int, Position, rank, rbind, Reduce, rep.int, rownames,\n",
        "    sapply, setdiff, sort, table, tapply, union, unique, unlist\n",
        "\n",
        "Loading required package: IRanges\n",
        "Loading required package: XVector\n",
        "\n",
        "Attaching package: \u2018GenomicRanges\u2019\n",
        "\n",
        "The following object is masked from \u2018package:data.table\u2019:\n",
        "\n",
        "    last\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!head /Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "C12764 119 120 + 3 0 1\r\n",
        "C12764 129 130 + 4 0 1\r\n",
        "C12764 132 133 + 4 0 1\r\n",
        "C12778 114 115 + 5 0 1\r\n",
        "C12802 119 120 + 4 0 1\r\n",
        "C12860 105 106 + 3 0 1\r\n",
        "C12860 113 114 + 4 0 1\r\n",
        "C12860 116 117 + 4 0 1\r\n",
        "C12860 121 122 + 3 0 1\r\n",
        "C12860 159 160 + 3 0 1\r\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R file.list <- list \n",
      "('/Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt',\n",
      " '/Volumes/web/cnidarian/YE_control_22sm_methkit3x.txt')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%R myobj=read(file.list,sample.id=list(\"mix\",\"control\"),assembly=\"v9\",treatment=c(0,1))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "<ListVector - Python:0x10a62c2d8 / R:0x10f94c588>\n",
        "[ListVector, ListVector]\n",
        "<ListVector - Python:0x10a62c2d8 / R:0x10f94c588>\n",
        "[ListVector, ListVector]\n",
        "<ListVector - Python:0x10a62c2d8 / R:0x10f94c588>\n",
        "[ListVector, ListVector]"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "meth<-unite(myobj)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "Error in dev.off() : \n",
        "  QuartzBitmap_Output - unable to open file '/var/folders/28/07h2q2ms7nd609f1ft3mj1240000gn/T/tmpHGwHzf/Rplots001.png'\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "PROBLEM WITH NAS"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "meth<-unite(myobj)\n",
      "head(meth)\n",
      "nrow(meth)\n",
      "getCorrelation(meth,plot=T)\n",
      "hc<- clusterSamples(meth, dist=\"correlation\", method=\"ward\", plot=T)\n",
      "PCA<-PCASamples(meth)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "        mix control\n",
        "mix       1      NA\n",
        "control  NA       1\n",
        "KernSmooth 2.23 loaded\n",
        "Copyright M. P. Wand 1997-2009\n",
        "Error in quantile.default(sds, sd.threshold) : \n",
        "  missing values and NaN's not allowed if 'na.rm' is FALSE\n"
       ]
      },
      {
       "ename": "RRuntimeError",
       "evalue": "Error in dev.off() : \n  QuartzBitmap_Output - unable to open file '/var/folders/28/07h2q2ms7nd609f1ft3mj1240000gn/T/tmpHGwHzf/Rplots001.png'\n",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[0;31mRRuntimeError\u001b[0m                             Traceback (most recent call last)",
        "\u001b[0;32m<ipython-input-10-b640763ed36e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun_cell_magic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mu'R'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34mu''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34mu'meth<-unite(myobj)\\nhead(meth)\\nnrow(meth)\\ngetCorrelation(meth,plot=T)\\nhc<- clusterSamples(meth, dist=\"correlation\", method=\"ward\", plot=T)\\nPCA<-PCASamples(meth)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
        "\u001b[0;32m/Users/sr320/anaconda/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc\u001b[0m in \u001b[0;36mrun_cell_magic\u001b[0;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[1;32m   2160\u001b[0m             \u001b[0mmagic_arg_s\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvar_expand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstack_depth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2161\u001b[0m             \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuiltin_trap\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2162\u001b[0;31m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmagic_arg_s\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcell\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2163\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2164\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/sr320/anaconda/lib/python2.7/site-packages/rpy2/ipython/rmagic.pyc\u001b[0m in \u001b[0;36mR\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n",
        "\u001b[0;32m/Users/sr320/anaconda/lib/python2.7/site-packages/IPython/core/magic.pyc\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(f, *a, **k)\u001b[0m\n\u001b[1;32m    191\u001b[0m     \u001b[0;31m# but it's overkill for just that one bit of state.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    192\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mmagic_deco\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 193\u001b[0;31m         \u001b[0mcall\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    194\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    195\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/sr320/anaconda/lib/python2.7/site-packages/rpy2/ipython/rmagic.pyc\u001b[0m in \u001b[0;36mR\u001b[0;34m(self, line, cell, local_ns)\u001b[0m\n\u001b[1;32m    634\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    635\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdevice\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m'png'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'svg'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 636\u001b[0;31m                 \u001b[0mro\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'dev.off()'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    637\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    638\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mtext_output\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/sr320/anaconda/lib/python2.7/site-packages/rpy2/robjects/__init__.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, string)\u001b[0m\n\u001b[1;32m    244\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    245\u001b[0m         \u001b[0mp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrinterface\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstring\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 246\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    247\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    248\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/sr320/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    164\u001b[0m                 \u001b[0mv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    165\u001b[0m                 \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mr_k\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 166\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mSignatureTranslatedFunction\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    167\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    168\u001b[0m \u001b[0mpattern_link\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mre\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr'\\\\link\\{(.+?)\\}'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;32m/Users/sr320/anaconda/lib/python2.7/site-packages/rpy2/robjects/functions.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     97\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     98\u001b[0m             \u001b[0mnew_kwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconversion\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpy2ri\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mFunction\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnew_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mnew_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconversion\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mri2ro\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mres\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    101\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
        "\u001b[0;31mRRuntimeError\u001b[0m: Error in dev.off() : \n  QuartzBitmap_Output - unable to open file '/var/folders/28/07h2q2ms7nd609f1ft3mj1240000gn/T/tmpHGwHzf/Rplots001.png'\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "</head>\n",
      "\n",
      "<body>\n",
      "<h1>DMR id with Methylkit</h1>\n",
      "\n",
      "<p>Install should only have to occur once.</p>\n",
      "\n",
      "<pre><code>source(&quot;http://methylkit.googlecode.com/files/install.methylKit.R&quot;)\n",
      "\n",
      "# installing the package with the dependencies\n",
      "install.methylKit(ver=&quot;0.9.2&quot;,dependencies=TRUE)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">library(methylKit)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">library(data.table)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">library(GenomicRanges)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code>## Loading required package: BiocGenerics\n",
      "## Loading required package: parallel\n",
      "## \n",
      "## Attaching package: &#39;BiocGenerics&#39;\n",
      "## \n",
      "## The following objects are masked from &#39;package:parallel&#39;:\n",
      "## \n",
      "##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,\n",
      "##     clusterExport, clusterMap, parApply, parCapply, parLapply,\n",
      "##     parLapplyLB, parRapply, parSapply, parSapplyLB\n",
      "## \n",
      "## The following object is masked from &#39;package:stats&#39;:\n",
      "## \n",
      "##     xtabs\n",
      "## \n",
      "## The following objects are masked from &#39;package:base&#39;:\n",
      "## \n",
      "##     anyDuplicated, append, as.data.frame, as.vector, cbind,\n",
      "##     colnames, duplicated, eval, evalq, Filter, Find, get,\n",
      "##     intersect, is.unsorted, lapply, Map, mapply, match, mget,\n",
      "##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,\n",
      "##     rbind, Reduce, rep.int, rownames, sapply, setdiff, sort,\n",
      "##     table, tapply, union, unique, unlist\n",
      "## \n",
      "## Loading required package: IRanges\n",
      "## Loading required package: XVector\n",
      "## \n",
      "## Attaching package: &#39;GenomicRanges&#39;\n",
      "## \n",
      "## The following object is masked from &#39;package:data.table&#39;:\n",
      "## \n",
      "##     last\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">file.list &lt;- list(&#39;/Volumes/web/cnidarian/YE_control_22sm_methkit3x.txt&#39;,&#39;/Volumes/web/cnidarian/YE_mix_22sm_methkit3x.txt&#39;)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">myobj=read(file.list,sample.id=list(&quot;mix&quot;,&quot;control&quot;),assembly=&quot;ensembl&quot;,treatment=c(1,0))\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">meth&lt;-unite(myobj)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">nrow(meth)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code>## [1] 309616\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">head(meth)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code>## methylBase object with 6 rows\n",
      "## --------------\n",
      "##   chr start end strand coverage1 numCs1 numTs1 coverage2 numCs2 numTs2\n",
      "## 1   5     6   6      +        11      0      0       110      0      1\n",
      "## 2   8     9   9      +         5      0      0         4      0      0\n",
      "## 3   8     9   9      +         5      0      0         8      0      0\n",
      "## 4   9    10  10      +         6      0      0        35      0      0\n",
      "## 5   9    10  10      +         6      0      0         3      0      0\n",
      "## 6  13    14  14      +         3      0      0         3      0      0\n",
      "## --------------\n",
      "## sample.ids: mix control \n",
      "## destranded FALSE \n",
      "## assembly: ensembl \n",
      "## context: CpG \n",
      "## treament: 1 0 \n",
      "## resolution: base\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">myDiff=calculateDiffMeth(meth)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">head(myDiff)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code>## methylDiff object with 6 rows\n",
      "## --------------\n",
      "##   chr start end strand pvalue qvalue meth.diff\n",
      "## 1   5     6   6      +      1      0         0\n",
      "## 2   8     9   9      +      1      0         0\n",
      "## 3   8     9   9      +      1      0         0\n",
      "## 4   9    10  10      +      1      0         0\n",
      "## 5   9    10  10      +      1      0         0\n",
      "## 6  13    14  14      +      1      0         0\n",
      "## --------------\n",
      "## sample.ids: mix control \n",
      "## destranded FALSE \n",
      "## assembly: ensembl \n",
      "## context: CpG \n",
      "## treament: 1 0 \n",
      "## resolution: base\n",
      "</code></pre>\n",
      "\n",
      "<pre><code class=\"r\">write.csv(myDiff, file=&quot;/Volumes/web/cnidarian/YE_mix_diff3.csv&quot;)\n",
      "</code></pre>\n",
      "\n",
      "<pre><code>## Warning: Setting class(x) to NULL; result will no longer be an S4 object\n",
      "</code></pre>\n",
      "\n",
      "<p><img src=\"http://eagle.fish.washington.edu/cnidarian/skitch/Screenshot_6_23_14_8_08_AM_19587A90.png\" alt=\"excel\"/></p>\n",
      "\n",
      "</body>\n",
      "\n",
      "</html>\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}