CoremicrobiotaAmplicon.Rmd

---
title: "Core microbiome analysis for Amplicon data"
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 - core}
  %\usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}  
-->

## Make phyloseq object

This tutorial is useful for analysis of output files from [(Mothur)](https://www.mothur.org/), [(QIIME or QIIME2)](https://qiime2.org/) or any tool that gives a biom file as output. There is also a simple way to read comma seperated (*.csv) files.  

Simple comma seperated files:  

```{r, read-simple-csv-otu-tables, warning=FALSE, message=FALSE, eval=FALSE}
library(microbiome)


otu.file <-
    system.file("extdata/qiita1629_otu_table.csv",
        package='microbiome')

tax.file <- system.file("extdata/qiita1629_taxonomy_table.csv",
        package='microbiome')

meta.file <- system.file("extdata/qiita1629_mapping_subset.csv",
        package='microbiome')

pseq.csv <- read_phyloseq(
          otu.file=otu.file, 
          taxonomy.file=tax.file, 
          metadata.file=meta.file, type = "simple")
```

Biom file:  

```{r, read-otu-biom, eval=FALSE}

# Read the biom file
biom.file <- 
  system.file("extdata/qiita1629.biom", 
              package = "microbiome")

# Read the mapping/metadata file
 meta.file <- 
  system.file("extdata/qiita1629_mapping.csv", 
              package = "microbiome")
# Make phyloseq object
pseq.biom <- read_phyloseq(otu.file = biom.file, 
                         metadata.file = meta.file, 
                         taxonomy.file = NULL, type = "biom")
```


Mothur shared OTUs and Consensus Taxonomy:  

```{r, read-otu-mothur, eval=FALSE}
otu.file <- system.file(
 "extdata/Baxter_FITs_Microbiome_2016_fit.final.tx.1.subsample.shared",
    package='microbiome')

tax.file <- system.file(
 "extdata/Baxter_FITs_Microbiome_2016_fit.final.tx.1.cons.taxonomy",
    package='microbiome')

meta.file <- system.file(
 "extdata/Baxter_FITs_Microbiome_2016_mapping.csv",
    package='microbiome')
 
pseq.mothur <- read_phyloseq(otu.file=otu.file,
        taxonomy.file =tax.file,
        metadata.file=meta.file, type = "mothur")
print(pseq.mothur)
```

Now, we proceed to core microbiota analysis.

## Core microbiota analysis  

Here the data from [Caporaso, J. Gregory, et al. "Moving pictures of the human microbiome." Genome biology 12.5 (2011): R50.](https://genomebiology.biomedcentral.com/articles/10.1186/gb-2011-12-5-r50?report=reader) will be used which is stored as example in [jeevanuDB](https://github.com/microsud/jeevanuDB) 

```{r, core-microbiota-amplicon-data, eval=TRUE}
# install
# install.packages("devtools")
# devtools::install_github("microsud/jeevanuDB")

# check the data 
library(jeevanuDB)
ps <- moving_pictures
table(meta(ps)$sample_type, meta(ps)$host_subject_id)
# Filter the data to include only gut samples from M3 subject
ps.m3 <- subset_samples(ps, sample_type == "stool" & host_subject_id == "M3") 
print(ps.m3)
# keep only taxa with positive sums
ps.m3 <- prune_taxa(taxa_sums(ps.m3) > 0, ps.m3)
print(ps.m3)

# Calculate compositional version of the data
# (relative abundances)
ps.m3.rel <- microbiome::transform(ps.m3, "compositional")
```

Output of deblur/dada2 will most likely have seqs as rownames instead of OTU ids or taxa names
```{r}

taxa_names(ps.m3.rel)[1:3]

```

We can change it to ASVIDs 

```{r}
library(Biostrings)
dna <- Biostrings::DNAStringSet(taxa_names(ps.m3.rel))
names(dna) <- taxa_names(ps.m3.rel)
ps.m3.rel <- merge_phyloseq(ps.m3.rel, dna)
taxa_names(ps.m3.rel) <- paste0("ASV", seq(ntaxa(ps.m3.rel)))
# now check again
taxa_names(ps.m3.rel)[1:3]
```

### Core microbiota analysis

If you only need the names of the core taxa, do as follows. This returns the taxa that exceed the given prevalence and detection thresholds. 

```{r core-members, message=FALSE, warning=FALSE, eval = FALSE}
core.taxa.standard <- core_members(ps.m3.rel, detection = 0.0001, prevalence = 50/100)

core.taxa.standard
```


A full phyloseq object of the core microbiota is obtained as follows:

```{r core-data, message=FALSE, warning=FALSE, eval=TRUE}
pseq.core <- core(ps.m3.rel, detection = 0.0001, prevalence = .5)
```


Retrieving the associated taxa names from the phyloseq object:

```{r core-taxa, message=FALSE, warning=FALSE, eval=TRUE}
core.taxa <- taxa(pseq.core)
class(core.taxa)
# get the taxonomy data
tax.mat <- tax_table(pseq.core)
tax.df <- as.data.frame(tax.mat)

# add the OTus to last column
tax.df$OTU <- rownames(tax.df)

# select taxonomy of only 
# those OTUs that are core memebers based on the thresholds that were used.
core.taxa.class <- dplyr::filter(tax.df, rownames(tax.df) %in% core.taxa)
knitr::kable(head(core.taxa.class))
```


## Core visualization

### Core line plots

Determine core microbiota across various abundance/prevalence
thresholds with the blanket analysis [(Salonen et al. CMI, 2012)](http://onlinelibrary.wiley.com/doi/10.1111/j.1469-0691.2012.03855.x/abstract) based on various signal and prevalences.

```{r core2, warning=FALSE, eval=TRUE}
# With compositional (relative) abundances
det <- c(0, 0.1, 0.5, 2, 5, 20)/100
prevalences <- seq(.05, 1, .05)

plot_core(ps.m3.rel, prevalences = prevalences, 
          detections = det, plot.type = "lineplot") + 
  xlab("Relative Abundance (%)") + 
  theme_bw()

```

### Core heatmaps

This visualization method has been used for instance in [Intestinal microbiome landscaping: Insight in community assemblage and implications for microbial modulation strategies](https://academic.oup.com/femsre/article/doi/10.1093/femsre/fuw045/2979411/Intestinal-microbiome-landscaping-insight-in#58802539). Shetty et al. _FEMS Microbiology Reviews_ fuw045, 2017.

Note that you can order the taxa on the heatmap with the order.taxa argument.

```{r core-example3, warning=FALSE, eval=TRUE}
# Core with compositionals:
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-4), log10(.2), length = 10)

# Also define gray color palette
gray <- gray(seq(0,1,length=5))
p1 <- plot_core(ps.m3.rel, 
                plot.type = "heatmap", 
                colours = gray,
                prevalences = prevalences, 
                detections = detections, min.prevalence = .5) +
    xlab("Detection Threshold (Relative Abundance (%))")
p1 <- p1 + theme_bw() + ylab("ASVs")
p1
```

Using viridis color palette  
```{r core-example3_plot, warning=FALSE, eval=TRUE}

library(viridis)
print(p1 + scale_fill_viridis())

```


As it can be seen, we see only OTu IDs and this may not be useful to interpret the data. We need to repreoccess this figure to include taxonomic information. We can do this as follows:  

```{r,core-example4, warning=FALSE, eval=TRUE}
library(RColorBrewer)
library(knitr)
# get the data used for plotting 
df <- p1$data 

# get the list of OTUs
list <- df$Taxa 

# check the OTU ids
# print(list) 

# get the taxonomy data
tax <- as.data.frame(tax_table(ps.m3.rel))

# add the ASVs to last column
tax$ASV <- rownames(tax)

# select taxonomy of only 
# those OTUs that are used in the plot
tax2 <- dplyr::filter(tax, rownames(tax) %in% list) 

# head(tax2)

# We will merege all the column into one except the Doamin as all is bacteria in this case
tax.unit <- tidyr::unite(tax2, Taxa_level,c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species", "ASV"), sep = "_;", remove = TRUE)

tax.unit$Taxa_level <- gsub(pattern="[a-z]__",replacement="", tax.unit$Taxa_level)

# add this new information into the plot data df

df$Taxa <- tax.unit$Taxa_level

# you can see now we have the taxonomic information
knitr::kable(head(df))

# replace the data in the plot object
p1$data <- df

plot(p1 + theme(axis.text.y = element_text(face="italic")))
```


## Genus level 

```{r}
ps.m3.rel.gen <- aggregate_taxa(ps.m3.rel, "Genus")
```

```{r fig.width=8, warning=FALSE, eval=TRUE}
library(RColorBrewer)
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-4), log10(.2), length = 10)

p1 <- plot_core(ps.m3.rel.gen, 
                plot.type = "heatmap", 
                colours = rev(brewer.pal(5, "RdBu")),
                prevalences = prevalences, 
                detections = detections, min.prevalence = .5) +
    xlab("Detection Threshold (Relative Abundance (%))")
p1 <- p1 + theme_bw() + ylab("ASVs")
p1
```