Unit2-RNAseq/3.5_DESeq2_expression_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.5 Gene expression analysis using DESeq2\n",
    "\n",
    "#### First use `git` to retrieve the gene counts files (total n=6). The path to files on the Github repo is:\n",
    "`Applied-Bioinformatics/Unit2-RNAseq/data/htseq_out/day5/`  \n",
    "  \n",
    "Then copy to a directory where you can find it (e.g. your homework folder).\n",
    "#### Locate counts files \n",
    "```R\n",
    "dir_counts <- \"/mnt/c/Users/Jerry/Applied-Bioinformatics-HW/data/htseq_out/day5\"\n",
    "counts_files <- list.files(dir_counts)\n",
    "counts_files\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Prepare a sample table in the format required by DESeq2\n",
    "Here we specify the group assignment manually. It is also possible to parse the sample data matrix from GEO, which can be found by following a link on the data series page: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE97919\n",
    "```\n",
    "samplesInfo <- as.data.frame(matrix(ncol=2, nrow=length(counts_files)))  \n",
    "samplesInfo$samplename <- counts_files  \n",
    "samplesInfo$filename <- counts_files  \n",
    "samplesInfo$group <- c(\"mock\", \"ZIKV\", \"mock\", \"ZIKV\", \"mock\", \"ZIKV\") \n",
    "#Remove the blank columns\n",
    "samplesInfo <- samplesInfo[,-c(1:2)]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check the sample table and make any corrections above if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>samplename</th><th scope=col>filename</th><th scope=col>group</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td>GSM2580321_counts.txt</td><td>GSM2580321_counts.txt</td><td>mock                 </td></tr>\n",
       "\t<tr><td>GSM2580322_counts.txt</td><td>GSM2580322_counts.txt</td><td>ZIKV                 </td></tr>\n",
       "\t<tr><td>GSM2580325_counts.txt</td><td>GSM2580325_counts.txt</td><td>mock                 </td></tr>\n",
       "\t<tr><td>GSM2580326_counts.txt</td><td>GSM2580326_counts.txt</td><td>ZIKV                 </td></tr>\n",
       "\t<tr><td>GSM2580329_counts.txt</td><td>GSM2580329_counts.txt</td><td>mock                 </td></tr>\n",
       "\t<tr><td>GSM2580330_counts.txt</td><td>GSM2580330_counts.txt</td><td>ZIKV                 </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|lll}\n",
       " samplename & filename & group\\\\\n",
       "\\hline\n",
       "\t GSM2580321\\_counts.txt & GSM2580321\\_counts.txt & mock                   \\\\\n",
       "\t GSM2580322\\_counts.txt & GSM2580322\\_counts.txt & ZIKV                   \\\\\n",
       "\t GSM2580325\\_counts.txt & GSM2580325\\_counts.txt & mock                   \\\\\n",
       "\t GSM2580326\\_counts.txt & GSM2580326\\_counts.txt & ZIKV                   \\\\\n",
       "\t GSM2580329\\_counts.txt & GSM2580329\\_counts.txt & mock                   \\\\\n",
       "\t GSM2580330\\_counts.txt & GSM2580330\\_counts.txt & ZIKV                   \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "samplename | filename | group | \n",
       "|---|---|---|---|---|---|\n",
       "| GSM2580321_counts.txt | GSM2580321_counts.txt | mock                  | \n",
       "| GSM2580322_counts.txt | GSM2580322_counts.txt | ZIKV                  | \n",
       "| GSM2580325_counts.txt | GSM2580325_counts.txt | mock                  | \n",
       "| GSM2580326_counts.txt | GSM2580326_counts.txt | ZIKV                  | \n",
       "| GSM2580329_counts.txt | GSM2580329_counts.txt | mock                  | \n",
       "| GSM2580330_counts.txt | GSM2580330_counts.txt | ZIKV                  | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  samplename            filename              group\n",
       "1 GSM2580321_counts.txt GSM2580321_counts.txt mock \n",
       "2 GSM2580322_counts.txt GSM2580322_counts.txt ZIKV \n",
       "3 GSM2580325_counts.txt GSM2580325_counts.txt mock \n",
       "4 GSM2580326_counts.txt GSM2580326_counts.txt ZIKV \n",
       "5 GSM2580329_counts.txt GSM2580329_counts.txt mock \n",
       "6 GSM2580330_counts.txt GSM2580330_counts.txt ZIKV "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samplesInfo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load DESeq2 package\n",
    "```\n",
    "library(\"DESeq2\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import data into R\n",
    "```\n",
    "dds1 <- DESeqDataSetFromHTSeqCount(sampleTable = samplesInfo, \n",
    "                                           directory = dir_counts, \n",
    "                                           design = ~ group)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check the data imported correctly\n",
    "```R\n",
    "colData(dds1)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "estimating size factors\n",
      "estimating dispersions\n",
      "gene-wise dispersion estimates\n",
      "mean-dispersion relationship\n",
      "final dispersion estimates\n",
      "fitting model and testing\n"
     ]
    }
   ],
   "source": [
    "#### Run the main DESeq normalization and analysis function\n",
    "dds1_deseq <- DESeq(dds1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run additional DESeq2 normalization functions\n",
    "```\n",
    "rld <- rlog(dds1_deseq)\n",
    "vsd <- varianceStabilizingTransformation(dds1_deseq)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot a dendrogram of sample transcriptomes\n",
    "```\n",
    "sampleDists <- dist(t(assay(vsd)))\n",
    "plot(hclust(sampleDists))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Principal component plot\n",
    "```\n",
    "plotPCA(rld, intgroup=\"group\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercises\n",
    "1) Repeat the hierarchical clustering using rlog-normalized counts instead of vsd.  \n",
    "2) Examine the normalized counts using `head(assay(vsd))` and `head(assay(rld))`. Now find the corresponding raw counts and plot raw vs normalized counts for the first 6 genes of the first sample. HINT: It is possible to extract elements from `dds1`, including raw counts. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}