edgeR-tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DE Analysis with EdgeR\n",
    "\n",
    "After acquiring a table/matrix of read counts per gene, we can begin conducting DE analysis to determine genes which are differentially expressed in different tissue populations. However, before we can get our list of DE genes, we need to process our data for usage with edgeR.\n",
    "\n",
    "EdgeR documentation can be found [here](https://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf). This resource has information on the different functions used in this tutorial, as well as a more in-depth explanation of how they work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: limma\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(edgeR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downloading Data\n",
    "\n",
    "First we need some data to analyze. Lets get some data from a [GEO Submission](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60450). Download the GSE60450_Lactation-GenewiseCounts.txt.gz file, unzip it, then load it into R."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rawcounts = read.table(\"GSE60450_Lactation-GenewiseCounts.txt\",header=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<caption>A data.frame: 6 × 14</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>EntrezGeneID</th><th scope=col>Length</th><th scope=col>MCL1.DG_BC2CTUACXX_ACTTGA_L002_R1</th><th scope=col>MCL1.DH_BC2CTUACXX_CAGATC_L002_R1</th><th scope=col>MCL1.DI_BC2CTUACXX_ACAGTG_L002_R1</th><th scope=col>MCL1.DJ_BC2CTUACXX_CGATGT_L002_R1</th><th scope=col>MCL1.DK_BC2CTUACXX_TTAGGC_L002_R1</th><th scope=col>MCL1.DL_BC2CTUACXX_ATCACG_L002_R1</th><th scope=col>MCL1.LA_BC2CTUACXX_GATCAG_L001_R1</th><th scope=col>MCL1.LB_BC2CTUACXX_TGACCA_L001_R1</th><th scope=col>MCL1.LC_BC2CTUACXX_GCCAAT_L001_R1</th><th scope=col>MCL1.LD_BC2CTUACXX_GGCTAC_L001_R1</th><th scope=col>MCL1.LE_BC2CTUACXX_TAGCTT_L001_R1</th><th scope=col>MCL1.LF_BC2CTUACXX_CTTGTA_L001_R1</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>   497097</td><td>3634</td><td>438</td><td>300</td><td> 65</td><td>237</td><td>354</td><td>287</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>100503874</td><td>3259</td><td>  1</td><td>  0</td><td>  1</td><td>  1</td><td>  0</td><td>  4</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>100038431</td><td>1634</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>    19888</td><td>9747</td><td>  1</td><td>  1</td><td>  0</td><td>  0</td><td>  0</td><td>  0</td><td> 10</td><td>  3</td><td> 10</td><td>  2</td><td>  0</td><td>  0</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>    20671</td><td>3130</td><td>106</td><td>182</td><td> 82</td><td>105</td><td> 43</td><td> 82</td><td> 16</td><td> 25</td><td> 18</td><td>  8</td><td>  3</td><td> 10</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>    27395</td><td>4203</td><td>309</td><td>234</td><td>337</td><td>300</td><td>290</td><td>270</td><td>560</td><td>464</td><td>489</td><td>328</td><td>307</td><td>342</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 6 × 14\n",
       "\\begin{tabular}{r|llllllllllllll}\n",
       "  & EntrezGeneID & Length & MCL1.DG\\_BC2CTUACXX\\_ACTTGA\\_L002\\_R1 & MCL1.DH\\_BC2CTUACXX\\_CAGATC\\_L002\\_R1 & MCL1.DI\\_BC2CTUACXX\\_ACAGTG\\_L002\\_R1 & MCL1.DJ\\_BC2CTUACXX\\_CGATGT\\_L002\\_R1 & MCL1.DK\\_BC2CTUACXX\\_TTAGGC\\_L002\\_R1 & MCL1.DL\\_BC2CTUACXX\\_ATCACG\\_L002\\_R1 & MCL1.LA\\_BC2CTUACXX\\_GATCAG\\_L001\\_R1 & MCL1.LB\\_BC2CTUACXX\\_TGACCA\\_L001\\_R1 & MCL1.LC\\_BC2CTUACXX\\_GCCAAT\\_L001\\_R1 & MCL1.LD\\_BC2CTUACXX\\_GGCTAC\\_L001\\_R1 & MCL1.LE\\_BC2CTUACXX\\_TAGCTT\\_L001\\_R1 & MCL1.LF\\_BC2CTUACXX\\_CTTGTA\\_L001\\_R1\\\\\n",
       "  & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int> & <int>\\\\\n",
       "\\hline\n",
       "\t1 &    497097 & 3634 & 438 & 300 &  65 & 237 & 354 & 287 &   0 &   0 &   0 &   0 &   0 &   0\\\\\n",
       "\t2 & 100503874 & 3259 &   1 &   0 &   1 &   1 &   0 &   4 &   0 &   0 &   0 &   0 &   0 &   0\\\\\n",
       "\t3 & 100038431 & 1634 &   0 &   0 &   0 &   0 &   0 &   0 &   0 &   0 &   0 &   0 &   0 &   0\\\\\n",
       "\t4 &     19888 & 9747 &   1 &   1 &   0 &   0 &   0 &   0 &  10 &   3 &  10 &   2 &   0 &   0\\\\\n",
       "\t5 &     20671 & 3130 & 106 & 182 &  82 & 105 &  43 &  82 &  16 &  25 &  18 &   8 &   3 &  10\\\\\n",
       "\t6 &     27395 & 4203 & 309 & 234 & 337 & 300 & 290 & 270 & 560 & 464 & 489 & 328 & 307 & 342\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 6 × 14\n",
       "\n",
       "| <!--/--> | EntrezGeneID &lt;int&gt; | Length &lt;int&gt; | MCL1.DG_BC2CTUACXX_ACTTGA_L002_R1 &lt;int&gt; | MCL1.DH_BC2CTUACXX_CAGATC_L002_R1 &lt;int&gt; | MCL1.DI_BC2CTUACXX_ACAGTG_L002_R1 &lt;int&gt; | MCL1.DJ_BC2CTUACXX_CGATGT_L002_R1 &lt;int&gt; | MCL1.DK_BC2CTUACXX_TTAGGC_L002_R1 &lt;int&gt; | MCL1.DL_BC2CTUACXX_ATCACG_L002_R1 &lt;int&gt; | MCL1.LA_BC2CTUACXX_GATCAG_L001_R1 &lt;int&gt; | MCL1.LB_BC2CTUACXX_TGACCA_L001_R1 &lt;int&gt; | MCL1.LC_BC2CTUACXX_GCCAAT_L001_R1 &lt;int&gt; | MCL1.LD_BC2CTUACXX_GGCTAC_L001_R1 &lt;int&gt; | MCL1.LE_BC2CTUACXX_TAGCTT_L001_R1 &lt;int&gt; | MCL1.LF_BC2CTUACXX_CTTGTA_L001_R1 &lt;int&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| 1 |    497097 | 3634 | 438 | 300 |  65 | 237 | 354 | 287 |   0 |   0 |   0 |   0 |   0 |   0 |\n",
       "| 2 | 100503874 | 3259 |   1 |   0 |   1 |   1 |   0 |   4 |   0 |   0 |   0 |   0 |   0 |   0 |\n",
       "| 3 | 100038431 | 1634 |   0 |   0 |   0 |   0 |   0 |   0 |   0 |   0 |   0 |   0 |   0 |   0 |\n",
       "| 4 |     19888 | 9747 |   1 |   1 |   0 |   0 |   0 |   0 |  10 |   3 |  10 |   2 |   0 |   0 |\n",
       "| 5 |     20671 | 3130 | 106 | 182 |  82 | 105 |  43 |  82 |  16 |  25 |  18 |   8 |   3 |  10 |\n",
       "| 6 |     27395 | 4203 | 309 | 234 | 337 | 300 | 290 | 270 | 560 | 464 | 489 | 328 | 307 | 342 |\n",
       "\n"
      ],
      "text/plain": [
       "  EntrezGeneID Length MCL1.DG_BC2CTUACXX_ACTTGA_L002_R1\n",
       "1    497097    3634   438                              \n",
       "2 100503874    3259     1                              \n",
       "3 100038431    1634     0                              \n",
       "4     19888    9747     1                              \n",
       "5     20671    3130   106                              \n",
       "6     27395    4203   309                              \n",
       "  MCL1.DH_BC2CTUACXX_CAGATC_L002_R1 MCL1.DI_BC2CTUACXX_ACAGTG_L002_R1\n",
       "1 300                                65                              \n",
       "2   0                                 1                              \n",
       "3   0                                 0                              \n",
       "4   1                                 0                              \n",
       "5 182                                82                              \n",
       "6 234                               337                              \n",
       "  MCL1.DJ_BC2CTUACXX_CGATGT_L002_R1 MCL1.DK_BC2CTUACXX_TTAGGC_L002_R1\n",
       "1 237                               354                              \n",
       "2   1                                 0                              \n",
       "3   0                                 0                              \n",
       "4   0                                 0                              \n",
       "5 105                                43                              \n",
       "6 300                               290                              \n",
       "  MCL1.DL_BC2CTUACXX_ATCACG_L002_R1 MCL1.LA_BC2CTUACXX_GATCAG_L001_R1\n",
       "1 287                                 0                              \n",
       "2   4                                 0                              \n",
       "3   0                                 0                              \n",
       "4   0                                10                              \n",
       "5  82                                16                              \n",
       "6 270                               560                              \n",
       "  MCL1.LB_BC2CTUACXX_TGACCA_L001_R1 MCL1.LC_BC2CTUACXX_GCCAAT_L001_R1\n",
       "1   0                                 0                              \n",
       "2   0                                 0                              \n",
       "3   0                                 0                              \n",
       "4   3                                10                              \n",
       "5  25                                18                              \n",
       "6 464                               489                              \n",
       "  MCL1.LD_BC2CTUACXX_GGCTAC_L001_R1 MCL1.LE_BC2CTUACXX_TAGCTT_L001_R1\n",
       "1   0                                 0                              \n",
       "2   0                                 0                              \n",
       "3   0                                 0                              \n",
       "4   2                                 0                              \n",
       "5   8                                 3                              \n",
       "6 328                               307                              \n",
       "  MCL1.LF_BC2CTUACXX_CTTGTA_L001_R1\n",
       "1   0                              \n",
       "2   0                              \n",
       "3   0                              \n",
       "4   0                              \n",
       "5  10                              \n",
       "6 342                              "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head(rawcounts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When analyzing read count data, its important that there are only as many columns as samples. So lets drop the length column and make the geneID column into the table index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill in the code below\n",
    "countdata = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining our groups"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can do any analysis, we have to differentiate the types of groups from each other in our RNA-seq experiment (treatment vs control or liver vs heart vs brain or....). To do so, create a vector of the different groups in order of the columns of your read count dataset.\n",
    "\n",
    "So, if you have a table with two heart columns followed by two liver columns, your group vector could be:\n",
    "\n",
    "`groups = c(\"hrt\",\"hrt\",\"liv\",\"liv\")`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the GEO data, there are 6 groups, according to this figure. \n",
    "\n",
    "![](https://galaxyproject.github.io/training-material/topics/transcriptomics/images/rna-seq-reads-to-counts/mouse_exp.png)\n",
    "\n",
    "\n",
    "However for simplicity's sake, lets assume there are only two groups, basal and luminal. In that case, our group vector should be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "mgroups = rep(c(0,1), each = 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>0</li><li>0</li><li>0</li><li>0</li><li>0</li><li>0</li><li>1</li><li>1</li><li>1</li><li>1</li><li>1</li><li>1</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 0\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\item 1\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 0\n",
       "2. 0\n",
       "3. 0\n",
       "4. 0\n",
       "5. 0\n",
       "6. 0\n",
       "7. 1\n",
       "8. 1\n",
       "9. 1\n",
       "10. 1\n",
       "11. 1\n",
       "12. 1\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " [1] 0 0 0 0 0 0 1 1 1 1 1 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mgroups"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a DGEList\n",
    "\n",
    "EdgeR analyzes data in a specific data object called a DGEList. Lets create one with our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = DGEList(counts=countdata,group=factor(mgroups))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The beauty of DGELists is that renaming the list does not overwrite the original data, just add more elements to it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering the Data\n",
    "\n",
    "Now lets get rid of genes which don't occur very frequently in either group. Let do this in two ways, one using a threshold we define, and one using edgeR's filtering method. \n",
    "\n",
    "For the first, only keep genes that have at least 5 reads in at least 4 samples each.\n",
    "\n",
    "For the second, use the `filterByExpr` function in the edgeR package.\n",
    "\n",
    "Lets save both of these filtered methods, and we can compare our list of DE genes using both methods at the end."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalizing the Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Its important to recognize that raw read count data is NOT directly equal to the in-vivo mRNA levels of transcripts in a population. Because of this, it's important to transform the read counts to more accurately reflect the in-vivo levels. By default, edgeR uses the TMM method to do so. You can read more about how TMM works [here](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2010-11-3-r25) and [here](https://www.frontiersin.org/articles/10.3389/fgene.2016.00164/full).\n",
    "\n",
    "In summary, TMM corrects for:\n",
    "\n",
    "* gene length (the longer the gene, the more reads it will have. Not very relevant for DE analysis)\n",
    "* library size (the more total reads in a sample, the more reads it will have per transcript)\n",
    "* transcriptome size (the more total DE genes in a sample, the less reads will be registered per other transcripts; think of reading power as constant, and having to be divided between number of transcripts in a cel)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets normalize both of our filtered datasets using `calcNormFactors`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill in the code below!\n",
    "y_thresh = calcNormFactors()\n",
    "y_edger = calcNormFactors()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What do these lists look like?\n",
    "y_edger"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "note that the raw values are not changed. Instead, normalization factors are saved, which will be used later in the calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EXTRA: Why not use other normalization methods?\n",
    "\n",
    "In addition to algorithms which correct for all of the above mentioned differences in samples/genes, some are much simpler.\n",
    "\n",
    "* cpm, or counts per million, corrects only for library size\n",
    "* RPKM/FPKM, or reads/fragments per kilobase per million, corrects for library size and gene size\n",
    "\n",
    "Implementations for both of these can be found [here](https://sites.google.com/site/wiki4metagenomics/pdf/definition/rpkm-calculation).\n",
    "\n",
    "Lets do a thought experiment on why we shouldn't use these values in DE analysis. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume you have a sample A with 50 million total reads. How many raw counts would a 40 kb transcript need to have an RPKM value equal to a 10 kb transcript? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even though these raw counts are very different, they still have the same RPKM value, making it unideal for comparing differential expression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the Data\n",
    "\n",
    "Now that we have normalized values, lets check our data to see how different our groups are from each other. We can do this with an MDS (multidimensional scaling) plot. It works via placing similar samples close to each other, and dissimilar samples far apart. Use the `plotMDS` function to generate the graph.\n",
    "\n",
    "Note: this is essentially the same as a PCA of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMDS()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the plot look? Does it appear that samples of the same group are clustered together? How do the different filtering methods look when their plots are compared to each other?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dispersion Estimation\n",
    "\n",
    "Recognizing how variable samples are from each other within a group is necessary for edgeR to properly calculate a list of DE genes. Use the `estimateDisp` function to estimate dispersions for your two DGELists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How does it work?\n",
    "\n",
    "It's difficult to properly measure dispersion between samples within a group when usually there are very few samples in a group for RNA-sequencing. To account for this issue, edgeR makes some assumptions about the distriution of gene read variances and uses genes with similar values and variances to fill in the necessary data to properly estimate dispersion (common dispersion). And to account for any broad assumptions made from calculating common dispersions, tagwise dispersions are also calculated using a [bayesian approach](https://academic.oup.com/bioinformatics/article/23/21/2881/372869)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DE Analysis\n",
    "\n",
    "Now that all our pre-reqs are finished, we can finally conduct DE Analysis. EdgeR uses a glm to generalize model creation for identifying DE genes. It also uses an exact test under specific conditions. Which one should we use in this case? Use it and create a list of DE genes for both DGELists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do these lists differ? Does the FDR differ between them for each gene? Why or why not? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heatmap Generation\n",
    "\n",
    "Congrats, we now have a list of DE genes between our groups! However, its always nice to have a visual representation of the differences in transcript levels between samples and between groups. To do this, we can generate a heatmap to indicate differences between samples. \n",
    "\n",
    "![](https://i.imgur.com/5DJt2wl.png)\n",
    "\n",
    "In this example, each row is a gene and each column is a sample. The data in each row is standardized, so a negative value indicates a value below the \"average reads\" of that row. We can see a strict separation of gene expressions between the two groups."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lets try making our own heatmaps. Create a list of the most differentially expressed genes you identified and create a heatmap. There are a variety of R packages to create heatmaps, including [pheatmap](https://davetang.org/muse/2018/05/15/making-a-heatmap-in-r-with-the-pheatmap-package/) and [heatmap2](https://davetang.org/muse/2010/12/06/making-a-heatmap-with-r/). Heatmaps can also be generated in Python.\n",
    "\n",
    "Remember to use **normalized** read counts in your heatmap (crucial whenever you are comparing samples) and to standardize your values. Though there are a variety of methods to standardize data, I recommend using z-scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does your heatmap look? Do the two groups separate nicely based on the DE genes? How does the heatmap look if you use another normalization method (cpm or rpkm or tmm)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following resources were used to help construct this tutorial:\n",
    "\n",
    "https://galaxyproject.github.io/training-material/topics/transcriptomics/tutorials/rna-seq-counts-to-genes/tutorial.html\n",
    "\n",
    "https://www.youtube.com/watch?v=5tGCBW3_0IA\n",
    "\n",
    "https://web.stanford.edu/class/bios221/labs/rnaseq/lab_4_rnaseq.html\n",
    "\n",
    "http://78.128.216.61:8081/manual/ngs-dea-edger-RNA.html\n",
    "\n",
    "https://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf"
   ]
  }
 ],
 "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": "4.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}