```{r include=FALSE, echo=FALSE, eval=TRUE}
options(width=300)
knitr::opts_chunk$set(echo = TRUE, eval=TRUE, fig.width = 7, fig.height = 5, fig.align = "center", size = "tiny", warning = FALSE, results = TRUE, message = FALSE, comment = "")
```

# Requirements

For people using their own laptop, install some R packages:

```{r eval=FALSE}
source("http://bioconductor.org/biocLite.R")
if (!require(DESeq2)) {
  biocLite("DESeq2", ask=FALSE)
}
if (!require(edgeR)) {
  biocLite("edgeR", ask=FALSE)
}
if (!require(gplots)) {
  biocLite("gplots", ask=FALSE)
}
```

# Context

 - Study of 48 WT yeast samples vs 48 Snf2 (KO) samples: 
 
 
    1. Gierliński,M., Cole,C., Schofield,P., Schurch,N.J., Sherstnev,A., Singh,V., Wrobel,N., Gharbi,K., Simpson,G., Owen-Hughes,T., et al. (2015) Statistical models for RNA-seq data derived from a two-condition 48-replicate experiment. Bioinformatics, 10.1093/bioinformatics/btv425. [pdf](http://bioinformatics.oxfordjournals.org/content/early/2015/08/17/bioinformatics.btv425.full.pdf)
    
    2. Schurch,N.J., Schofield,P., Gierliński,M., Cole,C., Sherstnev,A., Singh,V., Wrobel,N., Gharbi,K., Simpson,G.G., Owen-Hughes,T., et al. (2016) How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? [RNA, 10.1261/rna.053959.115](http://rnajournal.cshlp.org/content/early/2016/03/28/rna.053959.115).

 - RNA-Seq reads have been cleaned, mapped and counted to generated a count data matrix containing 7126 genes.

# Loading a data table

**R** enables to download data directly from the Web. Load the counts table containing one row per gene and one column per sample.

```{r}
URL.course <- "http://jvanheld.github.io/stats_avec_RStudio_EBA"

# Load the files content in an R data.frame
path.counts <- file.path(URL.course, "/practicals/yeast_2x48_replicates/data/counts.txt")
counts <- read.table(file=path.counts, sep="\t", header=TRUE, row.names=1)

path.expDesign <- "http://jvanheld.github.io/stats_avec_RStudio_EBA/practicals/yeast_2x48_replicates/data/expDesign.txt"
expDesign <- read.table(file=path.expDesign, sep="\t", header=TRUE)
```

# Checking the content of the count tables

```{r}
# look at the beginning of the counts and design table:
print(counts[1:4,1:4])
print(expDesign[1:4,])

# dimension of each table
dim(counts)
dim(expDesign)
```

```{r eval=FALSE}
View(counts)
View(expDesign)
```

# Factors and levels in R

Be careful to the reference level in the factor variables:

```{r}
print(expDesign$strain)
expDesign$strain <- relevel(expDesign$strain, "WT")
print(expDesign$strain)
```

Note the change: the first level is not WT. This genotype will thus be considered as reference for the differential analysis. 


# Basic description of the data: number of reads per sample

```{r}
barplot(colSums(counts)/1000000, main="Total number of reads per sample (million)")
```

# Basic description of the data: percentage of null counts per sample

```{r}
prop.null <- apply(counts, 2, function(x) 100*mean(x==0))
print(head(prop.null))
barplot(prop.null, main="Percentage of null counts per sample", las=1)
```

# Differential analysis with DESeq2

```{r message=FALSE}
# Load the DESeq2 R package
library(DESeq2)

# Build the DESeq2 main object with the count data + experimental design
dds0 <- DESeqDataSetFromMatrix(countData = counts, colData = expDesign, design = ~ strain)

# Print a summary of the data content
print(dds0)
```

# Get the results using two command lines

```{r message=TRUE}
# Detect differentially expressed genes
dds0 <- DESeq(dds0) 
res0 <- results(dds0)
print(res0)
print(summary(res0))
print(mcols(res0))
```

# How many replicates are required?

**Question**: how robust is the analysis to the number of samples? 

The original study contained 48 replicates per genotype, what happens if we select a smaller number?

**Test:** each attendee of this course could select a given number  (e.g. 3,4,5,10,15,20) and adapt the code on the next slide to run the analysis with that number of replicates per genotype. We then compare the results (number of genes, significance, ...).

# Sub-sampling: analysis using a few replicates

```{r}
nb.replicates <- 4 ## Each attendee chooses a number (3,4,5,10,15 or 20)
## Random sampling of the WT replicates (columns 1 to 48)
samples.WT <- sample(1:48, size=nb.replicates, replace=FALSE)
## Random sampling of the Snf2 replicates (columns 49 to 96)
samples.Snf2 <- sample(49:96, size=nb.replicates, replace=FALSE)
print(c(samples.WT, samples.Snf2))
dds <- DESeqDataSetFromMatrix(
  countData = counts[,c(samples.WT, samples.Snf2)], 
  colData = expDesign[c(samples.WT, samples.Snf2),], 
  design = ~ strain)
print(dds)
```

We now perform a differential analysis with DESeq2 step by step with some quality controls.

# Normalization

```{r fig.width=14}
dds <- estimateSizeFactors(dds)
print(sizeFactors(dds))
# effect of the normalization
par(mfrow=c(1,2))
boxplot(log2(counts(dds, normalized=FALSE)+1), main="Raw counts", col=rep(c("lightblue","orange"), each=nb.replicates))
boxplot(log2(counts(dds, normalized=TRUE)+1), main="Normalized counts", col=rep(c("lightblue","orange"), each=nb.replicates))
```

# Mean-variance relationship

```{r warning=FALSE}
mean.counts <- rowMeans(counts(dds))
variance.counts <- apply(counts(dds), 1, var)
plot(x=mean.counts, y=variance.counts, pch=16, cex=0.3, main="Mean-variance relationship", log="xy")
abline(a=0, b=1)
```

We observe over-dispersion in the data: the Poisson distribution is not adapted and we prefer the Negative-Binommial distribution.

# Principal Component Analysis (PCA)

```{r message=FALSE, fig.width=10}
# dispersions estimation
dds <- estimateDispersions(dds)
# make the data homoscedastic
rld <- rlog(dds) # alternative to the "Variance Stabilizing Transformation"
plotPCA(rld, intgroup="strain") # function from the DESeq2 package
```

# Statistical test for each gene

```{r}
dds <- nbinomWaldTest(dds)
res.DESeq2 <- results(dds, alpha=0.05, pAdjustMethod="BH")
print(head(res.DESeq2))
summary(res.DESeq2, alpha=0.05)
```

# Histogram of raw P-values

```{r}
hist(res.DESeq2$pvalue, col="lightblue", main="Histogram of raw P-values (DESeq2)", breaks=20, xlab="P-value")
```

# MA-plot

```{r}
plotMA(res.DESeq2, alpha=0.05) # function from the DESeq2 package
```

# Volcano-plot

Here we need to build the plot using R base functions:

```{r}
## Draw a volcano plot
plot(x=res.DESeq2$log2FoldChange, y=-log10(res.DESeq2$padj), 
     xlab="log2(Fold-Change)", ylab="-log10(adjusted P-value",
     col=ifelse(res.DESeq2$padj<=0.05, "red", "black"), main="Volcano plot")
grid() ## Add a grid to the plot
abline(v=0) ## plot the X axis
```

# Differential analysis using edgeR with a few replicates

```{r message=FALSE}
# load the edgeR R package
library(edgeR)
# create the edgeR main object: dge
dge <- DGEList(counts=counts[,c(samples.WT, samples.Snf2)], remove.zeros=FALSE)
dge$design <- model.matrix(~ strain, data=expDesign[c(samples.WT, samples.Snf2),])
print(dge)
```

# Normalization & dispersions estimation with edgeR

```{r}
# normalization
dge <- calcNormFactors(dge)
print(dge$samples$norm.factors)
# dispersions
dge <- estimateGLMCommonDisp(dge, dge$design)
dge <- estimateGLMTrendedDisp(dge, dge$design)
dge <- estimateGLMTagwiseDisp(dge, dge$design)
```

# Modeling and testing with edgeR

```{r}
fit <- glmFit(dge, dge$design)
print(dge$design)
lrt <- glmLRT(fit)
res.edgeR <- topTags(lrt,n=nrow(dge$counts),adjust.method="BH")$table
print(head(res.edgeR))
```

# Histogram of raw P-values

```{r}
hist(res.edgeR$PValue, col="lightblue", main="Histogram of raw P-values (edgeR)", breaks=20, xlab="P-value")
```

# Compare DESeq2 and edgeR results: normalization factors

```{r}
plot(x=sizeFactors(dds), y=dge$samples$norm.factors, xlab="DESeq2 size factors", ylab="edgeR normalization factors", pch=16)
```

The normalization/size factors computed by DESeq2 and edgeR are not comparable as they are used in a different manner in the statistical/mathematical models.

# Re-order the results according to the gene names

```{r}
print(head(res.DESeq2))
print(head(res.edgeR))
res.edgeR <- res.edgeR[order(rownames(res.edgeR)),]
res.DESeq2 <- res.DESeq2[order(rownames(res.DESeq2)),]
```

# Comparing log2(Fold-Change) estimations

```{r}
plot(x=res.edgeR$logFC, y=res.DESeq2$log2FoldChange, 
     pch=16, cex=0.4, xlim=c(-2,2), ylim=c(-2,2),
     xlab="edgeR log2(Fold-Change)", ylab="DESeq2 log2(Fold-Change)")
abline(a=0, b=1, col="red", lwd=4) # draw the y=x curve (y=a+b*x with a=0 and b=1)
abline(h=0, v=0) # horizontal and vertical line
```

# Comparing raw P-values

```{r}
plot(x=res.edgeR$PValue, y=res.DESeq2$pvalue, 
     pch=16, cex=0.4, xlab="edgeR raw P-value", ylab="DESeq2 raw P-value")
abline(a=0, b=1, col="red", lwd=4) # draw the y=x curve (y=a+b*x with a=0 and b=1)
```

# Number of differentially expressed genes

```{r}
# remember the number of replicates
print(nb.replicates)
# DESeq2
sum(res.DESeq2$padj <= 0.05, na.rm=TRUE)
# edgeR
sum(res.edgeR$FDR <= 0.05, na.rm=TRUE)
```

What's the behaviour of the number of differentially expressed genes according to the number of samples?

# Venn diagram

```{r message=FALSE}
library(gplots)
venn(list(DESeq2=rownames(res.DESeq2[which(res.DESeq2$padj <= 0.05),]), 
          edgeR=rownames(res.edgeR[which(res.edgeR$FDR <= 0.05),])))
```

Supplementary exercise: do the same plot for up- and down-regulated genes separately.

# What's edgeR or DESeq2-specific?

```{r}
DESeq2.genes <- rownames(res.DESeq2[which(res.DESeq2$padj <= 0.05),])
edgeR.genes <- rownames(res.edgeR[which(res.edgeR$FDR <= 0.05),])
# select a DESeq2 specific gene
spe.DESeq2 <- setdiff(DESeq2.genes, edgeR.genes)
summary(res.edgeR[spe.DESeq2,"FDR"])
# select a edgeR specific gene
spe.edgeR <- setdiff(edgeR.genes, DESeq2.genes)
summary(res.DESeq2[spe.edgeR,"padj"])
```

# DESeq2 results for one gene

```{r}
# plotCounts is a function from the DESeq2 R package
plotCounts(dds, gene="ycr017c", intgroup="strain", normalized=TRUE)
```

# Negative control

**Question:** how reliable are the genes declared "positive", i.e. differentially expressed?

In principle, a good DEG detection program should return a negative answer when the dataset contains no differentially expressed gene. 

**Test:** run DESeq2 and edgeR on a tricky dataset, built by selecting two subsets of replicates in the same genotype. Since the two sample sets are obtained in the same conditions, the analysis should return no DEG gene. 

This is a *negative control*: we build an experiment whose result is expected to be negative.

# Differential analysis under $H_0$

Here we perform a differential analysis in which we compare N WT samples vs N other WT samples.

```{r}
nb.replicates <- 10
samples.WT <- sample(1:48, size=2*nb.replicates, replace=FALSE)
print(samples.WT)
counts.H0 <- counts[,samples.WT]
expDesign.H0 <- expDesign[samples.WT,]
# add a fictive condition factor
expDesign.H0$condition <- factor(rep(c("A","B"), each=nb.replicates))
print(expDesign.H0)
dds.H0 <- DESeqDataSetFromMatrix(countData = counts.H0, colData = expDesign.H0, design = ~ condition)
```

# Differential analysis with DESeq2 under $H_0$

```{r message=FALSE}
dds.H0 <- DESeq(dds.H0)
res.H0 <- results(dds.H0)
summary(res.H0)
```

# sessionInfo

Here are the details of the R packages used to generate this document:

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