{
 "metadata": {
  "name": "",
  "signature": "sha256:56f4ea1a588a1856fd6ea5680da2c3e4ee336896e88b2b4a4bfbf806088a7e9f"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Reproducing results from Froelich et al"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext rpy2.ipython"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "from how to \n",
      "```\n",
      "Lab experiment zip file\n",
      "Contained in the laboratory zip file are CSV files of the normalized HIF-1\u03b1 data and associated R script. The Threshold and Time-course analyses and plots (Fig. 5 & 6) are saved in separate files.\n",
      "First, open the desired R script. Set the working directory so R can call the CSV files. Then simply run the code to view the results of the two-way ANOVA and recreation of the plots in the manuscript.\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pwd"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "u'/Users/sr320/Dropbox/Steven/ipython_nb'"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!cat /Users/sr320/Desktop/1190951/Lab\\ experiments/ANOVA_Threshold.r"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "#Set working directory\r",
        "\r\n",
        "#workdir <-\"C:/Users/\"\r",
        "\r\n",
        "setwd(workdir)\r",
        "\r\n",
        "\r",
        "\r\n",
        "##****************THRESHOLD EXPERIMENT DATA************************************\r",
        "\r\n",
        "threshold.data <-read.csv('Threshold_data.csv', header = TRUE)\r",
        "\r\n",
        "head(threshold.data)\r",
        "\r\n",
        "\r",
        "\r\n",
        "##****two-way ANOVA**************\r",
        "\r\n",
        "th.anova<-lm(NormalizedExp~Treatment + Size, data=threshold.data)\r",
        "\r\n",
        "summary(th.anova)\r",
        "\r\n",
        "anova(th.anova)\r",
        "\r\n",
        "\r",
        "\r\n",
        "##*********Plot Fig. 5 mean normalized and SE threshold data*****************\r",
        "\r\n",
        "library(Hmisc) #for errorbars\r",
        "\r\n",
        "#first calculate mean, sd, and count\r",
        "\r\n",
        "th.mean<-tapply(data1$NormalizedExp,data1$Treatment, mean)\r",
        "\r\n",
        "th.sd<-tapply(data1$NormalizedExp,data1$Treatment, sd)\r",
        "\r\n",
        "th.count<-tapply(data1$NormalizedExp,data1$Treatment, length)\r",
        "\r\n",
        "#calculate SE\r",
        "\r\n",
        "th.se<-th.sd/(sqrt(th.count))\r",
        "\r\n",
        "#Bind mean and SE for plotting\r",
        "\r\n",
        "threshold.plot<-cbind(th.mean,th.se)\r",
        "\r\n",
        "threshold.plot\r",
        "\r\n",
        "#Data was exported and relative expression and SE were claculated\r",
        "\r\n",
        "\r",
        "\r\n",
        "##Upload plot data\r",
        "\r\n",
        "th.plot <-read.csv('Mean_threshold_plot.csv', header = TRUE)\r",
        "\r\n",
        "head(th.plot )\r",
        "\r\n",
        "\r",
        "\r\n",
        "yplus<-th.plot$Relative+th.plot$RelativeSE\r",
        "\r\n",
        "yminus<-th.plot$Relative\r",
        "\r\n",
        "\r",
        "\r\n",
        "barx<-barplot(th.plot$Relative, las=1, ylab=\"Relative HIF-1a mRNA Expression\",\r",
        "\r\n",
        "  ylim=c(0,3),cex.lab=1.1, cex.axis=1.11,col=\"gray\",  axis.lty=1,\r",
        "\r\n",
        "  names.arg=c(\"Control\",\"Moderate\",\"Hypoxia\"), \r",
        "\r\n",
        "  cex=1.1)\r",
        "\r\n",
        "errbar(barx,th.plot$Relative, yplus,yminus, type='n', add=T)\r",
        "\r\n",
        "\r",
        "\r\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cd /Users/sr320/Desktop/1190951/Lab\\ experiments"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "/Users/sr320/Desktop/1190951/Lab experiments\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "threshold.data <-read.csv('Threshold_data.csv', header = TRUE)\n",
      "head(threshold.data)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "  Sample Treatment Size NormalizedExp\n",
        "1    14C   Control  128   0.000284707\n",
        "2    15C   Control  119   0.000094400\n",
        "3    16C   Control  108   0.000048100\n",
        "4    17C   Control  124   0.000103488\n",
        "5    18C   Control  114   0.000067300\n",
        "6    31C   Control  131   0.001671205\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "th.anova<-lm(NormalizedExp~Treatment + Size, data=threshold.data)\n",
      "summary(th.anova)\n",
      "anova(th.anova)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "Analysis of Variance Table\n",
        "\n",
        "Response: NormalizedExp\n",
        "          Df     Sum Sq    Mean Sq F value  Pr(>F)  \n",
        "Treatment  2 1.0011e-05 5.0057e-06  3.2856 0.04353 *\n",
        "Size       1 5.0000e-09 4.6000e-09  0.0030 0.95637  \n",
        "Residuals 67 1.0208e-04 1.5235e-06                  \n",
        "---\n",
        "Signif. codes:  0 \u2018***\u2019 0.001 \u2018**\u2019 0.01 \u2018*\u2019 0.05 \u2018.\u2019 0.1 \u2018 \u2019 1\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "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": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R \n",
      "library(Hmisc) #for errorbars\n",
      "#first calculate mean, sd, and count\n",
      "th.mean<-tapply(data1$NormalizedExp,data1$Treatment, mean)\n",
      "th.sd<-tapply(data1$NormalizedExp,data1$Treatment, sd)\n",
      "th.count<-tapply(data1$NormalizedExp,data1$Treatment, length)\n",
      "#calculate SE\n",
      "th.se<-th.sd/(sqrt(th.count))\n",
      "#Bind mean and SE for plotting\n",
      "threshold.plot<-cbind(th.mean,th.se)\n",
      "threshold.plot\n",
      "#Data was exported and relative expression and SE were claculated\n",
      "\n",
      "##Upload plot data\n",
      "th.plot <-read.csv('Mean_threshold_plot.csv', header = TRUE)\n",
      "head(th.plot )\n",
      "\n",
      "yplus<-th.plot$Relative+th.plot$RelativeSE\n",
      "yminus<-th.plot$Relative\n",
      "\n",
      "barx<-barplot(th.plot$Relative, las=1, ylab=\"Relative HIF-1a mRNA Expression\",\n",
      "  ylim=c(0,3),cex.lab=1.1, cex.axis=1.11,col=\"gray\",  axis.lty=1,\n",
      "  names.arg=c(\"Control\",\"Moderate\",\"Hypoxia\"), \n",
      "  cex=1.1)\n",
      "errbar(barx,th.plot$Relative, yplus,yminus, type='n', add=T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Error in tapply(data1$NormalizedExp, data1$Treatment, mean) : \n",
        "  object 'data1' not found\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "!cat /Users/sr320/Desktop/1190951/Lab\\ experiments/Mean_threshold_plot.csv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Treament,th.mean,th.se,Relative,RelativeSE\r",
        "\r\n",
        "Control,0.001018586,0.00024783,1,0.243307786\r",
        "\r\n",
        "Moderate,0.000942439,0.000210986,0.925242653,0.20713659\r",
        "\r\n",
        "Hypoxia,0.00177082,0.000290021,1.738507702,0.284729423\r",
        "\r\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Field Data"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "workdir <-\"/Users/sr320/Desktop/1190951/Field\\ data\"\n",
      "setwd(workdir)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "field.data <-read.csv('All_field_qPCR.csv', header = TRUE)\n",
      "head(field.data)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "   Sample Year  Month Site Size TrawlDepth NormalizedExp RelativeExp Elevated\n",
        "1 2012-42 2012 August    A  178       35.6   0.000449324  0.66562636        0\n",
        "2 2012-43 2012 August    A  195       35.6   0.000026700  0.03959079        0\n",
        "3 2012-44 2012 August    A  194       35.6   0.000029800  0.04408860        0\n",
        "4 2012-45 2012 August    A  195       35.6   0.001332981  1.97467024        1\n",
        "5 2012-46 2012 August    B  195       72.5   0.000909842  1.34783478        0\n",
        "6 2012-47 2012 August    B  171       72.5   0.000679925  1.00723675        0\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "library(lme4)\n",
      "#Year model\n",
      "yr.mod<- glm(Elevated ~ factor(Year) + Site + Size + TrawlDepth, family=binomial(link = \"logit\"), data=field.data)\n",
      "summary(yr.mod)\n",
      "anova(yr.mod, test=\"Chisq\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Error in library(lme4) : there is no package called \u2018lme4\u2019\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}