---
title: "Introduction using limma or edgeR"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{limma}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## pre-load to avoid load messages in report
library(Glimma)
library(limma)
library(edgeR)
```

## Introduction

In this vignette we present the basic features of Glimma. 
Glimma is an interactive R widget for creating plots for differential expression analysis, created using the Vega and htmlwidgets frameworks. The created plots can be embedded in R Markdown, or exported as standalone HTML documents.
The data presented here is slightly modified from the [RNAseq123](https://bioconductor.org/packages/release/workflows/html/RNAseq123.html) workflow and only a single contrast has been performed for simplicity. We can use either limma or edgeR to fit the models and they both share upstream steps in common.

To begin, the DGEList object from the workflow has been included with the package as internal data.

```{r}
library(Glimma)
library(limma)
library(edgeR)

dge <- readRDS(system.file("RNAseq123/dge.rds", package = "Glimma"))
```

## MDS Plot

The multidimensional scaling (MDS) plot is frequently used to explore differences in samples. When data has been MDS transformed, the first two dimensions explain the greatest variance between samples, and the amount of variance decreases monotonically with increasing dimension.

The Glimma MDS plot contains two main components:

1. a plot showing two MDS dimensions, and
2. a plot of the eigenvalues of each dimension

The Glimma MDS allows different dimensions to be plotted against each other, and for the colours of the points to be changed based on predefined factors. The grouping variables are taken from the `samples` component of `DGEList` objects used in `limma` and `edgeR`.

```{r}
glimmaMDS(dge)
```

### Interactions with the plot
In the plot above, try:

+ Scaling the points by library size (lib_size).
+ Changing the colour of points by group using the colour_by field.
+ Changing the colour scheme using to colour points using the colourscheme field.
+ Altering the shape of points by sample sequencing lane using the shape_by field.
+ Changing the dimensions plotted on the x-axis (x_axis) to dim2 and y-axis (y_axis) to dim3.
+ Saving the plots in either PNG or SVG formats using the "Save Plot" button.

### Modifications to the plot
Some customisations to the plot include:

+ `glimmaMDS(dge, width=1200, height=1200)`, which will adjust the dimensions in pixels of the created widget - default width and height is 920px.
+ `glimmaMDS(dge, continuous.color=TRUE)`, which specifies that continuous colour schemes should be used - useful for when a large number of differential selections are required.
+ `glimmaMDS(dge, groups=[vector or data frame])`, which allows changing the associated sample information such as experimental groups - this information is displayed in mouseover tooltips and can be used to adjust the plot using `scale_by`, `colour_by` and `shape_by`.


## MA Plot

The MA plot is a visualisation that plots the log-fold-change between experimental groups (M) against the mean expression across all the samples (A) for each gene.

The Glimma MA plot contains two main components:

1. a plot of summary statistics across all genes that have been tested, and
2. a plot of gene expression from individual samples for a given gene

The second plot shows gene expression from the last selected sample, which can be selected from the table or directly from the summary plot.

To create this plot we first need to run differential expression (DE) analysis for our data. We load in design and contrast matrices generated from the RNAseq123 workflow.

```{r}
design <- readRDS(
  system.file("RNAseq123/design.rds", package = "Glimma"))
contr.matrix <- readRDS(
  system.file("RNAseq123/contr.matrix.rds", package = "Glimma"))
```

### Using limma

We fit our DE analysis using `voom`, this leaves us with an object that contains test statistics for each gene.

```{r}
v <- voom(dge, design)
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
efit <- eBayes(vfit)
```

### Using edgeR

```{r}
dge <- estimateDisp(dge, design)
gfit <- glmFit(dge, design)
glrt <- glmLRT(gfit, design, contrast = contr.matrix)
```

The MA plot can then be created using the fitted object containing the statistics about the genes, and the `dge` object containing information about the samples and raw counts.

```{r}
glimmaMA(efit, dge = dge) # swap efit for glrt to use edgeR results
```

### Interactions with the plot
In the plot above, try:

+ Clicking points in the summary plot or rows in the table to plot the gene expression of the selection.
    + Clicking genes in the table after selecting individual points will remove the previous selection.
+ Searching for individual genes using the search box. The search results are displayed in the table.
    + If genes are currently selected, the search box will not function.
+ Setting a maximum value for the y-axis of the expression plot using the max_y_axis field. 
    + This allows for comparison of gene expression between genes on a comparable scale.
+ Saving the currently selected genes using the Save Data dropdown. 
    + From here, you can also choose to save the entire table.
+ Saving the summary plot or expression plot in either PNG or SVG formats, using the "Save Data" dropdown.

### Modifications to the plot
Some customisations to the plot include:

+ `glimmaMA(efit, dge=dge, width=1200, height=1200)`, which will adjust the dimensions in pixels of the created widget
    + Default width and height is 920px.
+ `glimmaMA(efit, dge=dge, continuous.color=TRUE)`, which specifies that continuous colour schemes should be used
    + Useful for when a large number of differential selections are required.
+ `glimmaMA(efit, dge=dge, groups=[vector or data frame])`, which allows changing the associated sample information such as experimental groups
    + This information is displayed in mouseover tooltips and can be used to adjust the plot using scale_by, colour_by and shape_by.
+ `glimmaMA(efit, dge=dge, status.cols=c("powderblue", "seashell", "salmon")`, which customises the colours associated with the status of each gene
    + These need to be valid CSS colour strings.
+ `glimmaMA(efit, dge=dge, sample.cols=colours)`, which colours each sample based on the character vector of valid CSS colour strings `colours`
    + This vector needs to be of length `ncol(dge$counts)` or `ncol(counts)` if specified.

## Saving widgets

The plots created are automatically embedded into Rmarkdown reports, but having many interactive plots can significantly slow down the page. It is instead recommended to save the plots and link to them via markdown hyperlinks. Plots can be saved either by providing the `html` argument a filename, or by using `htmlwidgets::saveWidget`, which also provides further customisation options.

```{r, eval = FALSE}
# creates ma-plot.html in working directory
# link to it in Rmarkdown using [MA-plot](ma-plot.html)
htmlwidgets::saveWidget(glimmaMA(efit, dge = dge), "ma-plot.html")
```

## Session Info

```{r}
sessionInfo()
```