---
title: 'DESeq2: Results'
output:
    html_document:
      number_sections: true
      toc: true
      theme: cosmo
      highlight: tango
---

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

library(DESeq2)
library(pheatmap)
library(genefilter)
```

# Import workspace

```{r eval=TRUE}
fcp = file.copy("DESEQ_WORKSPACE", "deseq.RData")
load("deseq.RData")
```

# Results {.tabset}

## Result table

```{r}
group = colnames(sample_table)[CONTRAST_GROUP]
res <- results(dds, contrast = c(group, 'TREATMENT_LEVEL', 'CONDITION_LEVEL'))
datatable(as.data.frame(res), style="bootstrap", filter = 'top',
          class="table-condensed", options = list(dom = 'tp', scrollX = TRUE))
```

## Result summary

```{r}
summary(res)
```


# MA-plot {.tabset}

## Shrinked with `lfcShrink()` function

```{r eval=FALSE}
shrink_res = DESeq2::lfcShrink(dds, contrast = c(group, 'TREATMENT_LEVEL', 'CONDITION_LEVEL'), res=res)
plotMA(shrink_res)
```

## Shrinked with Bayesian procedure

```{r}
plotMA(res)
```


# Histogram of p values

```{r}
hist(res$pvalue[res$baseMean > 1], breaks = 0:20/20,
     col = "grey50", border = "white", main = "",
     xlab = "Mean normalized count larger than 1")
```


# Gene clustering

```{r}
group_index = as.numeric(strsplit("CLUSTERING_GROUPS", ',')[[1]])
clustering_groups = colnames(sample_table)[group_index]

topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)
mat  <- assay(rld)[ topVarGenes, ]
mat  <- mat - rowMeans(mat)
annotation_col <- as.data.frame(colData(rld)[, clustering_groups])
colnames(annotation_col) = clustering_groups
rownames(annotation_col) = colnames(mat)
pheatmap(mat, annotation_col = annotation_col)
```