deseq2.Rmd

---
title: "DESeq2"
author: "Leo Lahti, Sudarshan Shetty et al. `r Sys.Date()`"
bibliography: 
- bibliography.bib
output:
  BiocStyle::html_document:
    number_sections: no
    toc: yes
    toc_depth: 4
    toc_float: true
    self_contained: true
    thumbnails: true
    lightbox: true
    gallery: true
    use_bookdown: false
    highlight: haddock
---
<!--
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{microbiome tutorial - comparisons}
  %\usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}  
-->


## Normalization and group-wise comparisons with DESeq2

Examples adapted from [Callahan et al. F1000 (2017)](https://f1000research.com/articles/5-1492/v2).

Load example data:

```{r deseq2-example, warning=FALSE, message=FALSE, fig.width=8, fig.height=5}
# Load libraries
library(microbiome)
library(ggplot2)
library(magrittr)
library(dplyr)
# Probiotics intervention example data 
data(dietswap) 

# Only check the core taxa to speed up examples
pseq <- core(dietswap, detection = 50, prevalence = 80/100)
```


Toy example, to be polished:

```{r deseq2-example2, warning=FALSE, message=FALSE, fig.width=8, fig.height=5}
library(phyloseq)
library(reshape2)
library(DESeq2)
library(knitr)
library(magrittr)
# Running the DESeq2 analysis
ds2 <- phyloseq_to_deseq2(pseq, ~ nationality)
dds <- DESeq(ds2)
res <- results(dds)
df <- as.data.frame(res)
df$taxon <- rownames(df)
df <- df %>% arrange(log2FoldChange, padj)

library(knitr)
print(head(kable((df))))
```


## Validating DESeq2 results

```{r deseq2-validate, warning=FALSE, message=FALSE, fig.width=8, fig.height=5, eval=FALSE}
# Identify top taxa based on standard ANOVA
source(system.file("extdata/check_anova.R", package = "microbiome"))
ano <- check_anova(pseq, "nationality");
ano$log2FC <- log2(ano$ave.AFR) - log2(ano$ave.AAM)
taxa.anova <- as.character(subset(ano, padj < 0.01 & abs(log2FC) > log2(2))$taxa)

# lowPick the top taxa based on DESEq2
taxa.deseq <- subset(res, padj < 0.01 & abs(log2FoldChange) > log2(2))$taxon

# Check overlap
# Most DESEq2 taxa are confirmed with ANOVA
library(gplots)

# Also the est p-values are well correlated (higher not so)
mf <- data.frame(df$padj, ano$padj)
p <- ggplot(mf, aes(x = log10(df$padj), y = log10(ano$padj))) +
       labs(x = 'DESeq2 adjusted p-value', y = 'ANOVA adjusted p-value') +
       geom_point()
print(p)
```


```{r plot_for_ven, echo=F, message=F, error=F, eval=FALSE}
library(venn) # Check UpSet plot instead
venn( list(ANOVA = taxa.anova,DESeq2 = taxa.deseq))
```