PlotDiversity.Rmd

---
title: "Plot alpha diversity"
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
---
<!--
  %\VignetteIndexEntry{microbiome tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}
-->

## Diversity plots  

Using the `alpha` function in `microbiome R packge` you can calculate a wide variaty of diversity indices. Comparison and visualising group based differecences or similarities is also important. Here, we show steps from calculating diversity indices using microbiome R package and visualising the differences and/or similarities between groups. A useful R package is [ggpubr](http://www.sthda.com/english/rpkgs/ggpubr/). If you have not installed it, please install it.     

Load libraries and data.

```{r Ad-plot, warning=FALSE, message=FALSE}

library(microbiome)
library(ggpubr)
library(knitr)
library(dplyr)

data(dietswap)
pseq <- dietswap
```


## Alpha diversity

This returns a table with selected diversity indicators. Check a separate page on [Alpha](Diversity.html) for other functions.

```{r Ad-plot1, warning=FALSE, message=FALSE, results="asis"}

ps1 <- prune_taxa(taxa_sums(pseq) > 0, pseq)

tab <- microbiome::alpha(ps1, index = "all")
kable(head(tab))

```


## Prepare data for vizualisation
Now, get the metadata (sample_data) from the `phyloseq` object

```{r Ad-plot2, warning=FALSE, message=FALSE, results="asis"}

ps1.meta <- meta(ps1)
kable(head(ps1.meta))

```

Add the diversity table to metadata

```{r Ad-plot3, warning=FALSE, message=FALSE}

ps1.meta$Shannon <- tab$diversity_shannon 
ps1.meta$InverseSimpson <- tab$diversity_inverse_simpson
```

Let's say we want to compare differences in Shannon index between bmi group of the study subjects.  

```{r Ad-plot4, warning=FALSE, message=FALSE}

# create a list of pairwise comaprisons
bmi <- levels(ps1.meta$bmi_group) # get the variables

# make a pairwise list that we want to compare.
bmi.pairs <- combn(seq_along(bmi), 2, simplify = FALSE, FUN = function(i)bmi[i])

print(bmi.pairs)
```

Sometimes that variables can be stored as characters and sometime as factors. In such a senario `levels()` may return an empty vector. A work around for this can be found [here](https://github.com/microbiome/microbiome/issues/143).   


## Violin plot  

Using `ggpubr` a violin plot will be created

```{r Ad-plot5, warning=FALSE, message=FALSE, fig.width=8, fig.height=6, eval=FALSE}
#ps1.meta$'' <- alpha(ps1, index = 'shannon')
p1 <- ggviolin(ps1.meta, x = "bmi_group", y = "Shannon",
 add = "boxplot", fill = "bmi_group", palette = c("#a6cee3", "#b2df8a", "#fdbf6f")) 
print(p1)
```
![violin plot](Rplotvio1.jpeg)

## Statistics  

Pairwise comparision using non-parametric test (Wilcoxon test).

```{r Ad-plot6, warning=FALSE, message=FALSE, eval=FALSE}
p1 <- p1 + stat_compare_means(comparisons = bmi.pairs) 
print(p1)
```
![violin for comparison](Rplotvio2.jpeg)


For more information and useful tips and suggestions check the [Statistical tools for high-throughput data analysis](http://www.sthda.com/english/rpkgs/ggpubr/).