2024_MergedRNAseqWorkshop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to RNAseq: experimental analysis and design \n",
    "\n",
    "\n",
    "# Abstract\n",
    "\n",
    "This workshop introduces RNA-seq analysis using R, including the (Ranged)SummarizedExperiment data class, exploratory data analysis and visualizations, data transformations, and experimental design consideratons.\n",
    "\n",
    "This workshop uses quotes and materials from [RNA-seq workflow: gene-level exploratory analysis and differential expression](https://www.bioconductor.org/help/workflows/rnaseqGene/) by Love, Anders, Kim, and Huber and the [Harvard Chan Bioinformatics Core (HBC)](https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html). The notebook is forked and modified from [two notebooks originally created by Levi Waldron](https://github.com/compbiocore/bioconductor-workshop-2).     \n",
    "\n",
    "If you need additional help with RNA-seq analysis, experimental design, or setting up your design matrix, the Computational Biology Core hosts weekly office hours. Please see the [CCV Calendar for an up-to-date schedule](https://events.brown.edu/ccv/week).     \n",
   "\n",
    "# Outline\n",
    "\n",
    "1. Introduction to RNA-seq Data\n",
    "    * What is RNAseq?\n",
    "        * Some more resources about sequencing technologies\n",
    "    * Computational workflow overview\n",
    "2. Standard data types\n",
    "    * Counts matrices\n",
    "    * SummarizedExperiment\n",
    "        * Creating a SummarizedExperiment object from scratch\n",
    "3. Exploring the Airway dataset\n",
    "    * Creating a DESeqDataSet\n",
    "    * Filtering genes not satisfying a minimum expression threshold\n",
    "    * Data transformation\n",
    "    * Exploring data at the sample level\n",
    "        * Boxplots\n",
    "        * Correlation plots\n",
    "        * Density plots\n",
    "        * Sample distances\n",
    "        * PCA\n",
    "        * MDS\n",
    "4. Differential Expression Analysis\n",
    "    * Running DESeq2\n",
    "        * Normalization\n",
    "        * Dispersion Estimation\n",
    "        * Hypothesis Testing\n",
    "    * Building the Results Table\n",
    "        * Statistical Terminology\n",
    "        * Multiple testing\n",
    "        * Other comparisons\n",
    "5. Diagnostic Plots\n",
    "    * MA Plot\n",
    "    * Volcano Plot\n",
    "6. Model Building and Analytic Approaches\n",
    "   * Model formulae and terminology\n",
    "   * The Design Matrix\t\n",
    "   * Contrasts\n",
    "7. Experimental Design   \n",
    "   * Single factor design\n",
    "   * Multi-factor design: two factors\n",
    "   * Blocking\n",
    "   * Interactions\n",
    "   * Repeat Measures/Time Series\n",
    "   * Likelihood ratio test\n",
    "   * Other R packages for differential expression analysis\t\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#  Introduction to RNAseq\n",
    "\n",
    "## What is RNAseq?\n",
    "\n",
    "The central dogma of biology refers to the concept that in general, biological information flows from DNA to RNA to proteins, which are the molecules that encode for various enzymes etc., that help a living organism maintain homeostasis and self-regulate. For this workshop, we are going to talk about analyzing transcriptome data, where genes are regulated at the level of DNA being transcribed into RNA.\n",
    "\n",
    "The general experimental process involves extracting RNA from your different experimental conditions, performing the library preparation steps (e.g., cDNA synthesis, fragmentation, and adding adapters). These libraries are sequenced on an Illumina FlowCell. The output from this process is these raw reads (.fastq or .fq files), which you can then align to the appropriate reference genome or transcriptome, where more reads will align to more highly expressed genes. This data can be processed to get back a counts matrix which will let you look at differential expression and perform other downstream analyses.\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/9/97/Journal.pcbi.1004393.g002.png/2560px-Journal.pcbi.1004393.g002.png\">\n",
    "\n",
    "### Some more resources about sequencing technologies:\n",
    "+ https://www.illumina.com/content/dam/illumina-marketing/documents/products/illumina_sequencing_introduction.pdf      \n",
    "+ Rory Stark, Marta Grezlak, and James Hadfield. NA sequencing: the teenage years Nature Reviews Genetics *Nature Reviews Genetics* **20,** 631-656 (2019)      \n",
    "\n",
    "## Computational workflow overview\n",
    "\n",
    "The computational analysis of an RNA-seq experiment begins from the raw FASTQ files that contain the nucleotide sequence of each read and a quality score at each position. These reads are **aligned** to a reference genome and the alignment files (.SAM or .BAM file) are turned into a counts matrix by programs such as HTSeq or featureCounts. Some newer and very fast alignment-free approaches can also estimate counts per transcript **without alignment**. Instead of directly aligning RNA-seq reads to a reference genome, they perform \"pseudo-alignment\" by assigning reads to transcripts that contain compatible k-mers and generate a table of read counts without the need for a counting step.\n",
    "\n",
    "<img src=\"Images/rnaseq_workflow.png\">\n",
    "\n",
    "There are many different programs for running alignment-based quantification. Each alignment program will make different assumptions about the data and will show some variation in performance. For example, not all aligners are 'splice-aware', which means that they will recognize intron-sized gaps that are spliced from transcripts. Some popular alignment-based tools include [STAR](https://github.com/alexdobin/STAR), [HISAT2](https://daehwankimlab.github.io/hisat2/), [GSNAP/GMAP](http://research-pub.gene.com/gmap/), among many others. Counts can be generated from SAM/BAM alignments using programs such as [featureCounts](http://bioinf.wehi.edu.au/featureCounts/) or [HTSeq](https://htseq.readthedocs.io/en/master/). \n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><b>Note:</b> The information you get back from approaches where you are aligning to a transcriptome (like Salmon and Kallisto) will be dependent on the quality of the annotations, which are not always great (especially for non-model organisms). </div>\n",
    "\n",
    "Alignment-free tools have gained popularity for their fast performance and low memory requirements while maintaining good accuracy.  The most popular of these tools are [Salmon](https://combine-lab.github.io) and [Kallisto](https://pachterlab.github.io/kallisto/). \n",
    "\n",
    "The differential expression analyses introduced here assume you are analyzing a raw counts table, where the rows are the genes or transcripts and the columns are your samples. It is always a good idea to read the documentation for any of these programs to make sure you understand what file formats should be used. The two most popular programs for  differential expression analyses are probably [DESeq2](http://bioconductor.org/packages/DESeq2) and [edgeR](http://bioconductor.org/packages/). Each program makes different assumptions and has different approaches to normalizing the data, and they will likely give you slightly different numbers of differentially expressed genes if you run both programs. \n",
    "\n",
    "The counts matrices can also be transformed so that they are compatible with common statistical methods for exploratory analysis of multidimensional data, for example clustering and principal components analysis (PCA), which work best for data that have the same variance at different mean values, i.e. the mean and variance are independent.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><b>Note:</b> The rest of this workshop assumes you are starting with the output of some sequence analysis pipeline (i.e. counts). </div> \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Standard data types\n",
    "\n",
    "## Counts matrices\n",
    "The counts matrix is a simple, standard way to represent RNAseq experiments. The rows correspond to genes, the columns correspond to each sample, and the data are the counts.\n",
    "\n",
    "<img src=\"Images/counts_matrix.png\" width=\"300\">\n",
    " \n",
    "## SummarizedExperiment\n",
    "\n",
    "The SummarizedExperiment is a data class for representing RNA-seq data. It builds on a simple counts matrix to include information about samples. It can also optionally include information about exons and genes, the experimental design, and additional metadata. It is becoming a very popular and standard way to represent sequencing data and metadata in R.\n",
    "**The component parts of a *SummarizedExperiment* object `se` are:** \n",
    "\n",
    "* `assay(se)` or `assays(se)$counts` contains the matrix of counts\n",
    "* `colData(se)` may contain data about the columns, e.g. patients or biological units\n",
    "* `rowData(se)` may contain data about the rows, e.g. genes or transcript\n",
    "* `rowRanges(se)` may contain genomic ranges for the genes/transcripts\n",
    "* `metadata(se)` may contain information about the experiment\n",
    "\n",
    "<img src=\"Images/summarizedexperiment.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the `airway`,  `SummarizedExperiment`, and other packages we will need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suppressPackageStartupMessages(library('airway'))\n",
    "suppressPackageStartupMessages(library('SummarizedExperiment'))\n",
    "suppressPackageStartupMessages(library('DESeq2'))\n",
    "suppressPackageStartupMessages(library('tidyr'))\n",
    "suppressPackageStartupMessages(library('ggplot2'))\n",
    "suppressPackageStartupMessages(library('patchwork'))\n",
    "suppressPackageStartupMessages(library('vsn'))\n",
    "suppressPackageStartupMessages(library('GGally'))\n",
    "suppressPackageStartupMessages(library('pheatmap'))\n",
    "suppressPackageStartupMessages(library('RColorBrewer'))\n",
    "suppressPackageStartupMessages(library('PoiClaClu'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the purposes of this tutorial, we will be using the `airway` data. This data is in the `airway` package, stored as a `RangedSummarizedExperiment` object -- a type of `SummarizedExperiment`. Four human airway smooth muscle cell lines were left untreated or treatreated with dexamethasone for 18 hours (Himes et al. 2014)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(airway)\n",
    "airway"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `assays` function to see the names of the assays stored in the object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assays(airway)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `$` accessor to view the `counts` assay. We will combine that with `head` and `assay` functions to access the first view rows of the counts matrix. The rows are genes, the columns are the samples, and the data are the counts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(assays(airway)$counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `colData` function to see more information about each sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(colData(airway))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a SummarizedExperiment object from scratch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can create `SummarizedExperiment` objects using the `SummarizedExperiment` function. The `assay` and `colData` arguments are required. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "airway_counts <- assay(airway)\n",
    "airway_samples <- colData(airway)\n",
    "new_se <- SummarizedExperiment(assays = airway_counts, \n",
    "                               colData = airway_samples)\n",
    "new_se"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use `?` to get more information about the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "?SummarizedExperiment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"><b>Exercise 1:</b> Based on what we have just discussed, how would you create a new SummarizedExperiment object that contains information about the genes in the airway data (e.g., the rows) and their location in the genome? What about additional metadata? </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploring the airway dataset\n",
    "\n",
    "Exploratory data analysis of RNAseq data is an important step for assessing data quality. You can look at the expression of genes of interest or relationships/similarities among sample groupings. Filtering and exploring your data is an iterative process; as you explore and visualize your data you can tweak your filtering parameters to see how the filtering impacts the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a `DESeqDataSet`\n",
    "\n",
    "Once we have our fully annotated *SummarizedExperiment* object,\n",
    "we can construct a *DESeqDataSet* object from it that will then form\n",
    "the starting point of the exploratory analysis. Here, `~ cell + dex` means that we are concerned with two covariates, the cell line (`cell`) and whether or not the samples were treated with dexamethasone (`dex`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dds <- DESeqDataSet(airway, design = ~ cell + dex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the `$` accessor to access the individual columns (like you would with a data frame). For example, we can look at the factor levels of our dex treatments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levels(colData(dds)$dex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can set `untrt` as the reference factor level using `relevel`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colData(dds)$dex <- relevel(colData(dds)$dex, ref = \"untrt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levels(colData(dds)$dex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also set all the factor levels more explicitly using `factor`  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colData(dds)$dex  <- factor(colData(dds)$dex, levels = c(\"untrt\",\"trt\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levels(colData(dds)$dex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering genes not satisfying a minimum expression threshold\n",
    "\n",
    "Our count matrix might contain rows with only zeroes or genes that are very lowly expressed in every condition tested. These genes are not likely to be informative, so to reduce the size of the object we can remove these rows.\n",
    "Here, we remove rows of the *DESeqDataSet* that\n",
    "have no counts, or only a single count across **all** samples. First, we can use `nrow` to see how many genes are in the dataset before we filter:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrow(dds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's re-assign `dds` to `dds.original` before we do any filtering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dds.original <- dds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the first few lines using `head`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(counts(dds.original))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's filter the expression matrix so we are including genes that have at least one fragment mapped in at least 1 sample (or `rowSums` greater than 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dds <- dds[ rowSums(counts(dds)) > 1, ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the first few rows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(counts(dds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the gene that had zero counts in every sample (`ENSG00000000005`) is gone, but `ENSG00000000938` remains because it has a rowSum greater than 1. How many genes do we have left after filtering?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrow(dds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only 46% of the genes remain after performing this filtering step. A more stringent method of filtering, such as removing all genes that have any 0 count entries at all, would yield an even smaller percentage of genes. Let's look at the count densities from before and after filtering:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we can use the `pivot_longer` function to turn our data into a longer format ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_melted_original <-  as.data.frame(counts(dds.original)) %>% \n",
    "pivot_longer(\n",
    "    cols = starts_with(\"SRR\"),\n",
    "    values_to = \"counts\",\n",
    "    names_to = \"library\",\n",
    "    values_drop_na = TRUE)\n",
    "\n",
    "counts_melted_filtered <- as.data.frame(counts(dds)) %>% \n",
    "pivot_longer(\n",
    "    cols = starts_with(\"SRR\"),\n",
    "    values_to = \"counts\",\n",
    "    names_to = \"library\",\n",
    "    values_drop_na = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can compare the density plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 <- ggplot(counts_melted_original, aes(x = counts)) + \n",
    "geom_density(aes(group = library, color = library), adjust = 5) + \n",
    "xlim(0, 200) + \n",
    "ylim(0, .6) +\n",
    "ggtitle('original')\n",
    "\n",
    "\n",
    "p2 <- ggplot(counts_melted_filtered, aes(x = counts)) + \n",
    "geom_density(aes(group = library, color = library), adjust = 5) + \n",
    "xlim(0, 200) + \n",
    "ylim(0, .6) + \n",
    "ggtitle('filtered')\n",
    "\n",
    "\n",
    "options(repr.plot.height=5, repr.plot.width=10)\n",
    "p1 | p2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The warning message just means that there are some data points that are not displayed because they are not within the axis ranges we have defined. We can see that our filtering has shifted our density plots -- before filtering nearly all of the area under the curve was near 0. These genes do not contribute meaningful information to our analysis (and can even bias some of our analytical techniques) so we remove them.\n",
    "\n",
    "**Note**: For differential expression analysis later, filtering is \n",
    "allowable but not necessary if using [Independent Hypothesis Weighting](http://www.bioconductor.org/packages/IHW),\n",
    "as is the default behavior of `DESeq2`. Independent hypothesis weighting (IHW) is a multiple testing procedure that increases power compared to the method of Benjamini and Hochberg by assigning data-driven weights to each hypothesis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"><b>Exercise 2:</b> How many genes pass a different filtering cutoff where at least 3 samples must have a count of 10 or higher? How does this filtering change the density plots? </div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Transformation\n",
    "\n",
    "Why transform the data? Many common statistical methods for *exploratory* analysis of\n",
    "multidimensional data, for example clustering and *principal\n",
    "components analysis* (PCA), **work best for data that generally has the\n",
    "same variance across different values of the mean**. When\n",
    "the expected amount of variance is approximately the same across\n",
    "different mean values, the data is said to be *homoskedastic*.      \n",
    "\n",
    "For RNA-seq count data, however, the expected variance grows with the mean. The problem with this mean-variance relationship is that **if one performs PCA directly on a matrix of counts or normalized counts (e.g. correcting for differences in sequencing depth), the resulting plot typically depends mostly on the genes with the *highest* counts because they show the largest absolute differences between samples**. Simply put, the genes with the highest counts (i.e., those that are more highly expressed) will dominate the analysis simply due to them having more variability.      \n",
    "\n",
    "A simple and often used strategy to address this problem is to take the logarithm of the normalized\n",
    "count values plus a pseudocount of 1 (this pseudocount of 1 accounts for 0-count observations, as the log of 0 is undefined and the log of 1 is always 0, irrespective of base); however, this approach can also introduce problems; depending on the choice of pseudocount, the genes with **the very *lowest* counts\n",
    "will now contribute a great deal of noise to the resulting plot, because\n",
    "taking the logarithm of small counts actually inflates their variance**.     \n",
    "\n",
    "As a solution, DESeq2 offers two additional transformations for count data that stabilize the variance across the mean: the **regularized-logarithm transformation or rlog** (Love, Huber, Anders, Genome Biology 2014); and the **variance stabilizing transformation or VST** for negative binomial data with a dispersion-mean trend (Anders and Huber, Genome Biology 2010)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare these different transformation methods. First, we will run the `rlog` and `vst` transformations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rld <- rlog(dds, blind = FALSE) #rlog\n",
    "vsd <- vst(dds, blind = FALSE) #vst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also run the standard `log2` transformation, where a pseudocount of 1 is added to each point before transforming. We need to first estimate *size factors* to\n",
    "account for sequencing depth, and then specify `normalized=TRUE`.\n",
    "Sequencing depth correction is done automatically for the `rlog`\n",
    "and the `vst`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsESF <- estimateSizeFactors(dds) #estimate size factors first\n",
    "logt <- data.frame(log2(counts(ddsESF, normalized=TRUE)[, 1:2] + 1)) #standard log2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can look at some scatter plots to compare the raw counts to the three different transformation methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 <- ggplot(data.frame(assay(dds)[, 1:2]), aes(x = SRR1039508, y = SRR1039509)) + \n",
    "geom_point() + \n",
    "geom_abline(intercept = 0, slope = 1, color = 'blue', linewidth = 2) +\n",
    "coord_fixed() +\n",
    "xlim(0,300000)+\n",
    "ylim(0,300000) +\n",
    "ggtitle('counts')\n",
    "\n",
    "p2 <- ggplot(logt, aes(x = SRR1039508, y = SRR1039509)) + \n",
    "geom_point() + \n",
    "geom_abline(intercept = 0, slope = 1, color = 'blue', linewidth = 2) +\n",
    "coord_fixed() +\n",
    "xlim(0,19)+\n",
    "ylim(0,19) +\n",
    "ggtitle('log2(x+1)')\n",
    "\n",
    "p3 <- ggplot(data.frame(assay(rld)[, 1:2]), aes(x = SRR1039508, y = SRR1039509)) + \n",
    "geom_point() + \n",
    "geom_abline(intercept = 0, slope = 1, color = 'blue', linewidth = 2) +\n",
    "coord_fixed() +\n",
    "ggtitle('rlog')\n",
    "\n",
    "p4 <- ggplot(data.frame(assay(vsd)[, 1:2]), aes(x = SRR1039508, y = SRR1039509)) + \n",
    "geom_point() + \n",
    "geom_abline(intercept = 0, slope = 1, color = 'blue', linewidth = 2) +\n",
    "coord_fixed() +\n",
    "xlim(0,19)+\n",
    "ylim(0,19) +\n",
    "ggtitle('vst')\n",
    "\n",
    "options(repr.plot.width=15, repr.plot.height=5)\n",
    "\n",
    "p1 | p2 | p3 | p4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some scatter plots showing the raw counts and the log2 (left), rlog (middle), and VST (right) transformed counts. The x and y axes are expression values in two samples from the `airway` data. Each black dot is a gene and its position on the plot indicates its expression levels in two of the `airway` data samples. The blue line is indicating the 1:1 line -- points closer to this line have similar expression values in sample x and in sample y. Genes in the bottom left corner are lowly expressed and genes in the upper right corner are highly expressed.    \n",
    "\n",
    "We can see how genes with low counts (bottom left-hand corner) are highly variable on the ordinary logarithmic scale (`log2`), while the `rlog` transform and `vst` compress differences for the low count genes. Since these low-count values are likely artifacts introduced by sequencing error or random chance, their exclusion helps to prevent false positives.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the standard deviation of each row (gene) against the mean is also helpful for understanding how data transformations impact the mean:variance relationship:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p5 <- meanSdPlot(counts(dds), ranks = FALSE, plot = FALSE)\n",
    "p6 <- meanSdPlot(log2(counts(ddsESF, normalized=TRUE)[, 1:2] + 1), ranks = FALSE, plot = FALSE)\n",
    "p7 <- meanSdPlot(rlog(counts(dds), blind = FALSE), ranks = FALSE, plot = FALSE)\n",
    "p8 <- meanSdPlot(vst(counts(dds), blind = FALSE), ranks = FALSE, plot = FALSE)\n",
    "\n",
    "options(repr.plot.width=15, repr.plot.height=5)\n",
    "\n",
    "p5$gg + ggtitle('counts')| p6$gg + ggtitle('log2(x+1)') | p7$gg + ggtitle('rlog') | p8$gg + ggtitle('vst')  + scale_x_continuous(limits = c(0, 19))  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On each of these plots, the x-axis is the mean and the y-axis is the standard deviation. The color of each point indicates how many data points have a given mean and standard deviation. Ligher blue colors indicate a greater number of data points at a particular mean and standard deviation, while darker blue indicates fewer data points.\n",
    "\n",
    "If we look at the `counts` plot on the far left, we see outliers in the upper right of the plot, each colored dark blue to show that very few genes fall within each point - these observations would cause a great deal of bias in the results if uncorrected.\n",
    "\n",
    "Moving to the `log2` plot, we see that transforming our data using the logarithm with a small pseudocount transformation (i.e., the standard log2 transformation) amplifies the standard deviation when the\n",
    "values are close to 0. As a result, the low count genes with low signal-to-noise\n",
    "ratio will now overly contribute to sample-sample distances and PCA\n",
    "plots.  \n",
    "\n",
    "For genes with high counts, the rlog and VST transformations will give results similar\n",
    "to that of the ordinary log2 transformation of normalized counts. **For genes\n",
    "with lower counts, the values are shrunken towards the genes'\n",
    "averages across all samples**. The rlog-transformed or VST data then\n",
    "becomes approximately homoskedastic, and can be used directly for\n",
    "computing distances between samples, making PCA plots, or as input to\n",
    "downstream methods which perform best with homoskedastic data.\n",
    "\n",
    "**Which transformation to choose?** The **rlog tends to work well on\n",
    "small datasets (n < 30)**, sometimes outperforming the VST when there is\n",
    "a large range of sequencing depth across samples (an order of\n",
    "magnitude difference). However, the VST is much faster to compute and is less\n",
    "sensitive to high count outliers than the rlog. Therefore, **we recommend\n",
    "the VST for large datasets (hundreds of samples)**. You can perform both\n",
    "transformations and compare the `meanSdPlot` or PCA plots generated.\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\"><b>Note:</b> Using transformed or normalized counts is only for visualizing and exploring data -- statistical analysis packages like DESeq2 and edgeR should be used with raw counts matrices </div>\n",
    "\n",
    "In the above function calls, we specified `blind = FALSE`, which means\n",
    "that differences between cell lines and treatment (the variables in\n",
    "the design) will not contribute to the expected variance-mean trend of\n",
    "the experiment. The experimental design is not used directly in the\n",
    "transformation, only in estimating the global amount of variability in\n",
    "the counts.  For a fully *unsupervised* transformation, one can set\n",
    "`blind = TRUE` (which is the default)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring data at the sample level\n",
    "\n",
    "A useful step in an RNA-seq analysis is often to assess overall similarity between samples: Which samples are similar to each other, which are different? Does this fit to the expectation from the experiment's design?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we get started, tidy the counts and transformed counts matrices so they are easier to visualize:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log2_tidy <- gather(as.data.frame(log2(assay(ddsESF)+1)), \n",
    "        key = 'library', \n",
    "        value = 'counts', \n",
    "        SRR1039508:SRR1039521)\n",
    "\n",
    "rld_tidy <- gather(as.data.frame((assay(rld))), \n",
    "        key = 'library', \n",
    "        value = 'counts', \n",
    "        SRR1039508:SRR1039521)\n",
    "\n",
    "vsd_tidy <- gather(as.data.frame((assay(vsd))), \n",
    "        key = 'library', \n",
    "        value = 'counts', \n",
    "        SRR1039508:SRR1039521)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boxplots\n",
    "\n",
    "Boxplots of the count distributions in each sample are a good way to understand the effects these transformations have at the level of individual subjects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bp1 <- ggplot(log2_tidy, aes(x = library, y = counts)) + \n",
    "geom_boxplot() +\n",
    "ggtitle('log2(x+1)') +\n",
    "theme(axis.text.x = element_text(angle = 45, hjust = 1))\n",
    "\n",
    "bp2 <- ggplot(rld_tidy, aes(x = library, y = counts)) + \n",
    "geom_boxplot() +\n",
    "ggtitle('rlog')+\n",
    "theme(axis.text.x = element_text(angle = 45, hjust = 1))\n",
    "\n",
    "bp3 <- ggplot(vsd_tidy, aes(x = library, y = counts)) + \n",
    "geom_boxplot() +\n",
    "ggtitle('vsd')+\n",
    "theme(axis.text.x = element_text(angle = 45, hjust = 1))\n",
    "\n",
    "options(repr.plot.width=15, repr.plot.height=5)\n",
    "\n",
    "bp1 | bp2 | bp3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In R's boxplots, points classified as outliers are shown on the plots as circles. Neither the standard log2 or the rlog transformations include any outliers, while the vst transformations do have outliers at the upper end of the distributions.  These points are classified as outliers because the vst method compresses lower end points more substantially (which can be seen on the scatter plots above). The actual magnitude of these points does not significantly exceed the magnitude seen with the other transformations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Density plots\n",
    "\n",
    "Density plots are another good way to visualize if the counts distributions agree (or disagree) across the samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp1 <- ggplot(log2_tidy, aes(x = counts, color = library)) + \n",
    "geom_density() +\n",
    "ggtitle('log2(x+1)') +\n",
    "ylim(0,.8)\n",
    "\n",
    "dp2 <- ggplot(rld_tidy, aes(x = counts, color = library)) + \n",
    "geom_density() +\n",
    "ggtitle('rlog') +\n",
    "ylim(0,.8)\n",
    "\n",
    "dp3 <- ggplot(vsd_tidy, aes(x = counts, color = library)) + \n",
    "geom_density() +\n",
    "ggtitle('vst') +\n",
    "ylim(0,.8)\n",
    "\n",
    "dp1 | dp2 | dp3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the different axis scales again -- the vst transformation axes do not start at 0 because this method shifts all the low counts data toward their mean across samples. We can also see the previously discussed issues with how the log2 transformation handles low counts data. You might also notice that some of the rlog counts are below 0 -- again, this is expected as rlog borrows information across samples and implements a log transformation (which will result in negative values for counts below 1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample distances\n",
    "\n",
    "We can also use the R function dist to calculate the \"Euclidean distance\" between samples; the Euclidean distance is simply the higher-dimensional analogue of the human-measurable distance between points in 2 or 3 dimensions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To ensure we have a roughly equal contribution from all genes, we calculate it from the rlog-transformed data. We need to transpose the matrix of values using t, because the dist function expects the different samples to be rows of its argument, and different dimensions (here, genes) to be columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleDists <- dist(t(assay(rld)))\n",
    "sampleDists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can make a heatmap to visualize the correlations. \n",
    "\n",
    "In order to plot the sample distance matrix with the rows/columns arranged by the distances in our distance matrix, we manually provide sampleDists to the clustering_distance argument of the pheatmap function. Otherwise the pheatmap function would assume that the matrix contains the data values themselves, and would calculate distances between the rows/columns of the distance matrix, which is not desired. We also manually specify a blue color palette using the colorRampPalette function from the RColorBrewer package (this coloring step is optional)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the sampleDists into a matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleDistMatrix <- as.matrix(sampleDists)\n",
    "head(sampleDistMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Give the rows more meaningful names and clean up the column names since they will be redundant in the heatmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rownames(sampleDistMatrix) <- paste( rld$dex, rld$cell, sep = \" - \" )\n",
    "colnames(sampleDistMatrix) <- NULL\n",
    "head(sampleDistMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set our color palette for the heatmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors <- colorRampPalette( rev(brewer.pal(9, \"Blues\")) )(255)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we make the heatmap, specifying that we have already calculated the clustering distances and stored them in `sampleDists`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.height=3, repr.plot.width=5)\n",
    "pheatmap(sampleDistMatrix,\n",
    "         clustering_distance_rows = sampleDists,\n",
    "         clustering_distance_cols = sampleDists,\n",
    "         col = colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another option for calculating sample distances is to use the Poisson Distance, implemented in the PoiClaClu CRAN package. This measure of dissimilarity between counts also takes the inherent variance structure of counts into consideration when calculating the distances between samples. The PoissonDistance function takes the original count matrix (not normalized) with samples as rows instead of columns, so we need to transpose the counts in dds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poisd <- PoissonDistance(t(counts(ddsESF)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "convert to a matrix and fix the row and column names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samplePoisDistMatrix <- as.matrix( poisd$dd )\n",
    "rownames(samplePoisDistMatrix) <- paste( rld$dex, rld$cell, sep=\" - \" )\n",
    "colnames(samplePoisDistMatrix) <- NULL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then create the heatmap using the same color palette previously defined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.height=3, repr.plot.width=5)\n",
    "pheatmap(samplePoisDistMatrix,\n",
    "         clustering_distance_rows = poisd$dd,\n",
    "         clustering_distance_cols = poisd$dd,\n",
    "         col = colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In both cases, the distance matrices are indicating that there are differences between cell types, but they generally cluster based on their treatment before cell type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the `ggpairs` function to look at trends across samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=10, repr.plot.height=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggpairs(as.data.frame(log2(counts(ddsESF, normalized=TRUE) + 1)), title = 'log2(x+1)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggpairs(as.data.frame(rlog(counts(dds))), title = 'rlog')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggpairs(as.data.frame(vst(counts(dds))), title = 'vst')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting in the lower left corner of the ggpairs plot, we see several scatter plots that compare each of the samples. The scatter plots are just like the 1:1 plots shown earlier in this workshop, where each point is a gene and points that are closer to the 1:1 line show similar expression levels in the two samples compared. Again, we can see that the  log2 transformation has higher variance when the counts are closer to 0. The rlog transformation compresses the variation, while the vst transformation shifts counts away from 0 (note the axes that start at 4).   \n",
    "\n",
    "Distribution plots are shown on the diagonal -- we can see how each of these transformations influences the distribution of the counts -- note the axis ranges on the vst plots, which do not start at 0. The scatterplots are prone to over-plotting so it can be difficult to tell how many points are in each section of the plot -- the density plots help convey that information.    \n",
    "\n",
    "The upper right corner shows the Pearson correlations -- a concise way to make sure your samples are behaving the ways you anticipate (for example, that your replicates agree).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA\n",
    "\n",
    "Another way to visualize sample-to-sample distances is a\n",
    "principal components analysis (PCA). In this ordination method, the\n",
    "data points (here, the samples) are projected onto the 2D plane\n",
    "such that they spread out in the two directions that explain most of\n",
    "the differences (figure below). The x-axis is the direction (i.e., principal component) that separates the data\n",
    "points the most. The values of the samples in this direction are\n",
    "written *PC1* (\"principal component 1\"). The y-axis is a direction (another principal component, which must be *orthogonal* - perpendicular - to\n",
    "the first direction) that separates the data the second most. The\n",
    "values of the samples in this direction are written *PC2*.\n",
    "The percent of the total variance that is contained in the direction\n",
    "is printed in the axis label. Note that these percentages do not add to\n",
    "100%, because there are more dimensions - often many of them - that contain the remaining\n",
    "variance (although each of these remaining dimensions will by definition explain\n",
    "less than the two that we see)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colData(rld)\n",
    "options(repr.plot.width=7, repr.plot.height=7)\n",
    "plotPCA(rld, intgroup = c(\"dex\", \"cell\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have used the function `plotPCA()` that comes with *DESeq2* to plot the rlog transformed data.\n",
    "The two terms specified by `intgroup` are the interesting groups for\n",
    "labeling the samples; they tell the function to use them to choose\n",
    "colors. Each unique combination of treatment and cell line is given its own color.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\"><b>Exercise 3:</b> If you look at `?plotPCA`, you'll see that the default behavior of `plotPCA` is to plot the top 500 genes. Does the PCA change if you change those settings and set `ntop = 50`?  </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also build the PCA plot from scratch using the\n",
    "[ggplot2](https://cran.r-project.org/web/packages/ggplot2/index.html) package.\n",
    "This is done by asking the `plotPCA()` function\n",
    "to return the data used for plotting rather than building the plot.\n",
    "See the *ggplot2* [documentation](http://docs.ggplot2.org/current/)\n",
    "for more details on using `ggplot()`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run `plotPCA` with `returnData = TRUE`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcaData <- plotPCA(rld, intgroup = c( \"dex\", \"cell\"), returnData = TRUE)\n",
    "pcaData\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "percentVar <- round(100 * attr(pcaData, \"percentVar\"))\n",
    "percentVar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then use this structure to build up a second plot in a figure below, specifying that the\n",
    "color of the points should reflect dexamethasone treatment and the\n",
    "shape should reflect the cell line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot(pcaData, aes(x = PC1, y = PC2, color = dex, shape = cell)) +\n",
    "  geom_point(size =3) +\n",
    "  xlab(paste0(\"PC1: \", percentVar[1], \"% variance\")) +\n",
    "  ylab(paste0(\"PC2: \", percentVar[2], \"% variance\")) +\n",
    "  coord_fixed()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the PCA plot, we see that the differences between cells (the\n",
    "different plotting shapes) are considerable, though not stronger than the differences due to\n",
    "treatment with dexamethasone (red vs blue color). This shows why it will be important to\n",
    "account for this in differential testing by using a paired design\n",
    "(\"paired\", because each dex treated sample is paired with one\n",
    "untreated sample from the *same* cell line). We are already set up for\n",
    "this design by assigning the formula `~ cell + dex` earlier when specifying our model.\n",
    "\n",
    "**A note on outlier detection:** there are no hard, fast rules for identifying outliers; however, PCA can provide you with useful information for making decisions about which samples may be outliers. One general approach is to visually inspect the PCA bi-plot produced by the `plotPCA` function (i.e., the plot visualizing the samples along the two principal components explaining the most variance). Samples that are outliers will be many of orders of magnitude away from the main group of samples along principal component 1; the difference should be easy to identify visually and usually will not require statistical justification. A second approach is to quantify outliers using z-score standardization. Specifically, you take the principal component 1 values from your PCA object and convert them to z-scores using the following standardization: \n",
    "\n",
    "$x_{gi}$ $\\leftarrow$ $\\frac{(x_{gi} - \\bar{x}_g)}{s_g}$  \n",
    "\n",
    "In our context, this just means we take each value on the first principal component, subtract it from the mean of all the values on the first principal component, and then divide this difference by the standard deviation of all the values on the first principal component. To do so, we would access the first principal component from our PCA object using `pcaData$PC1` and standardize these scores using `scale(pcaData$PC1)`. Once you have the standardized scores, outliers can be identified as z-scores greater than 3 or 6 in absolute value, depending on how conservative one wishes to be.\n",
    "\n",
    "**Note:** It is important to note that the two aforementioned approaches are just general guidelines and identification of outliers requires statistical knowledge as well as detailed knowledge of experimental details, lab techniques used, etc. It is advised that you consult with a statistician when making such decisions. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MDS\n",
    "Another plot, very similar to the PCA plot, can be made using the \n",
    "*multidimensional scaling* (MDS) function in base R. MDS is a useful technique for examining sample-to-sample distances that is **used when we don't have the original data, but instead have only a matrix of distances**. Here we compute the MDS for the distances calculated from the *rlog*\n",
    "transformed counts and plot these in a figure below.\n",
    "\n",
    "This method, unlike PCA, does not show the percent of variance explained by the individual components."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**MDS plot using rlog-transformed values.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleDists <- dist(t(assay(rld)))\n",
    "sampleDistMatrix <- as.matrix( sampleDists )\n",
    "mds <- cbind(as.data.frame(colData(rld)), cmdscale(sampleDistMatrix))\n",
    "ggplot(mds, aes(x = `1`, y = `2`, color = dex, shape = cell)) +\n",
    "  geom_point(size = 3) + coord_fixed()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**MDS plot using the Poisson Distance.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(PoiClaClu)\n",
    "poisd <- PoissonDistance(t(counts(ddsESF)))\n",
    "samplePoisDistMatrix <- as.matrix( poisd$dd )\n",
    "rownames(samplePoisDistMatrix) <- paste( rld$dex, rld$cell, sep=\" - \" )\n",
    "\n",
    "mdsPois <- cbind(as.data.frame(colData(ddsESF)),\n",
    "   cmdscale(samplePoisDistMatrix))\n",
    "ggplot(mdsPois, aes(x = `1`, y = `2`, color = dex, shape = cell)) +\n",
    "  geom_point(size = 3) + coord_fixed()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A key difference between PCA and MDS is that PCA preserves the covariance of the data, while MDS preserves the distance between the data points. Depending on what type of distance you use for your MDS analysis, these two results can be the same (if Euclidean distance is used) or different.  Here, we have used the Poisson distance (since it takes the structure of the variance into account), so the MDS and PCA results are not exactly the same.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Differential expression analysis\n",
    "\n",
    "Now that we have completed a thorough examination of the data and are happy with the quality, we are ready to proceed with the differential expression analysis.\n",
    "## Running the differential expression pipeline\n",
    "\n",
    "As a reminder the initial input for the analysis was the count matrix and experiment design (i.e., model formula) into the DESeqDataSet() function `dds <- DESeqDataSet(airway, design = ~ cell + dex)` and the factor levels must also be set `dds$dex <- relevel(dds$dex, ref = 'untrt')`.\n",
    "We can now run the differential expression pipeline on the raw counts with a single call to the function *DESeq*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dds <- DESeq(dds) # run differential expression analysis pipeline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** This function will print out a message for the various steps it\n",
    "performs. These steps are: \n",
    "\n",
    "1.) the estimation of size factors (normalization); \n",
    "\n",
    "2.) the estimation of dispersion values for each gene; \n",
    "\n",
    "3.) and fitting a generalized linear model (GLM) with Wald tests.\n",
    "\n",
    "A *DESeqDataSet* is returned that contains all the fitted\n",
    "parameters within it, and the following section describes how to\n",
    "extract out results tables of interest from this object.\n",
    "\n",
    "The steps performed the DESeq() function are described in more detail below and in the manual page for\n",
    "*DESeq*, which can be accessed by typing `?DESeq`.\n",
    "\n",
    "### Normalization: A Note on Counts in DESeq2 \n",
    "\n",
    "The starting point of the DESeq2 analysis is a count matrix. Given the count matrix, the first step in a differential expression analysis workflow is to perform count normalization so that we can make accurate comparisons of gene expression between samples. Normalization is the process of scaling raw count values to account for the factors we are not interested in. These are factors such as sequencing depth, gene length, and RNA composition. By accounting for such factors, we ensure that expression levels are more comparable between and/or within samples. \n",
    "\n",
    "It is important to note that DESeq2 automatically handles normalization internally when performing its analysis. It does so in an attempt to avoid the effects of count magnitude on differential expression status - otherwise, nonsignificant genes with very large counts would appear as false positives and significant genes with very small counts would appear as false negatives.\n",
    "\n",
    "Common methods of normalization include Trimmed Mean of M-Values (TMM), Relative Log Expression (RLE), and Median Ratio Normalization (MRN). TMM, used by the edgeR package, and RLE, yield very similar results. DESeq2 uses the RLE approach, which accounts for sequencing depth and RNA composition. Given that tools for differential expression analysis are comparing the counts between sample groups for the same gene, gene length is a factor that does not need to be accounted for with this design.\n",
    "\n",
    "**A word of caution.** When obtaining a result from the DESeq() function highlighting a list of genes, people will often look back at their original count table using a program like Excel and examine the counts tabulated therein. Doing so is inadvisable and misleading, as these non-normalized counts are not proper indicators of significance.\n",
    "\n",
    "Here, we will show how different normalized and non-normalized counts can be, illustrating the importance of normalization. Note that these plots are on the log scale, so a linear distance represents an exponential change. Very low counts are unaffected by normalization, as the normalization process adds a 0.5 pseudocount, but higher counts change substantially. Here, one of the two lower counts in the treatment group are almost halved and almost doubled in magnitude, respectively - an enormous change that would be misleading if viewed without factoring in the normalization.\n",
    "\n",
    "#### Counts plot\n",
    "A quick way to visualize the counts for a particular gene is to use the plotCounts function that takes as arguments the DESeqDataSet, a gene name, and the group over which to plot the counts (figure below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotCounts(dds, gene = \"ENSG00000132518\", intgroup=c(\"dex\"), main = \"Normalized\")\n",
    "plotCounts(dds, gene = \"ENSG00000132518\", intgroup=c(\"dex\"), normalized = F, main = \"Non-normalized\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dispersion Estimation\n",
    "\n",
    "The next step in differential expression analysis is the estimation of gene-wise dispersions. **What is dispersion and what does it represent?** Simply put, dispersion is a measure of spread or variability in the data. The measure of dispersion used with DESeq2 is related to both the mean and variance of the data. Specifically, The **DESeq2 dispersion estimates are inversely related to the mean and directly related to variance**. Based on this relationship, the dispersion is higher for small mean counts and lower for large mean counts. The dispersion estimates for genes with the same mean will differ only based on their variance. Therefore, the dispersion estimates reflect the variance in gene expression for a given mean value. \n",
    "\n",
    "The plot of mean versus variance in count data below shows the variance in gene expression increases with the mean expression (each black dot is a gene). Notice that the relationship between mean and variance is linear on the log scale, and for higher means, we could predict the variance relatively accurately given the mean. However, **for low mean counts, the variance estimates have a much larger spread; therefore, the dispersion estimates will differ much more between genes with small means**.\n",
    "\n",
    "<img src=\"Images/DispersionEstimation1.png\">\n",
    "\n",
    "**Why does dispersion matter for our modelling process?** To accurately model sequencing counts, we need to generate accurate estimates of within-group variation (variation between replicates of the same sample group) for each gene. With only a few (3-6) replicates per group, the estimates of variation for each gene are often unreliable (due to the large differences in dispersion for genes with similar means).\n",
    "\n",
    "To address this problem, DESeq2 shares information across genes to generate more accurate estimates of variation based on the mean expression level of the gene using a method called ‘shrinkage’. DESeq2 assumes that genes with similar expression levels have similar dispersion.\n",
    "\n",
    "**Estimating the dispersion for each gene separately:** To model the dispersion based on expression level (mean counts of replicates), the dispersion for each gene is estimated using maximum likelihood estimation. Maximum likelihood is an approach for estimating a parameter(s) by maximizing what's known as a likelihood function. The likelihood function is simply a function representing the support the data provides for each value of the parameter (e.g., dispersion). Thus, given the count values of the replicates, the most likely estimate of dispersion is calculated.\n",
    "\n",
    "**Fit curve to gene-wise dispersion estimates:** The next step in the workflow is to fit a curve to the dispersion estimates for each gene. The idea behind fitting a curve to the data is that different genes will have different scales of biological variability, but, over all genes, there will be a distribution of reasonable estimates of dispersion.\n",
    "\n",
    "This curve is displayed as a red line in the figure below, which plots the estimate for the expected dispersion value for genes of a given expression strength. Each black dot is a gene with an associated mean expression level and maximum likelihood estimation (MLE) of the dispersion (Step 1).\n",
    "\n",
    "<img src=\"Images/DispersionEstimation2.png\">\n",
    "\n",
    "\n",
    "**Shrink gene-wise dispersion estimates toward the values predicted by the curve**\n",
    "The next step in the workflow is to shrink the gene-wise dispersion estimates toward the expected dispersion values.\n",
    "\n",
    "The curve allows for more accurate identification of differentially expressed genes when sample sizes are small, and the strength of the shrinkage for each gene depends on :\n",
    "\n",
    "how close gene dispersions are from the curve\n",
    "\n",
    "sample size (more samples = less shrinkage)\n",
    "\n",
    "This shrinkage method is particularly important to reduce false positives in the differential expression analysis. Genes with low dispersion estimates are shrunken towards the curve, and the more accurate, higher shrunken values are output for fitting of the model and differential expression testing.\n",
    "\n",
    "Dispersion estimates that are slightly above the curve are also shrunk toward the curve for better dispersion estimation; however, genes with extremely high dispersion values are not. This is due to the likelihood that the gene does not follow the modeling assumptions and has higher variability than others for biological or technical reasons. Shrinking the values toward the curve could result in false positives, so these values are not shrunken. These genes are shown surrounded by blue circles below.\n",
    "\n",
    "<img src=\"Images/DispersionEstimation3.png\">\n",
    "\n",
    "This is a good plot to examine to ensure your data is a good fit for the DESeq2 model. You expect your data to generally scatter around the curve, with the dispersion decreasing with increasing mean expression levels. If you see a cloud or different shapes, then you might want to explore your data more to see if you have contamination (mitochondrial, etc.) or outlier samples. Note how much shrinkage you get across the whole range of means in the plotDispEsts() plot for any experiment with low degrees of freedom.\n",
    "\n",
    "Examples of worrisome dispersion plots are shown below:\n",
    "\n",
    "<img src=\"Images/DispersionEstimation4.png\">\n",
    "\n",
    "<img src=\"Images/DispersionEstimation5.png\">\n",
    "\n",
    "\n",
    "### Hypothesis Testing: Wald test\n",
    "\n",
    "The count data we work with for differential expression is assumed to follow a negative binomial distribution (Poisson-like), and thus the DESeq2 authors model the data as such. Specifically, after estimating dispersion, DESeq2 fits a negative binomial regression model to each gene to make inferences about differential expression. To make these inferences, we perform hypothesis testing on the regression coefficients in the negative binomial regression model using the Wald test. \n",
    "\n",
    "**Hypothesis testing**\n",
    "\n",
    "The first step in hypothesis testing is to set up a null hypothesis for each gene. In our case, the null hypothesis is that there is no differential expression across sample groups (i.e., regression coefficient == 0). Notice that we can do this without observing any data, because it is based on a thought experiment. Subsequently, we use a statistical test to determine if, based on the observed data, we should accept the null hypothesis as being true. \n",
    "\n",
    "**Wald test**\n",
    "\n",
    "With DESeq2, the Wald test is the default used for hypothesis testing when comparing two groups. The Wald test is a test of hypothesis usually performed on parameters that have been estimated via maximum likelihood. Maximum likelihood is just an estimation method for deriving the parameters to a statistical model. In our case, we are testing each gene model coefficient (i.e., regression coefficient) to see if it is significantly different from zero.  \n",
    "\n",
    "**DESeq2 implements the Wald test by:**\n",
    "\n",
    "1.) Regression coefficients of the model are estimated for each sample group along with their standard error. Note that in DE the coefficents are the estimates for the log2 foldchanges for each sample group.\n",
    "\n",
    "2.) Take the regression coefficient and divide it by its standard error, which results in a z-statistic\n",
    "\n",
    "3.) The z-statistic is then compared to a standard normal distribution, and from this a p-value is derived. which reports the probability of observing a z-statistic at least as extreme as the one we observed, given that the null hypothesis is true\n",
    "\n",
    "4.) If the p-value is small, we reject the null hypothesis and state that there is evidence against the null (i.e. the gene is differentially expressed) \n",
    "\n",
    "<div class=\"alert alert-block alert-warning\"><b>Note:</b> The Wald test and the previously discussed Likelihood Ratio Test (LRT) are essentially the same thing and test the same hypothesis. The big difference is their finite sample performance. Specifically, in smaller to moderate sample sizes, the LRT is to be preferred, as the Wald test becomes less reliable. </div>\n",
    "\n",
    "## Building the results table\n",
    "\n",
    "Calling `results` without any arguments will extract the estimated log2 fold changes and p-values **for the last variable in the design **formula**. If there are more than 2 levels for this variable, `results` will extract the results table for a comparison of the last level over the first level.\n",
    "\n",
    "Here, our model was `~ cell + dex`, so the default comparison printed at the top of the output is `dex trt vs dex untrt`.  As the `airway` dataset's `dex` variable only has two levels, there is no need to worry about explicitly specifying level comparisons in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a results table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res <- results(dds)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Note on Statistical Terminology \n",
    "\n",
    "* For large samples, violations of distributional assumptions have diminishing importance\n",
    "* A sampling distribution is a probability distribution of a statistic (e.g., the mean) obtained from repeated random sampling\n",
    "* _Standard Error_ is the standard deviation of a normally distributed sampling distribution\n",
    "* The null hypothesis is an assumption about the form of a sampling distribution\n",
    "\n",
    "We have thus far defined the p-value with respect to the null hypothesis, but have yet to state what the null hypothesis actually means.  Generally, in a biological experiment, the null hypothesis is that there is no difference between two (or more) different groups according so some metric being studied e.g. gene expression.  The alternative hypothesis is that there is some difference.  To determine which of these hypotheses is more likely to be correct, we calculate a p-value.\n",
    "\n",
    "**p-value:** the probability of getting a result at least as extreme as the one observed *given* that the null hypothesis is true.\n",
    "\n",
    "Thus, a low p-value means it's unlikely that you'd see what you're seeing if the null hypothesis is true, providing evidence that the null hypothesis is not true - in classical statistics, if that probability is lower than .05 (5%), we *reject the null hypothesis*.  Note that this is not strictly identical to as *accepting the alternative hypothesis*, but with the null rejected, we only have the alternative hypothesis left (under this framework)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Note On Multiple testing\n",
    "\n",
    "In high-throughput biology, we are careful to not use the p-values\n",
    "directly as evidence against the null, but instead to correct for\n",
    "*multiple testing*. The reason we do this is because we are running multiple (sometimes thousands of) hypothesis tests when performing differential expression analysis, and this increases the chances of observing false positives (Type I error). To deal with this, *DESeq2* uses the [Benjamini-Hochberg (BH)](https://www.jstor.org/stable) adjustment as implemented in\n",
    "the base R `p.adjust` function; in brief, this method calculates, for\n",
    "each gene, an adjusted p-value that answers the following question:\n",
    "if one called significant all genes with an adjusted p-value less than or\n",
    "equal to this gene's adjusted p-value threshold, what would be the fraction\n",
    "of false positives (the *false discovery rate*, FDR) among them, in\n",
    "the sense of the calculation outlined above? These values, called the\n",
    "BH-adjusted p-values, are given in the column `padj` of the `res`\n",
    "object.\n",
    "\n",
    "\n",
    "The BH correction is appropriate in the sense that it is not as overlyconservative as other correction approaches and maintains the desired coverage rate; however, it is important to note that this aproach isn't robust to violations of indepednence. This can be problematic in our context because the multiple hypothesis tests performed in differential expression analysis may very well be dependent. To deal with this, we can use a Bonferroni correction appraoch, which simply involves dividing the conventional p-value (usually 0.05) by the number of hypothesis tests being performed. This then creates a new nominal p-value against which the observed p-values are compared for signficiance. This can be done in DESeq2 by using the `p.adjust=\"Bonferroni\"` argument. This approach is more robust to violations of independence than the BH correction; however, it is also important to note that this approach is known to be overlyconservative and thus it is harder to achieve signficance with. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## Other comparisons\n",
    "\n",
    "The results for a comparison of any two levels of a variable can be extracted using the `contrast` argument to *results*. The user should specify three values: \n",
    "- the name of the variable   \n",
    "- the name of the level for the numerator    \n",
    "- and the name of the level for the denominator       \n",
    "\n",
    "Here we extract results for the log2 of thefold change of one cell line over another. We will explicitly specify which cell lines to compare since there are more than two levels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results(dds, contrast = c(\"cell\", \"N061011\", \"N61311\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are additional ways to build results tables for certain comparisons after running *DESeq* once. If results for an interaction term are desired, the `name` argument of `results` should be used. Please see the help page for the `results` function for details on the additional ways to build results tables. In particular, the **Examples** section of the help page for `results` gives some pertinent examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We subset the results table to these genes and then sort it by the\n",
    "log2 fold change estimate to get the significant genes with the\n",
    "strongest downregulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resSig <- subset(res, padj < 0.1)\n",
    "head(resSig[ order(resSig$log2FoldChange), ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and with the strongest upregulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(resSig[ order(resSig$log2FoldChange, decreasing = TRUE), ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diagnostic plots\n",
    "\n",
    "## MA-plot\n",
    "\n",
    "An *MA-plot* ([Dudoit et al., Statistica Sinica, 2002](http://wwwf.imperial.ac.uk/~das01/BioinformaticsMSc/Papers/sinica.final.pdf)) provides a useful overview for\n",
    "the distribution of the estimated coefficients in the model,\n",
    "i.e. the comparisons of interest, across all genes.\n",
    "On the y-axis, the \"M\" stands for \"minus\" --\n",
    "subtraction of log values is equivalent to the log of the ratio -- and\n",
    "on the x-axis, the \"A\" stands for \"average\". You may hear this plot\n",
    "also referred to as a \"mean-difference plot\", or a \"Bland-Altman plot\".\n",
    "\n",
    "Before making the MA-plot, we use the\n",
    "`lfcShrink` function to shrink the log2 fold-changes for the\n",
    "comparison of dex treated vs untreated samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res <- results(dds)\n",
    "res <- DESeq2::lfcShrink(dds, coef = \"dex_trt_vs_untrt\", res = results(dds))\n",
    "DESeq2::plotMA(res, ylim=c(-5,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**An MA-plot of changes induced by treatment.**\n",
    "The log2 fold-change for a particular comparison is plotted on the y-axis and the average of the counts normalized by size factor is shown on the x-axis.\n",
    "Each gene is represented with a dot. Genes with an adjusted *p* value below a threshold (here 0.1, the default) are shown in blue.\n",
    "\n",
    "The *DESeq2* package uses a Bayesian procedure to moderate (or\n",
    "\"shrink\") log2 fold changes from genes with very low counts and highly\n",
    "variable counts, as can be seen by the narrowing of the vertical\n",
    "spread of points on the left side of the MA-plot. As shown above, the\n",
    "`lfcShrink` function performs this operation.  For a detailed\n",
    "explanation of the rationale of moderated fold changes, please see the\n",
    "*DESeq2* paper ([Love and Huber, Genome Biology, 2014](http://wwwf.imperial.ac.uk/~das01/BioinformaticsMSc/Papers/sinica.final.pdf)).\n",
    "\n",
    "If we had not used statistical moderation to shrink the noisy log2\n",
    "fold changes, we would have instead seen the following plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.noshr <- results(dds)\n",
    "plotMA(res.noshr, ylim = c(-5, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can label individual points on the MA-plot as well. Here we use the `with` R function to plot a circle and text for a selected row of the results object. Within the `with` function, only the `baseMean` and `log2FoldChange` values for the selected rows of `res` are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMA(res, ylim = c(-5,5))\n",
    "topGene <- rownames(res)[which.min(res$padj)]\n",
    "with(res[topGene, ], {\n",
    "  points(baseMean, log2FoldChange, col=\"dodgerblue\", cex=2, lwd=2)\n",
    "  text(baseMean, log2FoldChange, topGene, pos=2, col=\"dodgerblue\")\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Another useful diagnostic plot is the histogram of the p-values (figure below). This plot is best formed by excluding genes with very small counts, which otherwise generate spikes in the histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist(res$pvalue[res$baseMean > 1], breaks = 0:20/20,\n",
    "     col = \"grey50\", border = \"white\", main = \"p-value Histogram\", xlab = \"p-value\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Histogram of *p* values for genes with mean normalized count larger than 1.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Volcano plot\n",
    "\n",
    "A typical volcano plot is a scatterplot of $-log_{10}$(p-value) vs. $log_2$(fold-change). It allows you to visualize the difference in expression between groups against the statistical significance of that difference. Some things to look for:\n",
    "* are diffentially expressed genes skewed towards up or down-regulation?\n",
    "* are there a lot of significant but very low fold-change genes? This can result from a confounded batch effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "volcanoPlot <- function(res, lfc=2, pval=0.01){\n",
    "    tab = data.frame(logFC = res$log2FoldChange, negLogPval = -log10(res$pvalue))\n",
    "    plot(tab, pch = 16, cex = 0.6, xlab = expression(log[2]~fold~change), ylab = expression(-log[10]~pvalue))\n",
    "    signGenes = (abs(tab$logFC) > lfc & tab$negLogPval > -log10(pval))\n",
    "    points(tab[signGenes, ], pch = 16, cex = 0.8, col = \"red\") \n",
    "    abline(h = -log10(pval), col = \"green3\", lty = 2) \n",
    "    abline(v = c(-lfc, lfc), col = \"blue\", lty = 2) \n",
    "    mtext(paste(\"pval =\", pval), side = 4, at = -log10(pval), cex = 0.8, line = 0.5, las = 1) \n",
    "    mtext(c(paste(\"-\", lfc, \"fold\"), paste(\"+\", lfc, \"fold\")), side = 3, at = c(-lfc, lfc), cex = 0.8, line = 0.5)\n",
    "}\n",
    "volcanoPlot(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Building and Analytic Approaches \n",
    "\n",
    "This section covers some model building basics and elementary experimental design for bulk sequence data. \n",
    "\n",
    "## Model formulae and terminology\n",
    "\n",
    "To begin, we must familiarize ourselves with the \"model formula\" interface in R. For important details complementary to this workshop, please explore the following link: \n",
    "\n",
    "[Model formulae tutorial](http://ww2.coastal.edu/kingw/statistics/R-tutorials/formulae.html)\n",
    "\n",
    "* It is important to know the model formula interface because many regression functions in R use it.\n",
    "* The formula determines the model that will be built (and tested) by the R procedure. The basic format is:\n",
    "\n",
    "`> response variable ~ explanatory variables`\n",
    "\n",
    "* The tilde means \"is modeled by\" or \"is modeled as a function of.\" \n",
    "\n",
    "### Terminology\n",
    "\n",
    "<img src=\"Images/Terminology.png\">\n",
    "\n",
    "Note: order generally doesn't matter (u+v OR v+u)\n",
    "\n",
    "*Examples*: `~treatment`, `~treatment + time`, `~treatment * time`\n",
    "\n",
    "\n",
    "## The Design Matrix\n",
    "\n",
    "To build models, we must decide upon a design matrix. To illustrate this matrix, we will use an example of a multiple linear regression model:\n",
    "\n",
    "$y_i = \\beta_0 + \\beta_1 x_{1i} + \\beta_2 x_{2i} + ... + \\beta_p x_{pi} + \\epsilon_i$\n",
    "\n",
    "- $x_{ji}$ is the value of predictor $x_j$ for observation $i$\n",
    "\n",
    "\n",
    "Matrix notation for the multiple linear regression model:\n",
    "\n",
    "$$\n",
    "\\,\n",
    "\\begin{pmatrix}\n",
    "Y_1\\\\\n",
    "Y_2\\\\\n",
    "\\vdots\\\\\n",
    "Y_N\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "1&x_1\\\\\n",
    "1&x_2\\\\\n",
    "\\vdots\\\\\n",
    "1&x_N\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\beta_0\\\\\n",
    "\\beta_1\n",
    "\\end{pmatrix} +\n",
    "\\begin{pmatrix}\n",
    "\\varepsilon_1\\\\\n",
    "\\varepsilon_2\\\\\n",
    "\\vdots\\\\\n",
    "\\varepsilon_N\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "where: \n",
    "* $\\mathbf{Y}$ is the observed expression for all read counts (from all conditions) for a given gene \n",
    "* $\\mathbf{X}$ encodes the condition (for RNA-seq, this is very often a discrete factor, e.g., “WT” or “mutant”) \n",
    "* $\\beta_0$ is the intercept, which represents the expression values of the baseline condition\n",
    "* $\\beta_1$ is the regression coefficient, which captures the difference between Y from samples of different conditions \n",
    "* $\\boldsymbol{\\varepsilon}$ captures the error or uncertainty, i.e. the difference of the regression estimates from the observed expression values. \n",
    "\n",
    "or more simply: \n",
    "\n",
    "$$\n",
    "\\mathbf{Y}=\\mathbf{X}\\boldsymbol{\\beta}+\\boldsymbol{\\varepsilon}\n",
    "$$\n",
    "\n",
   "* The design matrix is $\\mathbf{X}$\n",
    "    * which the computer will take as a given when solving for $\\boldsymbol{\\beta}$ by minimizing the sum of squares of residuals $\\boldsymbol{\\varepsilon}$.\n",
    "    \n",
    "\n",
    "It's important to note that there are multiple possible and reasonable design matrices for a given study design\n",
    "* the model formula is used to encode a model matrix, for example, we can use the model formula in the code below to set up a design matrix with one factor variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group <- factor( c(1, 1, 2, 2) )\n",
    "model.matrix(~ group)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we forgot to code group as a factor? The result is that R will enter the variable into the matrix as a single numeric variable, with different numbers representing group membership, as shown in the example below. This is different than the result one gets when using the factor() function, which gives k-1 dummy variables for indicating group membership, where k represents the number of groups. When performing a regression, this can make a signficiant difference for your results, so it is important to be mindful of this coding!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, the factor variable can have more than two levels or groupings. For example, if we have a single factor predictor variable representing three groups, the design matrix would be set up with k-1 (3-1) dummy variables as follows:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group <- factor(c(1,1,2,2,3,3)) # factor variable with three groupings \n",
    "model.matrix(~ group)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contrasts             \n",
    "It's important to note that using dummy variables in our design matrix involves leaving one group out for comparative purposes. The group left out of the design matrix is known as the baseline group (or the reference level) and is the one that other groups are contrasted against. By default, R will choose the baseline group based on either numerical or alphabetical order. However, you may be interested in specific comparisons, and so will have a need to change the baseline group. You can change the baseline group manually as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group <- factor(c(1,1,2,2,3,3)) # factor variable with three groups \n",
    "group <- relevel(x=group, ref=3) # use the ref argument to make group 3 our baseline group  \n",
    "model.matrix(~ group)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Experimental Design        \n",
    "\n",
    "For further information, we encourage you to explore [this vignette to approriately set up a design matrix to suit your experimental designs](https://bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/designmatrices.html)  \n",
    "\n",
    "## Single factor design\n",
    "\n",
    "In single factor designs, you have a single factor predictor variable of interest, such as treatment condition. So, for example, you may have three mice who receive a treatment and three mice who do not and therefore serve as a control. This would look as follows: \n",
    "\n",
    "<img src=\"Images/Single_factor_mice.png\">\n",
    "\n",
    "Under this design, you would be interested in the average group difference on the response variable between the two condition types. To examine something like this, we could code our design matrix as follows: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Single factor design \n",
    "condition <- factor( c(\"treatment\", \"treatment\", \"treatment\", \"control\", \"control\", \"control\") )\n",
    "model.matrix(~ condition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-factor design: Two factors \n",
    "\n",
    "The model matrix can be extended to represent multiple factor predictor variables, i.e. multiple (linear/log-linear/logistic) regression. This is useful, for example, in experiments where there are more than one factor influencing the outcome. In our previous mice example, let's say that in addition to condition, we introduced a variable capturing mutation. We would then want to account for this factor by adding it to the design and thus control for the additional variation in the outcome due to mice mutation/color. This illustration illustrates such a design: \n",
    "\n",
    "<img src=\"Images/Two_factor_mice.png\">\n",
    "\n",
    "Another example is if you have experimental condition samples (i.e., treated vs untreated) that are balanced across experimental batches, you can include batch as a second factor to the design. A batch variable might be something like the day that you extracted the RNA from your sample or the cages the mice were kept in during the course of the experiment. In either case, the addition of factors to your design serves to control for additional variation in the outcome (i.e., counts) and increases the sensitivity for finding differences due to experimental condition. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Blocking \n",
    "\n",
    "Releated to multi-factor designs is the concept of blocking. Blocking is the arranging of experimental units into groups that are similar to one another. The differences between the \"blocks\" of individuals is \"nuisance\" variability that we want to control. In other words, we use this strategy to partial out (essentially remove) variability due to irrelavant factors in order to improve the accuracy and sensitivity of our estimates for variables that are of primary interest. In biology experiments, we typically use blocking to apply an experimental randomization to more homogeneous \"blocks\" of individuals.         \n",
    "\n",
    "For example, we could use blocking to assign treated and control specimens, or case & control specimens, equally and randomly to batches. This would be referred to as **blocking on batch** and would be included as part of your experimental design. However, if such a consideration is not taken during the experimental design phase of your study, you can still control for nuisance variables, e.g., batch, afterwards by including the block as a factor in your regression analysis, as discussed and illustrated in the multi-factor design example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of controlling for batch in regression analysis \n",
    "treatment <- factor( c( \"treatment\", \"treatment\", \"treatment\", \"control\", \"control\", \"control\") )\n",
    "batch <- factor (c (\"Cage_1\", \"Cage_1\", \"Cage_2\", \"Cage_1\", \"Cage_2\", \"Cage_2\") )\n",
    "model.matrix(~ condition + batch) #here we are blocking on batch, which is cage mice were kept in "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above example includes batch as a factor in our regression analysis, thereby controlling nuisance variance. We now can get a better estimate of the effect associated with experimental condition. While this is an acceptable approach, note that it is always advised that you include blocking into your experimental design whenever possible. Incorporating a randomized block design is preferable to using blocks as controls in a regression because it often simplifies statistical analyses and allows for stronger, more credible inferences. \n",
    "\n",
    "## Interactions\n",
    "\n",
    "When multiple factors are being considered, you are able to specify interaction terms via the : operator. Generally speaking, interaction terms tell us if the effect of one predictor on the response variable varies across the levels (or different values) of another predictor variable (see bar chart below for visual representation of this). \n",
    "\n",
    "<img src=\"Images/Interaction.png\">\n",
    "\n",
    "We can specify a regression with two predictor variables and an interaction term as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code for examining an interaction effect, assuming 'low dose' is control and 'high dose' is treatment\n",
    "condition <- factor( c( \"treatment\", \"treatment\", \"treatment\", \"control\", \"control\", \"control\") )\n",
    "sex <- factor (c (\"female\", \"female\", \"male\", \"female\", \"male\", \"male\") )\n",
    "model.matrix(~ condition + sex + condition:sex) #here we use the : operator to include interaction effect "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repeated Measures/Time Series\n",
    "\n",
    "It is important to briefly speak about data that includes multiple time points, as analyzing data arising from this design is not as straightforward. There are a number of ways to analyze time-series experiments, depending on the biological question of interest. In order to test for any differences over multiple time points, one can use a design including the time factor, and then test for time effects using the likelihood ratio test as described in the following section, where the time factor is removed in the reduced formula. For a control and treatment time series, one can use a design formula containing the condition factor, the time factor, and the interaction of the two (as discussed above). In this case, using the likelihood ratio test with a reduced model which does not contain the interaction term between condition and time will test whether the condition induces a change in gene expression at any time point after the reference level time point (time 0). An example of the later analysis is provided in our RNA-seq workflow.        \n",
    "\n",
    "## Likelihood ratio test\n",
    "\n",
    "DESeq2 offers the likelihood ratio test (LRT), which is a statistical test that assesses the goodness of fit between two competing models. This test is particularly useful for time-series experiments, as mentioned above. The hypothesis under this test examines two models for the counts, a full model with a certain number of terms and a reduced model, in which some of the terms of the full model are removed. The test determines if the increased likelihood of the data using the extra terms in the full model is more than expected if those extra terms are truly zero.\n",
    "\n",
    "The LRT is therefore useful for testing multiple terms at once, for example, testing if there is a significant time by condition effect in a time series experiment. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.\n",
    "\n",
    "The likelihood ratio test can be performed by specifying test=\"LRT\" when using the DESeq function, and providing a reduced design formula, e.g. one in which a number of terms from design(dds) are removed (e.g., time x condition). The degrees of freedom for the test is obtained from the difference between the number of parameters in the two models. A likelihood ratio test for the control and treatment time-series experiment example, where the full design was ~condition + time + condition:time, would look like:\n",
    "\n",
    "`dds <- DESeq(dds, test=\"LRT\", reduced=~condition + time)` \n",
    "\n",
    "`res <- results(dds)` \n",
    "\n",
    "**Important Note:** The above is just one example of how to deal with time series data. In general, data involving measurements taken over time requires careful consideration, as repeated measurements create a statistical dependency in your data that needs to be accounted for. If this dependency structure goes unaccounted for or is inappropriately handled, your estimates, and thus the inferences you make, can become biased. Therefore, it is recommended that you consult with a statistician when analyzing repeated measure data.\n"
   ]
  },

  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R Packages for performing differential expression analyses \n",
    "  \n",
    "We could use a variety of Bioconductor packages for more exploration and differential expression analysis of the count data, including:\n",
    "* [DESeq2](http://www.bioconductor.org/packages/DESeq2)\n",
    "* [edgeR](http://www.bioconductor.org/packages/edgeR)\n",
    "* [limma](http://www.bioconductor.org/packages/limma) (use Voom method for RNA-seq)\n",
    "\n",
    "[Schurch *et al.* 2016](https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/27022035/) compared performance of different statistical methods\n",
    "for RNA-seq using a large number of biological replicates and can help\n",
    "users to decide which tools make sense to use, and how many\n",
    "biological replicates are necessary to obtain a certain sensitivity. \n",
    "\n",
    "###  Comparison of edgeR / DESeq2 / limma-voom\n",
    "\n",
    "**DESeq2**, **edgeR**, and **limma-voom** are all well-documented packages. They each implement similar key features:        \n",
    "\n",
    "* Empirical Bayes \"regularization\" or \"shrinkage\" of expression variance\n",
    "   + borrows prior expectation for variance from all genes, reduces false positives in small samples\n",
    "* Use of model formula, model matrices, and contrasts for flexible differential expression analysis      \n",
    "* Reporting of log fold-change and False Discovery Rate \n",
    "\n",
    "Some ways they differ are:\n",
    "* Error terms used: **DESeq2** and **edgeR** use log-linear Generalized Linear Model with negative binomial error term , **limma-voom** uses linear regression on transformed counts     \n",
    "    \n",
    "* **DESeq2** \"automagically\" handles independent filtering, using Independent Hypothesis Weighting ([IHW package](http://www.bioconductor.org/packages/IHW))\n",
    "* Speed (**limma-voom** is much faster and should definitely be used if you have hundreds of samples)\n",
    "* Approaches to hierarchical modeling, e.g. samples nested within groups\n",
    "    + limma-voom has a built-in feature for sample correlations (called *duplicateCorrelation*)\n",
    "* Re-use of core Bioconductor structures (DESeq2 extends *SummarizedExperiment*, edgeR and limma create their own)\n",
    "\n",
    "DEseq2 is generally a practical and simple pipeline to run.  However, for large datasets, especially if doing permutations for gene set enrichment analysis, limma with its voom method is far more practical. It is worth noting that for large sample sizes, it hardly matters which error model you use, owing to the Central Limit Theorem. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}