PERMANOVA.Rmd

---
title: "Group-wise comparisons of microbiome composition"
author: "Leo Lahti, Sudarshan Shetty et al."
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}  
-->


## PERMANOVA for community-level multivariate comparisons

PERMANOVA quantifies multivariate community-level differences between
groups.


Load example data:

```{r boxplot-example, warning=FALSE, message=FALSE}
# Load libraries
library(microbiome)
library(ggplot2)
library(dplyr)

# Probiotics intervention example data 
data(peerj32) # Source: https://peerj.com/articles/32/
pseq <- peerj32$phyloseq # Rename the example data

# Pick relative abundances (compositional) and sample metadata 
pseq.rel <- microbiome::transform(pseq, "compositional")
otu <- abundances(pseq.rel)
meta <- meta(pseq.rel)
```


## Visualize microbiome variation

Visualize the population density and highlight sample groups (probiotic treatment LGG  vs Placebo):

```{r comparisons_permanova_visu, error=FALSE, message=FALSE, warnings=FALSE}
p <- plot_landscape(pseq.rel, method = "NMDS", distance = "bray", col = "group", size = 3)
print(p)
```


## PERMANOVA significance test for group-level differences

Now let us evaluate whether the group (probiotics vs. placebo) has a
significant effect on overall gut microbiota composition. Perform PERMANOVA: 

```{r comparisons_permanova_analyse, message=FALSE, warnings=FALSE}
# samples x species as input
library(vegan)
permanova <- adonis(t(otu) ~ group,
               data = meta, permutations=99, method = "bray")

# P-value
print(as.data.frame(permanova$aov.tab)["group", "Pr(>F)"])
```


## Checking the homogeneity condition

Check that variance homogeneity assumptions hold (to ensure the reliability of the results):


```{r comparisons-permanova2, message=FALSE, warnings=FALSE}
# Note the assumption of similar multivariate spread among the groups
# ie. analogous to variance homogeneity
# Here the groups have signif. different spreads and
# permanova result may be potentially explained by that.
dist <- vegdist(t(otu))
anova(betadisper(dist, meta$group))
```

## Investigate the top factors

Show coefficients for the top taxa separating the groups

```{r permanova_top, fig.width=5, fig.height=5, message=FALSE, error=FALSE, warnings=FALSE}
coef <- coefficients(permanova)["group1",]
top.coef <- coef[rev(order(abs(coef)))[1:20]]
par(mar = c(3, 14, 2, 1))
barplot(sort(top.coef), horiz = T, las = 1, main = "Top taxa")
```