---
title: "Practical session about Normalization and differential analysis of RNAseq data - Effect of normalization and differential tests"
author: "Nathalie Vialaneix"
date: "Biostatistics PF training - March 29-30, 2023"
output:
  html_document:
    numbered_sections: yes
    theme: journal
    toc: yes
    toc_float:
      collapsed: no
    self_contained: yes
---

This file illustrates the normalization and differential analysis of RNAseq data
on a real life dataset. The data are provided by courtesy of the transcriptomic
platform of IPS2 and the name of the genes were shuffled. Hence the results are
not biologically interpretable. They are composed of three files, all included 
in the directory ```data```:

* ```D1-counts.txt``` contains the raw counts of the experiment (13 columns: 
the first one contains the gene names, the others correspond to 12 different
samples);

* ```D1-genesLength.txt``` contains information about gene lengths;

* ```D1-targets.txt``` contains information about the sample and the 
experimental design. 

In the rest of this file, we:

1/ first import and describe the data;

2/ then perform and compare different normalizations;

3/ finally perform a differential analysis using different methods and models.

We advice that you use the file by reading the **R** code while trying to make
sense of it. Then, you run it and analyze the results. The packages 
```devtools```, ```ggplot2```, ```gridExtra```, ```reshape2```, ```mixOmics```,
```RColorBrewer```, ```edgeR``` and ```VennDiagram``` are required to compile 
this file. Information about my system (including package and **R** versions) 
are provided at the end of this file (HTML format), in Section "Session 
information".

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

The packages are loaded with:

```{r loadPackages, message=FALSE, warning=FALSE}
library(ggplot2)
library(gridExtra)
library(reshape2)
library(mixOmics)
library(RColorBrewer)
library(FactoMineR)
library(factoextra)
library(edgeR)
library(VennDiagram)
library(devtools)
library(plotly)
```

## Settings

This command line is used to create a subdirectory to store results of the 
analyses:
```{r createDir}
dir.create("results")
```


## Data description and importation

The data used in this practical session correspond to 12 samples of RNAseq 
obtained from microdissections on plants. Plants have four different genotypes:
the wild type ("wt") plant and three types of mutants and are observed three 
times in the same conditions. for information Mutants 1 and 2 have a similar 
phenotype (more complicated leaves than the wild type) whereas mutant 3 has the 
opposite phenotype (simpler leaves than the wild type).

The datasets are included in the repository ```data``` located at the root of 
the material for this practical session. They can be loaded using:

```{r importData, message=FALSE}
raw_counts <- read.table("data/D1-counts.txt", header = TRUE, row.names = 1)
raw_counts <- as.matrix(raw_counts)
design <- read.table("data/D1-targets.txt", header = TRUE)
gene_lengths <- scan("data/D1-genesLength.txt")
```

### Raw counts

The dimensions of the raw count matrix are equal to:
```{r basicCountDesc}
dim(raw_counts)
```
which shows that the data contains `r ncol(raw_counts)` columns (samples) and 
`r nrow(raw_counts)` rows (genes). And a quick look is obtained by:
```{r headCounts}
head(raw_counts)
```
which shows that the row names are gene names and that the data are made of 
counts. Many of these counts are equal to 0. A basic exploratory analysis of the
count data is provided in the next section.

### Experimental design

The experimental design is described in the object ```design```:
```{r basicDesignDesc}
design
```
that contains 3 columns: the first one ```labels``` is the sample  name, the 
second one ```group``` is the group of the sample (wild type or one of three
types of mutant) and the last one ```replicat``` is the biological replicate 
number (that we will suppose to be corresponding from one group to another).


## Basic exploratory analysis of raw counts

We start the analysis by filtering out the genes for which no count has been 
found:
```{r filterCounts}
raw_counts_wn <- raw_counts[rowSums(raw_counts) > 0, ]
dim(raw_counts_wn)
```
The number of genes which have been filtered out is thus equal to 
`r nrow(raw_counts) - nrow(raw_counts_wn)`. 

It is often useful, to visualize the count distribution, to compute "pseudo 
counts", which are log-transformed counts:
```{r pseudoCounts}
pseudo_counts <- log2(raw_counts_wn + 1)
head(pseudo_counts)

df_raw <- melt(pseudo_counts, id = rownames(raw_counts_wn))
names(df_raw)[1:2]<- c("id", "sample")
df_raw$method <- rep("Raw counts", nrow(df_raw))
head(df_raw)
```
The object ```df_raw``` will be used later to compare the effect of different
normalization on the count distribution in the different samples.

### Count distribution

First, we provide an overview of the distribution of the first sample (without
the genes that have been filtered out) by plotting the histograms of raw counts
and pseudo counts:
```{r histoCount, warning=FALSE, message=FALSE}
df <- data.frame(rcounts = raw_counts_wn[ ,1], prcounts = pseudo_counts[ ,1])

p <- ggplot(data=df, aes(x = rcounts, y = ..density..)) +
  geom_histogram(fill = "lightblue") + theme_bw() +
  ggtitle(paste0("count distribution '", design$labels[1], "'")) +
  xlab("counts")

p2 <- ggplot(data=df, aes(x = prcounts, y = ..density..)) +
  geom_histogram(fill = "lightblue") + theme_bw() +
  ggtitle(paste0("count distribution - '", design$labels[1], "'")) +
  xlab(expression(log[2](counts + 1)))

grid.arrange(p, p2, ncol = 2)
```

This figure shows that the count distribution is highly skewed with a large 
proportion of genes with a count equal to 0 and a few number of genes with an
extreme count (more that $2^{15}$).

```{r saveHist, echo=FALSE, fig.show='hide', message=FALSE, results='hide'}
grid.arrange(p, p2, ncol = 2)
dev.print(png, file = "results/count-hist.png", width = 1000, height = 600, 
          res = 150)
```

An alternative representation can be obtained with boxplots (all samples, 
colored by group)
```{r boxplotCount, warning=FALSE, message=FALSE}
df <- cbind(design, t(pseudo_counts))
df <- melt(df, id.vars = c("labels", "group", "replicat"))

p <- ggplot(data = df, aes(x = labels, y = value, fill = group)) + 
  geom_boxplot() + theme_bw() + ggtitle("count distributions") + 
  xlab("sample") + ylab("counts") + 
  theme(axis.text.x = element_text(angle = 90))
p
```
```{r saveBoxplots, echo=FALSE, fig.show='hide', message=FALSE, results='hide'}
print(p)
dev.print(png, file = "results/count-boxplots.png", width = 1000, height = 600, 
          res = 150)
```


### Reproducibility between samples (MA plots)

```{r MAplots, message=FALSE, warning=FALSE}
df <- data.frame("A" = (pseudo_counts[ ,3] + pseudo_counts[ ,4])/2,
                 "M" = pseudo_counts[ ,3] - pseudo_counts[ ,4])
p <- ggplot(data = df, aes(x = A, y = M)) + geom_point(alpha = 0.2) + 
  theme_bw() + geom_smooth(colour = "red")
p
```
```{r saveMA, echo=FALSE, fig.show='hide', message=FALSE, results='hide'}
print(p)
dev.print(png, file = "results/count-MA.png", width = 1000, height = 600,
          res = 150)
```


### Heatmap

```{r cimPseudoCounts, fig.width=10, fig.height=10}
df <- as.matrix(dist(t(pseudo_counts)))
df
cim_color <- colorRampPalette(rev(brewer.pal(9, "Reds")))(16)
cim(df, color = cim_color, symkey = FALSE)
```

```{r saveCim, echo=FALSE, fig.show='hide', message=FALSE, results='hide'}
png(file = "results/count-heatmap.png", width = 800, height = 800, res = 150)
cim(df, color = cim_color, symkey = FALSE)
dev.off()
```


### PCA

Finally, a PCA is performed to identify the main sources of variability in the
dataset. Pseudo-counts are used to perform this PCA.
```{r PCA, cache = TRUE}
resPCA <- PCA(t(pseudo_counts), ncp = 12, graph = FALSE)
fviz_eig(resPCA)
```

The first component explains more variance than the other components which 
almost all equally reproduce variance. The projection of the samples on the 
first two PCs is obtained with:
```{r plotIndPCA, warning = FALSE}
p <- fviz_pca_ind(resPCA, col.ind = design$group)
p
plotly::ggplotly(p)
```

This figure shows that mutants 1 and 3 are well separated from the other wild
type but that wild type is not well separated from mutant 2. 

```{r savePCA, echo=FALSE, fig.show='hide', message=FALSE, results='hide'}
fviz_pca_ind(resPCA, col.ind = design$group)
dev.print(png, file = "results/raw_pca.png", width = 800, height = 800, 
          res = 150)
```


## Normalization

This section performs different normalization approaches for RNAseq data. These
normalization are performed using edgeR (TC, RPKM, UQ, TMM). The final count
distribution is compared to the raw count distribution in the last part of this
section.

### edgeR

In this section, the package edgeR is perform to compare the other normalization
approaches. To do so, a DGEList object has to be created from the count data:
```{r prepareDataEdgeR}
dge2 <- DGEList(raw_counts_wn)
dge2
```

The library sizes and the normalization factors for all samples are obtained by:
```{r edgeRLibSize}
dge2$samples
```

#### Total count

Normalization by Total Count (TC) is obtained directly by the function 
```cpm```. Pseudo-counts ($\log_2$ transformed counts) are computed and stored
in a data frame for later use.
```{r normTC}
norm_lib <- cpm(dge2)
pseudo_TC <- log2(cpm(dge2) + 1)

df_TC <- melt(pseudo_TC, id = rownames(raw_counts_wn))
names(df_TC)[1:2] <- c ("id", "sample")
df_TC$method <- rep("TC", nrow(df_TC))
```

#### RPKM

RPKM needs information about gene lengths which have to be passed to the 
function ```rpkm```. Again, pseudo-counts are computed and stored in a data 
frame for further use.
```{r RPKM}
gene_lengths_wn <- gene_lengths[rowSums(raw_counts) > 0]
pseudo_RPKM <- log2(rpkm(dge2, gene.length = gene_lengths_wn) + 1)

df_RPKM <- melt(pseudo_RPKM, id = rownames(raw_counts_wn))
names(df_RPKM)[1:2] <- c ("id", "sample")
df_RPKM$method <- rep("RPKM", nrow(df_RPKM))
```

#### Upper quartile

Upper quartile normalization is obtained with the function 
```calcNormFactors```. This function changes the variable ```norm.factors``` in
the DGEList object and thus also the way the functions ```cpm``` and ```rpkm```
are computing counts: now, a normalized library size, equal to the original 
library size multiplied by the normalization factor, is used to compute 
normalized counts: $\tilde{K}_{gj} = \frac{K_{gj}}{\mbox{modified LS}}*10^6$.
The normalization factors are computed using:
```{r normUQ}
dge2 <- calcNormFactors(dge2, method = "upperquartile")
dge2$samples
```

The previous formula is compared to the result of the function ```cpm``` in the
following command lines:
```{r compareUQ}
test_normcount <- sweep(dge2$counts * 10^6, 2,
                        dge2$samples$lib.size * dge2$samples$norm.factors,
                        "/")
range(as.vector(test_normcount - cpm(dge2)))
```
which shows no difference. Finally, pseudo counts are obtained and stored in a
data frame:
```{r pseudoUQ}
pseudo_UQ <- log2(cpm(dge2) + 1)

df_UQ <- melt(pseudo_UQ, id = rownames(raw_counts_wn))
names(df_UQ)[1:2] <- c ("id", "sample")
df_UQ$method <- rep("UQ", nrow(df_UQ))
```

```{r RLE}
dge2 <- calcNormFactors(dge2, method = "RLE")
dge2$samples
pseudo_RLE <- log2(cpm(dge2) + 1)

df_RLE <- melt(pseudo_RLE, id = rownames(raw_counts_wn))
names(df_RLE)[1:2] <- c ("id", "sample")
df_RLE$method <- rep("RLE", nrow(df_RLE))
```


#### TMM

TMM normalization works similarly as UQ and updates the value of the variable
```norm.factors```:
```{r TMM}
dge2 <- calcNormFactors(dge2, method = "TMM")
dge2$samples
```

Normalized pseudo counts are obtained with the function ```cpm``` and stored in
a data frame:
```{r pseudoTMM}
pseudo_TMM <- log2(cpm(dge2) + 1)

df_TMM <- melt(pseudo_TMM, id = rownames(raw_counts_wn))
names(df_TMM)[1:2] <- c ("id", "sample")
df_TMM$method <- rep("TMM", nrow(df_TMM))
```

### Comparison

This last section shows the comparison between all normalization methods (and
original raw data) for this data set. First, the distributions of all samples 
and all normalization methods are compared with boxplots and density plots:

```{r boxCompareNorm}
df_allnorm <- rbind(df_raw, df_TC, df_RPKM, df_UQ, df_TMM, df_RLE)
df_allnorm$method <- factor(df_allnorm$method,
                            levels = c("Raw counts", "TC", "RPKM", "TMM", "UQ", "RLE"))

p <- ggplot(data=df_allnorm, aes(x=sample, y=value, fill=method)) +
  geom_boxplot() + theme_bw() +
  ggtitle("Boxplots of normalized pseudo counts\n for all samples by normalization methods") + facet_grid(. ~ method) +
  ylab(expression(log[2] ~ (normalized ~ count + 1))) + xlab("") +
  theme(title = element_text(size=10), axis.text.x = element_blank(), 
        axis.ticks.x = element_blank())
print(p)
```

```{r saveBoxComparePlot, echo=FALSE, fig.show='hide', message=FALSE, results='hide'}
print(p)
dev.print(png, file = "results/norm_comparison.png", width = 1000, height = 600,
          res = 150)
```

```{r densityCompareNorm, fig.width=10}
p <- ggplot(data = df_allnorm, aes(x = value, colour = sample)) + 
  geom_density() + theme_bw() +
  ggtitle("Density of normalized pseudo counts\n for all samples by normalization methods") +
  facet_grid(. ~ method) + ylab("density") + xlab("") + 
  theme(title = element_text(size=10), axis.text.x = element_blank(), 
        axis.ticks.x = element_blank())
print(p)
```

Visually, DESeq, UQ and TMM seem to behave similarly and provide comparable 
distributions between samples. Finally, a PCA is performed on the pseudo counts
obtained after TMM normalization:

```{r normPCA, cache = TRUE}
resPCA <- PCA(t(pseudo_TMM), ncp  = 12, graph = FALSE)
fviz_eig(resPCA)
```

The first component explains more variance than the other components which 
almost all equally reproduce variance. The projection of the samples on the 
first two PCs is obtained with:
```{r plotIndNormPCA, warning = FALSE}
fviz_pca_ind(resPCA, col.ind = design$group)
```

This figure shows that a better separation between wild type samples and mutant
2 samples. On the contrary mutants 2 and 1 are less well separated than in raw
count data. This is a biologically plausible result since the organization of
the different mutants are concordant with their phenotypes. The normalization 
has thus improved the data quality.


```{r heatmapTMM, fig.width=10, fig.height=10}
df <- as.matrix(dist(t(pseudo_TMM)))
df
cim_color <- colorRampPalette(rev(brewer.pal(9, "Blues")))(16)
cim(df, color = cim_color, symkey = FALSE)
```


## Differential analysis

This section will compare the results of different types of approach to obtain
genes which are differentially expressed between the wild type plants and each 
of the mutants:

1. a standard NB exact test between two conditions;

2. a GLM with only the genotype effect;

3. a GLM with the plant (individual) and genotype effects.

All analyses are performed with the **R** package edgeR. A last section compares
the results obtained using the different methods.

### First approach: pairwise exact test

The following command lines loop over the different mutant types (stored in the
vector ```all_conditions```) and perform the following operations:

* first, a DGEList object is created with the wild type samples and the current
mutant type samples. The "wt" condition is chosen as the reference level. 
Normalization is performed using TMM;

* then, the functions ```estimateCommonDisp``` and ```estimateTagwiseDisp```
respectively estimate a common dispersion parameter and uses it as a prior to 
estimate gene specific dispersions $\phi_g$ for all $g$. The "Biological 
Coefficient of Variation" plot is also displayed which displays $\phi_g$ versus
the average (normalized) count for all genes;

* finally, an exact test is performed with the function ```exactTest``` to find
genes which are differentially expressed between the two conditions (wild type
and mutant). The output of the function is printed and p-values and adjusted
p-values (Benjamini and Hochberg approach) are stored in
```pvals_pairwiseExact```. The function ```topTags``` is used to print the genes
with the smallest p-values and the function ```decideTestsDGE``` is used to
extract genes whose adjusted p-value is smaller than 5%. The names of those
genes as well as the fact that they are up/down regulated in the mutant type are
saved in ```DGE_pairwiseExact```.

```{r twoGroupTest, cache=TRUE}
cur_counts <- raw_counts_wn[ ,1:6]
cur_design <- design[1:6, ]
cur_dge <- DGEList(cur_counts, group = cur_design$group, remove.zeros = TRUE)
cur_dge <- calcNormFactors(cur_dge, method = "TMM")
pseudo_counts <- log2(cpm(cur_dge) + 1)

# estimate dispersion
cur_dge <- estimateDisp(cur_dge)
plotBCV(cur_dge, main = paste0("BCV plot"))

# perform test
res_et <- exactTest(cur_dge)
p_adjust <- p.adjust(res_et$table$PValue, method = "BH")
df <- data.frame("pvalue" = c(res_et$table$PValue, p_adjust),
                 "type" = rep(c("raw", "adjusted"), each = nrow(res_et$table)))
df$type <- factor(df$type, levels = c("raw", "adjusted"))
p <- ggplot(df, aes(x = pvalue)) + geom_histogram() + theme_bw() +
  facet_wrap(. ~ type) + xlab("p-value")
p

# extracting results
to_export <- topTags(res_et, n = nrow(res_et$table))
write.table(to_export, file = "results/mesjolisresultats.csv", sep = ";",
            row.names = TRUE, col.names = TRUE, dec = ",")
topTags(res_et, adjust.method = "bonferroni")
cur_res <- decideTestsDGE(res_et, adjust.method = "BH", p.value = 0.05)
cur_res
sel_deg <- which(cur_res[ ,1] != 0)
length(sel_deg)
DEG_pairwiseExact <- rownames(cur_res)[sel_deg]
head(DEG_pairwiseExact)

# explore results
df <- data.frame("logCPM" = res_et$table$logCPM,
                 "logFC" = res_et$table$logFC,
                 "selected" = unname(ifelse(cur_res[ ,1] != 0, "DEG", "not DE")))
p <- ggplot(df, aes(x = logCPM, y = logFC, colour = selected)) + 
  geom_point(alpha = 0.5) + theme_bw() + 
  scale_colour_manual(name = "test result", values = c("red", "black"))
p

df$adjpval <- p.adjust(res_et$table$PValue, "BH")
p <- ggplot(df, aes(x = logFC, y = -log10(adjpval), colour = selected)) +
  geom_point(alpha = 0.5) + theme_bw() +
  scale_colour_manual(name = "test result", values = c("red", "black"))
p
```

```{r CIMdiff, fig.width=10, fig.height=10}
df <- pseudo_counts[sel_deg, ]
df <- scale(t(as.matrix(df)))
cim_color <- colorRampPalette(rev(brewer.pal(11, "RdYlGn")))(256)
res_cim <- cim(df, color = cim_color, cluster = "col", col.names = FALSE, 
               dist.method = c("euclidean", "correlation"))
```
Extracting dendrogram associated with these genes are groups is easy with:
```{r HACgenes}
res_hac <- as.hclust(res_cim$ddc)
plot(res_hac, labels = FALSE, xlab = "", sub = "", main = "Groups of genes")
rect.hclust(res_hac, k = 3, border = "red")
genes_group <- cutree(res_hac, k = 3)
head(genes_group)
table(genes_group)
```


### Second approach: overall GLM tests with group effect

In this section, we fit a generalized linear model to account for the group
effect. The model writes $K_{gj} \sim \mbox{NB}(\mu_{gj}, \phi_g)$ with
$\mathbb{E}(\log K_{gj}) = \log(s_j) + \log(\lambda_{gj})$. Here, $j$ denotes
one of the four genotypes (wild type or mutant 1, 2 and 3) and $\lambda_{gj}$ 
is explained by: $\log(\lambda_{g,\textrm{mutant}}) = \lambda_{g0} + \beta_{g1} 
\mathbb{I}_{\textrm{mutant}_1} + \beta_{g2} \mathbb{I}_{\textrm{mutant}_2} +
\beta_{g3} \mathbb{I}_{\textrm{mutant}_3}$.  In this model, the wild type plant
is the reference genotype and thus, testing if the second (resp., the third, 
the fourth) coefficient in the model is equal to zero is equivalent to find DEGs
between the wild type and mutant 1 (resp. 2, 3).

First, a DGEList object is created with all samples and is normalized by TMM:
```{r GLM1}
# create dataset with 'wt' and current mutant and normalize (TMM)
cur_dge <- DGEList(raw_counts_wn, group = design$group)
cur_dge <- calcNormFactors(cur_dge, method = "TMM")
```

Then, a design matrix is created: the group condition is stored in a variable 
```group``` in which the reference level is set to "wt". Then, this variable is
used in the function ```model.matrix``` to obtain the design matrix:
```{r GLM1Fit, cache=TRUE}
# design matrix
group <- relevel(as.factor(design$group), ref = "wt")
design_matrix <- model.matrix(~ group)
design_matrix

# estimate dispersion
cur_dge <- estimateDisp(cur_dge, design_matrix)
plotBCV(cur_dge, main = paste0("BCV plot"))

# fit model and perform test
fit <- glmFit(cur_dge, design_matrix)
res_GLM1 <- glmLRT(fit, coef = 2)
table_res <- res_GLM1$table
p_adjust <- p.adjust(table_res$PValue, method = "BH")
table_res$padj <- p_adjust

df <- data.frame("pvalue" = c(res_GLM1$table$PValue, p_adjust),
                 "type" = rep(c("raw", "adjusted"), each = nrow(res_GLM1$table)))
df$type <- factor(df$type, levels = c("raw", "adjusted"))
p <- ggplot(df, aes(x = pvalue)) + geom_histogram() + theme_bw() +
  facet_wrap(. ~ type) + xlab("p-value")
p

# extracting results
cur_res <- decideTestsDGE(res_GLM1, adjust.method = "BH", p.value = 0.05)
cur_res
sel_deg <- which(cur_res[ ,1] != 0)
length(sel_deg)
DEG_GLM1 <- rownames(cur_res)[sel_deg]
head(DEG_GLM1)

# explore results
df <- data.frame("logCPM" = res_GLM1$table$logCPM,
                 "logFC" = res_GLM1$table$logFC,
                 "selected" = unname(ifelse(cur_res[ ,1] != 0, "DEG", "not DE")))
p <- ggplot(df, aes(x = logCPM, y = logFC, colour = selected)) + 
  geom_point(alpha = 0.5) + theme_bw() + 
  scale_colour_manual(name = "test result", values = c("red", "black"))
p

df$adjpval <- p.adjust(res_GLM1$table$PValue, "BH")
p <- ggplot(df, aes(x = logFC, y = -log10(adjpval), colour = selected)) +
  geom_point(alpha = 0.5) + theme_bw() +
  scale_colour_manual(name = "test result", values = c("red", "black"))
p
```

Compared to what was found with the Exact test, the number of obtained DEGs is 
of the same order (only a bit larger). A more precise comparison is made in the 
last part of this section.

The difference between mutant 1 and mutant 3 can be assessed using contrasts
on the same fitted model:

```{r GLM1Contrasts}
# fit model and perform test
contrasts <- c(0, 1, 0, -1)
res <- glmLRT(fit, contrast = contrasts)

# extracting results
cur_res <- decideTestsDGE(res, adjust.method = "BH", p.value = 0.05)
sel_deg <- which(cur_res[ ,1] != 0)
length(sel_deg)
DEG_GLM1c <- rownames(cur_res)[sel_deg]
head(DEG_GLM1c)
```


### Third approach: overall GLM tests with group and replicate effects

The previous model does not take into account all information about the 
experimental design. Only the group effect has been tested whereas an additional
effect might influence the expression level: the replicat effect. In this 
section we fit a model with an additive contribution of both effects. This 
problem is equivalent to introduce a (fixed) individual effect in the model. In 
this situation, the average level of $\log$ normalized counts is expressed as:
$\log(\lambda_{g,\textrm{rep},\textrm{mutant}}) = \lambda_{g0} + \beta_{g1}
\mathbb{I}_{\textrm{rep}_1} + \beta_{g2} \mathbb{I}_{\textrm{rep}_2} + 
\beta_{g3} \mathbb{I}_{\textrm{mutant}_1} + \beta_{g4}
\mathbb{I}_{\textrm{mutant}_2} + \beta_{g5} \mathbb{I}_{\textrm{mutant}_3}$ and 
DEGs for each of the three mutant types are obtained by testing coefficients 
$\beta_{g3}$, $\beta_{g4}$ and $\beta_{g5}$.

The method is performed similarly as before: a DGEList is created and 
normalization factors are obtained with the TMM approach. Then a design matrix
that incorporates the two factors in an additive model is created:
```{r GLM2, cache = TRUE}
# create dataset with 'wt' and current mutant and normalize (TMM)
cur_dge <- DGEList(raw_counts_wn, group = group)
cur_dge <- calcNormFactors(cur_dge, method = "TMM")

# create design matrix
design_matrix <- model.matrix(~ design$replicat + group)
design_matrix
```

The dispersion is then estimated using the same approach as in the previous 
section:
```{r GLM2Disp, cache=TRUE}
cur_dge <- estimateDisp(cur_dge, design_matrix)
plotBCV(cur_dge, main = paste0("BCV plot"))
```

And finally, the model is fitted with ```glmFit``` and tests are performed on 
the coefficients of interest with ```glmLRT```. Results extracted with `topTags`
and ```decideTestsDGE``` are saved in ```pvals_GLM2``` and ```listDEG_GLM2```:
```{r GLM2Fit}
# GLM fit
fit <- glmFit(cur_dge, design_matrix)
res_GLM2 <- glmLRT(fit, coef = 4)
df <- data.frame("pvalue" = c(res_GLM2$table$PValue, 
                               p.adjust(res_GLM2$table$PValue, method = "BH")),
                 "type" = rep(c("raw", "adjusted"), each = nrow(res_GLM2$table)))
df$type <- factor(df$type, levels = c("raw", "adjusted"))
p <- ggplot(df, aes(x = pvalue)) + geom_histogram() + theme_bw() +
  facet_wrap(. ~ type) + xlab("p-value")
p

# extracting results
cur_res <- decideTestsDGE(res_GLM2, adjust.method = "BH", p.value = 0.05)
cur_res
sel_deg <- which(cur_res[ ,1] != 0)
length(sel_deg)
DEG_GLM2 <- rownames(cur_res)[sel_deg]
head(DEG_GLM2)

# explore results
df <- data.frame("logCPM" = res_GLM2$table$logCPM,
                 "logFC" = res_GLM2$table$logFC,
                 "selected" = unname(ifelse(cur_res[ ,1] != 0, "DEG", "not DE")))
p <- ggplot(df, aes(x = logCPM, y = logFC, colour = selected)) + 
  geom_point(alpha = 0.5) + theme_bw() + 
  scale_colour_manual(name = "test result", values = c("red", "black"))
p

df$adjpval <- p.adjust(res_GLM2$table$PValue, "BH")
p <- ggplot(df, aes(x = logFC, y = -log10(adjpval), colour = selected)) +
  geom_point(alpha = 0.5) + theme_bw() +
  scale_colour_manual(name = "test result", values = c("red", "black"))
p
```


### Comparison

In this last part, a Venn diagram is displayed for the results (unique genes for
mutant type 1 versus reference) of the three approaches:
```{r venn}
vd <- venn.diagram(x=list("Exact test" = DEG_pairwiseExact,
                          "GLM\n group" = DEG_GLM1, 
                          "GLM\n group + replicate" = DEG_GLM2), 
                   fill = brewer.pal(3, "Set3"), 
                   cat.col = c("darkgreen", "black", "darkblue"),
                   cat.cex = 1.5, fontface = "bold", filename = NULL)
grid.draw(vd)
```

This chart yields to several conclusions:

* first, a large number of genes ($201$ out of approximately 400) are in common 
between all three approaches. As expected, the two more consistent methods are 
the two GLM approaches (with an additional number of $57$ DEGs in common and 
$46+28+9+19 = 102$ DEGs that are specific of a given model). The Exact Test has 
42 DEGs that are not found by one of the other two approaches and only 
$19+9=28$ DEGs in common with one of the two GLM models only;

* the 46 DEGs that were found exclusively by GLM with two additive factors might
be interesting because they are DEGs which can not be found if the genotype 
effect is not included in the model. 

```{r saveVDPlot, echo=FALSE, fig.show='hide', message=FALSE, results='hide'}
png(file = "results/venn_deg.png", width = 600, height = 600, res = 150)
grid.draw(vd)
dev.off()
system("rm VennDiagram*.log")
```

## Session information

To obtain results similar to mines, make sure that your session has identical
features than the one used to compile this document:

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