NLM_Reproducibility_Drosophila.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reproducibility\n",
    "\n",
    "## Paper\n",
    "\n",
    "### “Drosophila anti-nematode and antibacterial immune regulators revealed by RNA-Seq,” BMC Genomics \n",
    "### https://doi.org/10.1186/s12864-015-1690-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "# get public key and log on to AWS\n",
    "\n",
    "##Get Data\n",
    "\n",
    "#RNA-seq Data\n",
    "#https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token-szetgcuwtvsvhyf&acc=GSE61466\n",
    "\n",
    "# Sample: XGWIE_Symbiotic30h\n",
    "fastq-dump -I --split-files -O datasets --gzip SRR1574318\n",
    "\n",
    "# Sample: XGWIE_Uninfectedcontrol\n",
    "fastq-dump -I --split-files -O datasets --gzip SRR1573961\n",
    "\n",
    "# Reference Genome\n",
    "#ftp://ftp.ensembl.org/pub/release-78/fasta/drosophila_melanogaster/dna/\n",
    " #wget -r --no-parent ftp://ftp.ensembl.org/pub/release-78/fasta/drosophila_melanogaster/dna/\n",
    "wget ftp://ftp.ensembl.org/pub/release-78/fasta/drosophila_melanogaster/dna/Drosophila_melanogaster.BDGP5.dna.toplevel.fa.gz\n",
    "    \n",
    "# .gtf file\n",
    "wget ftp://ftp.ensembl.org/pub/release-78/gtf/drosophila_melanogaster/Drosophila_melanogaster.BDGP5.78.gtf.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "##Get Software\n",
    "# Bioconda\n",
    "# Steps to set up bioconda\n",
    "# 1. Download miniconda\n",
    "wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh\n",
    "chmod +x Miniconda3-latest-Linux-x86_64.sh\n",
    "./Miniconda3-latest-Linux-x86_64.sh\n",
    "\n",
    "# 2. log out the server\n",
    "# 3. set up channels\n",
    "conda config --add channels defaults\n",
    "conda config --add channels bioconda\n",
    "conda config --add channels conda-forge\n",
    "\n",
    "## to install packages\n",
    "conda install -y r\n",
    "\n",
    "## to search for packages\n",
    "conda search bwa\n",
    "\n",
    "## Reference: https://bioconda.github.io\n",
    "\n",
    "# Install SRAtoolkit, HTseq, R, DEseq, cufflinks from Bioconda\n",
    "conda install -y sra-tools htseq r deseq cufflinks\n",
    "\n",
    "#FastQC\n",
    "conda install -y fastqc #(Do we need fastqc??)\n",
    "\n",
    "#Bowtie\n",
    "wget https://sourceforge.net/projects/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip\n",
    "unzip bowtie2-2.3.4.3-linux-x86_64.zip\n",
    "sudo mv bowtie2-2.3.4.3-linux-x86_64 /usr/local/bin\n",
    "PATH=$PATH:/usr/local/bin/bowtie2-2.3.4.3-linux-x86_64\n",
    "export PATH\n",
    "\n",
    "rm bowtie2-2.3.4.3-linux-x86_64.zip\n",
    "\n",
    "#Tophat\n",
    "wget http://ccb.jhu.edu/software/tophat/downloads/tophat-2.1.1.Linux_x86_64.tar.gz\n",
    "tar xvfz tophat-2.0.0.Linux_x86_64.tar.gz\n",
    "sudo mv tophat-2.0.0.Linux_x86_64 /usr/local/bin\n",
    "PATH=$PATH:/usr/local/bin/tophat-2.0.0.Linux_x86_64\n",
    "export PATH\n",
    "\n",
    "rm tophat-2.0.0.Linux_x86_64.tar.gz\n",
    "# Sam tools\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Alignment reads and coverage analysis\n",
    "**requires that gtf file be unzipped and tophat, bowtie, and samtools be installed**\n",
    "***\n",
    "\n",
    "### Generate bowtie index files for Drosophila genome "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "bowtie2-build --large-index \\\n",
    "~/datasets/Drosophila/Drosophila_melanogaster.BDGP5.dna.toplevel.fa.gz \\\n",
    "~/datasets/Drosophila_melanogaster.BDGP5.dna.toplevel "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform alignment using genome index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "  \n",
    "tophat -r 101 \\\n",
    "       --segment-mismatches 2 \\\n",
    "       --segment-length 30 \\\n",
    "       -x 25 \\\n",
    "       -G ~/datasets/Drosophila/Drosophila_melanogaster.BDGP5.78.gtf \\\n",
    "       -p 2 ~/datasets/Drosophila_melanogaster.BDGP5.dna.toplevel \\\n",
    "       ~/datasets/SRR1573961_1.fastq.gz ~/datasets/SRR1573961_2.fastq.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Differential gene expression analysis\n",
    "\n",
    "#Calculate RPKM\n",
    "#R Version 3.3.1\n",
    "\n",
    "#Generate data frame with Data\n",
    "Gene_Name = c(\"A\",\"B\",'C','D')\n",
    "Length_in_KB = c(2,4,1,10)\n",
    "Rep1_Counts = c(10, 20, 5, 0)\n",
    "Rep2_Counts = c(12, 25, 8, 0)\n",
    "Rep3_Counts = c(30, 60, 15, 1)\n",
    "df = data.frame(Gene_Name,Length_in_KB,Rep1_Counts,Rep2_Counts,Rep3_Counts)\n",
    "\n",
    "#Find RPM and RPKM\n",
    "scalingFactor = colSums(df[3:5])/10^6 \n",
    "RPM = cbind(df[,1:2],t(t(df[3:5])/scalingFactor))\n",
    "RPKM = cbind(RPM[1:2],RPM[,3:5]/RPM[,2])\n",
    "\n",
    "#based on https://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Transcript analysis using CUFFLINKS\n",
    "\n",
    "\n",
    "%%bash\n",
    "cufflinks -g ~/datasets/Drosophila/Drosophila_melanogaster.BDGP5.78.gtf -b ~/datasets/Drosophila_melanogaster.BDGP5.dna.toplevel.fa ~/tophat_SRR1574318/accepted_hits.bam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Differential transcript analysis using Cuffdiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Gene ontology analysis (GO)\n",
    "\n",
    "#Website\n",
    "#\"https://david.ncifcrf.gov/summary.jsp\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Differential gene expression analysis using DESeq and General Linear Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# generate read counts using HTSeq\n",
    "htseq-count --stranded=no -f bam Ax_12h_accepted_hits.bam \\ \n",
    "                                 Ax_30h_accepted_hits.bam \\ \n",
    "                                 Sym_12h_accepted_hits.bam \\ \n",
    "                                 Sym_30h_accepted_hits.bam \\ \n",
    "                                 Ph_12h_accepted_hits.bam \\ \n",
    "                                 Ph_30h_accepted_hits.bam \\ \n",
    "                                 Drosophila_melanogaster.BDGP5.78.gtf.gtf > result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "require(DESeq)\n",
    "require(Biobase)\n",
    "\n",
    "# Differential expression analysis\n",
    "countTable <- read.table('result', header=FALSE, row.names=1)\n",
    "colnames(countTable) <- c('Ax_12h', 'Ax_30h', 'Sym_12h', 'Sym_30h', 'Ph_12h', 'Ph_30h')\n",
    "condition <- c('Ax', 'Ax', 'Sym', 'Sym', 'Ph', 'Ph')\n",
    "condition <- factor(condition)\n",
    "# get the countDataSet instance\n",
    "cds <- newCountDataSet(countTable, condition)\n",
    "# normalization\n",
    "cds <- estimateSizeFactors(cds)\n",
    "# variance estimation \n",
    "cds <- estimateDispersions(cds, method='pooled', sharingMode='fit-only')\n",
    "# differential expression tests to get log2 fold changes\n",
    "Ax_Vs_Sym <- nbinomTest(cds, 'Ax', 'Sym')\n",
    "Sym_Vs_Ph <- nbinomTest(cds, 'Sym', 'Ph')\n",
    "Ax_Vs_Ph <- nbinomTest(cds, 'Ax', 'Ph')\n",
    "\n",
    "# GLM tests\n",
    "# model A\n",
    "countTable_A <- read.table('result', header=FALSE, row.names=1)[:,c(1:6)]\n",
    "colnames(countTable_A) <- c('Ax_12h', 'Ax_30h', 'Sym_12h', 'Sym_30h', 'Ph_12h', 'Ph_30h')\n",
    "condition <- c('Ax', 'Ax', 'Sym', 'Sym', 'Ph', 'Ph')\n",
    "condition <- factor(condition)\n",
    "# get the countDataSet instance\n",
    "cds_A <- newCountDataSet(countTable_A, condition)\n",
    "# normalization\n",
    "cds_A <- estimateSizeFactors(cds_A)\n",
    "# variance estimation \n",
    "cds_A <- estimateDispersions(cds_A, method='pooled', sharingMode='fit-only')\n",
    "# model fitting and testing\n",
    "fit1_A <- fitNbinomGLMs(cds_A, count ~ condition)\n",
    "fit0_A <- fitNbinomGLMs(cds_A, count ~ 1)\n",
    "pvalsGLM_A <- nbinomGLMTest(fit1_A, fit0_A)\n",
    "padjGLM_A = p.adjust(pvalsGLM, method='BH')\n",
    "\n",
    "# model B\n",
    "countTable_B <- read.table('result', header=FALSE, row.names=1)[:,c(1,2,5,6)]\n",
    "colnames(countTable_B) <- c('Ax_12h', 'Ax_30h', 'Ph_12h', 'Ph_30h')\n",
    "condition <- c('Ax', 'Ax', 'Ph', 'Ph')\n",
    "condition <- factor(condition)\n",
    "# get the countDataSet instance\n",
    "cds_B <- newCountDataSet(countTable_B, condition)\n",
    "# normalization\n",
    "cds_B <- estimateSizeFactors(cds_B)\n",
    "# variance estimation \n",
    "cds_B <- estimateDispersions(cds_B, method='pooled', sharingMode='fit-only')\n",
    "# model fitting and testing\n",
    "fit1_B <- fitNbinomGLMs(cds_B, count ~ condition)\n",
    "fit0_B <- fitNbinomGLMs(cds_B, count ~ 1)\n",
    "pvalsGLM_B <- nbinomGLMTest(fit1_B, fit0_B)\n",
    "padjGLM_B = p.adjust(pvalsGLM, method='BH')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Docker instructions -- Thanks to Steve Tsang\n",
    "#### To run the Docker image "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "docker run -it -v `pwd`:`pwd` -w `pwd` flytools /bin/bash"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}