Unit2-RNAseq/3.4_DESeq2_import_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RNAseq 3.4: Importing gene counts and metadata into DESeq2 for expression analysis\n",
    "This jupyter notebook describes the **minimal** import process of read count data into DESeq2. It provides the backbone workflow which can be extended by adding more samples, factors, etc.\n",
    "\n",
    "#### Start by doing `git pull` in the course directory then copying the counts files to your preferred location.\n",
    "#### Locate gene counts files in a dedicated directory (no other files should be present)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>'GSM2580319_counts.txt'</li>\n",
       "\t<li>'GSM2580320_counts.txt'</li>\n",
       "\t<li>'GSM2580323_counts.txt'</li>\n",
       "\t<li>'GSM2580324_counts.txt'</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'GSM2580319\\_counts.txt'\n",
       "\\item 'GSM2580320\\_counts.txt'\n",
       "\\item 'GSM2580323\\_counts.txt'\n",
       "\\item 'GSM2580324\\_counts.txt'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'GSM2580319_counts.txt'\n",
       "2. 'GSM2580320_counts.txt'\n",
       "3. 'GSM2580323_counts.txt'\n",
       "4. 'GSM2580324_counts.txt'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"GSM2580319_counts.txt\" \"GSM2580320_counts.txt\" \"GSM2580323_counts.txt\"\n",
       "[4] \"GSM2580324_counts.txt\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "list.files(\"/mnt/c/Users/Jerry/Applied-Bioinformatics-HW/data/htseq_out/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```R\n",
    "dir_counts <- \"/mnt/c/Users/Jerry/Applied-Bioinformatics-HW/data/htseq_out\"\n",
    "```\n",
    "\n",
    "#### Then read the filenames into a character vector\n",
    "```R\n",
    "counts_files <- list.files(dir_counts)\n",
    "counts_files\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### This step allows you to subset your samples if necessary (in this example we use all files)\n",
    "```R\n",
    "counts_files[1:4]\n",
    "counts_files <- counts_files[1:4]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create data frame for metadata\n",
    "Follow the specification from DESeq2's help pages:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`?DESeqDataSetFromHTSeqCount`  \n",
    "sampleTable \n",
    "\n",
    "for htseq-count: a data.frame with three or more columns. Each row describes one sample. The first column is the sample name, the second column the file name of the count file generated by htseq-count, and the remaining columns are sample metadata which will be stored in colData."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "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\") \n",
    "#Remove the blank columns\n",
    "samplesInfo <- samplesInfo[,-c(1:2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check the sample table (double check this info against the SRA records)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "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>GSM2580319_counts.txt</td><td>GSM2580319_counts.txt</td><td>mock                 </td></tr>\n",
       "\t<tr><td>GSM2580320_counts.txt</td><td>GSM2580320_counts.txt</td><td>ZIKV                 </td></tr>\n",
       "\t<tr><td>GSM2580323_counts.txt</td><td>GSM2580323_counts.txt</td><td>mock                 </td></tr>\n",
       "\t<tr><td>GSM2580324_counts.txt</td><td>GSM2580324_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 GSM2580319\\_counts.txt & GSM2580319\\_counts.txt & mock                   \\\\\n",
       "\t GSM2580320\\_counts.txt & GSM2580320\\_counts.txt & ZIKV                   \\\\\n",
       "\t GSM2580323\\_counts.txt & GSM2580323\\_counts.txt & mock                   \\\\\n",
       "\t GSM2580324\\_counts.txt & GSM2580324\\_counts.txt & ZIKV                   \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "samplename | filename | group | \n",
       "|---|---|---|---|\n",
       "| GSM2580319_counts.txt | GSM2580319_counts.txt | mock                  | \n",
       "| GSM2580320_counts.txt | GSM2580320_counts.txt | ZIKV                  | \n",
       "| GSM2580323_counts.txt | GSM2580323_counts.txt | mock                  | \n",
       "| GSM2580324_counts.txt | GSM2580324_counts.txt | ZIKV                  | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  samplename            filename              group\n",
       "1 GSM2580319_counts.txt GSM2580319_counts.txt mock \n",
       "2 GSM2580320_counts.txt GSM2580320_counts.txt ZIKV \n",
       "3 GSM2580323_counts.txt GSM2580323_counts.txt mock \n",
       "4 GSM2580324_counts.txt GSM2580324_counts.txt ZIKV "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "samplesInfo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load the DESeq2 package"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```R\n",
    "library(\"DESeq2\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import read counts, sample table and specify the design factor (there is only one factor in this example: group)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```R\n",
    "dds1 <- DESeqDataSetFromHTSeqCount(sampleTable = samplesInfo, \n",
    "                                           directory = dir_counts, \n",
    "                                           design = ~ group)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next we perform the main DESeq2 analysis step which includes statistical model fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "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": [
    "dds1_deseq <- DESeq(dds1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "'DESeqDataSet'"
      ],
      "text/latex": [
       "'DESeqDataSet'"
      ],
      "text/markdown": [
       "'DESeqDataSet'"
      ],
      "text/plain": [
       "[1] \"DESeqDataSet\"\n",
       "attr(,\"package\")\n",
       "[1] \"DESeq2\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Formal class 'DESeqDataSet' [package \"DESeq2\"] with 8 slots\n",
      "  ..@ design            :Class 'formula'  language ~group\n",
      "  .. .. ..- attr(*, \".Environment\")=<environment: R_GlobalEnv> \n",
      "  ..@ dispersionFunction:function (q)  \n",
      "  .. ..- attr(*, \"coefficients\")= Named num [1:2] 0.00454 5.57153\n",
      "  .. .. ..- attr(*, \"names\")= chr [1:2] \"asymptDisp\" \"extraPois\"\n",
      "  .. ..- attr(*, \"fitType\")= chr \"parametric\"\n",
      "  .. ..- attr(*, \"varLogDispEsts\")= num 1.08\n",
      "  .. ..- attr(*, \"dispPriorVar\")= num 0.496\n",
      "  ..@ rowRanges         :Formal class 'CompressedGRangesList' [package \"GenomicRanges\"] with 5 slots\n",
      "  .. .. ..@ unlistData     :Formal class 'GRanges' [package \"GenomicRanges\"] with 7 slots\n",
      "  .. .. .. .. ..@ seqnames       :Formal class 'Rle' [package \"S4Vectors\"] with 4 slots\n",
      "  .. .. .. .. .. .. ..@ values         : Factor w/ 0 levels: \n",
      "  .. .. .. .. .. .. ..@ lengths        : int(0) \n",
      "  .. .. .. .. .. .. ..@ elementMetadata: NULL\n",
      "  .. .. .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. .. .. ..@ ranges         :Formal class 'IRanges' [package \"IRanges\"] with 6 slots\n",
      "  .. .. .. .. .. .. ..@ start          : int(0) \n",
      "  .. .. .. .. .. .. ..@ width          : int(0) \n",
      "  .. .. .. .. .. .. ..@ NAMES          : NULL\n",
      "  .. .. .. .. .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. .. .. .. .. ..@ elementMetadata: NULL\n",
      "  .. .. .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. .. .. ..@ strand         :Formal class 'Rle' [package \"S4Vectors\"] with 4 slots\n",
      "  .. .. .. .. .. .. ..@ values         : Factor w/ 3 levels \"+\",\"-\",\"*\": \n",
      "  .. .. .. .. .. .. ..@ lengths        : int(0) \n",
      "  .. .. .. .. .. .. ..@ elementMetadata: NULL\n",
      "  .. .. .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. .. .. ..@ seqinfo        :Formal class 'Seqinfo' [package \"GenomeInfoDb\"] with 4 slots\n",
      "  .. .. .. .. .. .. ..@ seqnames   : chr(0) \n",
      "  .. .. .. .. .. .. ..@ seqlengths : int(0) \n",
      "  .. .. .. .. .. .. ..@ is_circular: logi(0) \n",
      "  .. .. .. .. .. .. ..@ genome     : chr(0) \n",
      "  .. .. .. .. ..@ elementMetadata:Formal class 'DataFrame' [package \"S4Vectors\"] with 6 slots\n",
      "  .. .. .. .. .. .. ..@ rownames       : NULL\n",
      "  .. .. .. .. .. .. ..@ nrows          : int 0\n",
      "  .. .. .. .. .. .. ..@ listData       : Named list()\n",
      "  .. .. .. .. .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. .. .. .. .. ..@ elementMetadata: NULL\n",
      "  .. .. .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. ..@ elementMetadata:Formal class 'DataFrame' [package \"S4Vectors\"] with 6 slots\n",
      "  .. .. .. .. ..@ rownames       : NULL\n",
      "  .. .. .. .. ..@ nrows          : int 58735\n",
      "  .. .. .. .. ..@ listData       :List of 22\n",
      "  .. .. .. .. .. ..$ baseMean                        : num [1:58735] 1000.43 3.36 636.36 468.68 369.45 ...\n",
      "  .. .. .. .. .. ..$ baseVar                         : num [1:58735] 5665.49 7.02 1467.37 771.31 382.05 ...\n",
      "  .. .. .. .. .. ..$ allZero                         : logi [1:58735] FALSE FALSE FALSE FALSE FALSE FALSE ...\n",
      "  .. .. .. .. .. ..$ dispGeneEst                     : num [1:58735] 1.00e-08 7.45e-01 3.68e-03 2.93e-03 1.00e-08 ...\n",
      "  .. .. .. .. .. ..$ dispGeneIter                    : num [1:58735] 11 10 12 36 2 15 3 4 15 5 ...\n",
      "  .. .. .. .. .. ..$ dispFit                         : num [1:58735] 0.0101 1.6648 0.0133 0.0164 0.0196 ...\n",
      "  .. .. .. .. .. ..$ dispersion                      : num [1:58735] 0.00701 1.33737 0.01043 0.01246 0.01338 ...\n",
      "  .. .. .. .. .. ..$ dispIter                        : num [1:58735] 7 10 8 8 6 9 9 10 8 8 ...\n",
      "  .. .. .. .. .. ..$ dispOutlier                     : logi [1:58735] FALSE FALSE FALSE FALSE FALSE FALSE ...\n",
      "  .. .. .. .. .. ..$ dispMAP                         : num [1:58735] 0.00701 1.33737 0.01043 0.01246 0.01338 ...\n",
      "  .. .. .. .. .. ..$ Intercept                       : num [1:58735] 10.05 1.74 9.33 8.86 8.58 ...\n",
      "  .. .. .. .. .. ..$ group_ZIKV_vs_mock              : num [1:58735] -0.1731 0.0238 -0.0275 0.0233 -0.111 ...\n",
      "  .. .. .. .. .. ..$ SE_Intercept                    : num [1:58735] 0.0912 1.3096 0.112 0.1237 0.1294 ...\n",
      "  .. .. .. .. .. ..$ SE_group_ZIKV_vs_mock           : num [1:58735] 0.129 1.844 0.158 0.174 0.183 ...\n",
      "  .. .. .. .. .. ..$ WaldStatistic_Intercept         : num [1:58735] 110.23 1.33 83.3 71.64 66.33 ...\n",
      "  .. .. .. .. .. ..$ WaldStatistic_group_ZIKV_vs_mock: num [1:58735] -1.3403 0.0129 -0.1739 0.1337 -0.6064 ...\n",
      "  .. .. .. .. .. ..$ WaldPvalue_Intercept            : num [1:58735] 0 0.183 0 0 0 ...\n",
      "  .. .. .. .. .. ..$ WaldPvalue_group_ZIKV_vs_mock   : num [1:58735] 0.18 0.99 0.862 0.894 0.544 ...\n",
      "  .. .. .. .. .. ..$ betaConv                        : logi [1:58735] TRUE TRUE TRUE TRUE TRUE TRUE ...\n",
      "  .. .. .. .. .. ..$ betaIter                        : num [1:58735] 2 4 2 2 2 4 3 3 2 2 ...\n",
      "  .. .. .. .. .. ..$ deviance                        : num [1:58735] 43.6 19.5 42.2 40.4 38.3 ...\n",
      "  .. .. .. .. .. ..$ maxCooks                        : logi [1:58735] NA NA NA NA NA NA ...\n",
      "  .. .. .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. .. .. ..@ elementMetadata:Formal class 'DataFrame' [package \"S4Vectors\"] with 6 slots\n",
      "  .. .. .. .. .. .. ..@ rownames       : NULL\n",
      "  .. .. .. .. .. .. ..@ nrows          : int 22\n",
      "  .. .. .. .. .. .. ..@ listData       :List of 2\n",
      "  .. .. .. .. .. .. .. ..$ type       : chr [1:22] \"intermediate\" \"intermediate\" \"intermediate\" \"intermediate\" ...\n",
      "  .. .. .. .. .. .. .. ..$ description: chr [1:22] \"mean of normalized counts for all samples\" \"variance of normalized counts for all samples\" \"all counts for a gene are zero\" \"gene-wise estimates of dispersion\" ...\n",
      "  .. .. .. .. .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. .. .. .. .. ..@ elementMetadata: NULL\n",
      "  .. .. .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. ..@ elementType    : chr \"GRanges\"\n",
      "  .. .. ..@ metadata       : list()\n",
      "  .. .. ..@ partitioning   :Formal class 'PartitioningByEnd' [package \"IRanges\"] with 5 slots\n",
      "  .. .. .. .. ..@ end            : int [1:58735] 0 0 0 0 0 0 0 0 0 0 ...\n",
      "  .. .. .. .. ..@ NAMES          : chr [1:58735] \"ENSG00000000003\" \"ENSG00000000005\" \"ENSG00000000419\" \"ENSG00000000457\" ...\n",
      "  .. .. .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. .. .. ..@ elementMetadata: NULL\n",
      "  .. .. .. .. ..@ metadata       : list()\n",
      "  ..@ colData           :Formal class 'DataFrame' [package \"S4Vectors\"] with 6 slots\n",
      "  .. .. ..@ rownames       : chr [1:4] \"GSM2580319_counts.txt\" \"GSM2580320_counts.txt\" \"GSM2580323_counts.txt\" \"GSM2580324_counts.txt\"\n",
      "  .. .. ..@ nrows          : int 4\n",
      "  .. .. ..@ listData       :List of 2\n",
      "  .. .. .. ..$ group     : Factor w/ 2 levels \"mock\",\"ZIKV\": 1 2 1 2\n",
      "  .. .. .. ..$ sizeFactor: Named num [1:4] 0.841 1.054 1.11 1.02\n",
      "  .. .. .. .. ..- attr(*, \"names\")= chr [1:4] \"GSM2580319_counts.txt\" \"GSM2580320_counts.txt\" \"GSM2580323_counts.txt\" \"GSM2580324_counts.txt\"\n",
      "  .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. ..@ elementMetadata:Formal class 'DataFrame' [package \"S4Vectors\"] with 6 slots\n",
      "  .. .. .. .. ..@ rownames       : NULL\n",
      "  .. .. .. .. ..@ nrows          : int 2\n",
      "  .. .. .. .. ..@ listData       :List of 2\n",
      "  .. .. .. .. .. ..$ type       : chr [1:2] \"input\" \"intermediate\"\n",
      "  .. .. .. .. .. ..$ description: chr [1:2] \"\" \"a scaling factor for columns\"\n",
      "  .. .. .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. .. .. ..@ elementMetadata: NULL\n",
      "  .. .. .. .. ..@ metadata       : list()\n",
      "  .. .. ..@ metadata       : list()\n",
      "  ..@ assays            :Reference class 'ShallowSimpleListAssays' [package \"SummarizedExperiment\"] with 1 field\n",
      "  .. ..$ data: NULL\n",
      "  .. ..and 14 methods.\n",
      "  ..@ NAMES             : NULL\n",
      "  ..@ elementMetadata   :Formal class 'DataFrame' [package \"S4Vectors\"] with 6 slots\n",
      "  .. .. ..@ rownames       : NULL\n",
      "  .. .. ..@ nrows          : int 58735\n",
      "  .. .. ..@ listData       : Named list()\n",
      "  .. .. ..@ elementType    : chr \"ANY\"\n",
      "  .. .. ..@ elementMetadata: NULL\n",
      "  .. .. ..@ metadata       : list()\n",
      "  ..@ metadata          :List of 1\n",
      "  .. ..$ version:Classes 'package_version', 'numeric_version'  hidden list of 1\n",
      "  .. .. ..$ : int [1:3] 1 20 0\n"
     ]
    }
   ],
   "source": [
    "class(dds1_deseq)\n",
    "str(dds1_deseq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<ol class=list-inline>\n",
       "\t<li>58735</li>\n",
       "\t<li>4</li>\n",
       "</ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 58735\n",
       "\\item 4\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 58735\n",
       "2. 4\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] 58735     4"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dim(dds1_deseq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see from the above that `dds1_deseq` is a DESeqDataSet object with a complicated structure but you need no special knowledge of this object class to use DESeq2. \n",
    "If you can reproduce the above output you have successfully imported and fitted your RNA-seq reads to DESeq2.\n",
    "#### From here, many analyses can be performed including overall transcriptome comparisons, variance analysis, principal component analysis, differential expression analysis etc. These will be explored in subsequent classes."
   ]
  }
 ],
 "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
}