Core.Rmd

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

See also related functions for the analysis of rare and variable taxa (rare_members; rare_abundance; rare_members; rare_abundance; low_abundance).  


```{r setup, message = FALSE, warning = FALSE, results = 'hide'}
library("devtools")
#install_github("microbiome/microbiome")
```



# HITChip Data

Load example data:


```{r core-prevalence, warning=FALSE, message=FALSE}
# Load data
library(microbiome)
data(peerj32)

# Rename the data
pseq <- peerj32$phyloseq

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


## Prevalence of taxonomic groups

Relative population frequencies; at 1% compositional abundance threshold:

```{r core-prevalence2}
head(prevalence(pseq.rel, detection = 1/100, sort = TRUE))
```


Absolute population frequencies (sample count):

```{r core-prevalence2b}
head(prevalence(pseq.rel, detection = 1/100, sort = TRUE, count = TRUE))
```


## 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(pseq.rel, detection = 0, prevalence = 50/100)
```


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

```{r core-data, message=FALSE, warning=FALSE}
pseq.core <- core(pseq.rel, detection = 0, prevalence = .5)
```

We can also collapse the rare taxa into an "Other" category

```{r core_aggregate_rare, message=FALSE, warning=FALSE}
pseq.core2 <- aggregate_rare(pseq.rel, "Genus", detection = 0, prevalence = .5)
```

Retrieving the core taxa names from the phyloseq object:

```{r core-taxa, message=FALSE, warning=FALSE}
core.taxa <- taxa(pseq.core)
```



## Core abundance and diversity

Total core abundance in each sample (sum of abundances of the core members):

```{r core-ab, message=FALSE, warning=FALSE}
core.abundance <- sample_sums(core(pseq.rel, detection = .01, prevalence = .95))
```


## 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, fig.width=9, fig.heigth=6, out.width="400px", warning=FALSE}
# With compositional (relative) abundances
det <- c(0, 0.1, 0.5, 2, 5, 20)/100
prevalences <- seq(.05, 1, .05)
 #ggplot(d) + geom_point(aes(x, y)) + scale_x_continuous(trans="log10", limits=c(NA,1))


plot_core(pseq.rel, prevalences = prevalences, detections = det, plot.type = "lineplot") + xlab("Relative Abundance (%)")
```


## 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 taxa.order argument.

```{r core-example3, fig.width=8, fig.heigth=18, fig.show='hold', warning=FALSE, out.width = '200px', warning=FALSE}
# Core with compositionals:

library(RColorBrewer)
library(reshape)
prevalences <- seq(.05, 1, .05)
detections <- 10^seq(log10(1e-3), log10(.2), length = 10)

# Also define gray color palette
gray <- gray(seq(0,1,length=5))
#pseq.rel<- microbiome::transform(pseq, 'compositional')
p <- plot_core(pseq.rel, plot.type = "heatmap", colours = gray,
       prevalences = prevalences, detections = detections) +
    labs(x = "Detection Threshold (Relative Abundance (%))")
print(p)    

# Core with absolute counts and horizontal view:
# and minimum population prevalence (given as percentage)
detections <- 10^seq(log10(1), log10(max(abundances(pseq))/10), length = 10)

p <- plot_core(pseq, plot.type = "heatmap", 
       		 prevalences = prevalences,
       		 detections = detections,
		 colours = rev(brewer.pal(5, "Spectral")),
		 min.prevalence = .2, horizontal = TRUE) 
print(p)
```

# Core Microbiota using Amplicon data

## 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 [(Halfvarson et al. Nature Microbiology 2, 2017)](http://www.nature.com/articles/nmicrobiol20174) will be used and only healthy samples will be included.  

```{r, core-microbiota-amplicon-data, eval=FALSE}
# check the data 
print(pseq.biom) 

# Filter the data to include only healthy subjects
pseq.1 <- subset_samples(pseq.biom, ibd_subtype == "HC" & timepoint == "1") 
print(pseq.1)
# keep only taxa with positive sums
pseq.2 <- prune_taxa(taxa_sums(pseq.1) > 0, pseq.1)

print(pseq.2)

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


### Prevalence of taxonomic groups

Relative population frequencies; at 1% compositional abundance threshold:

```{r core-prevalenceA, eval=FALSE}
head(prevalence(pseq.rel, detection = 1, sort = TRUE))
```

We can see that only OTU ids are listed with no taxonomic information.
Absolute population frequencies (sample count):

```{r core-prevalenceB, eval=FALSE}
head(prevalence(pseq.rel, detection = 1, sort = TRUE, count = TRUE))
```


### 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-membersA, message=FALSE, warning=FALSE, eval = FALSE}
core.taxa.standard <- core_members(pseq.rel, detection = 0, prevalence = 50/100)
```


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

```{r core-dataA, message=FALSE, warning=FALSE}
pseq.core <- core(pseq.rel, detection = 0, prevalence = .5)
```


Retrieving the associated taxa names from the phyloseq object:

```{r core-taxaB, message=FALSE, warning=FALSE, eval=FALSE}
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 abundance and diversity

Total core abundance in each sample (sum of abundances of the core members):

```{r core-abc, message=FALSE, warning=FALSE, eval=FALSE}
core.abundance <- sample_sums(core(pseq.rel, detection = .001, prevalence = .95))
```


## 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 coreC, fig.width=9, fig.heigth=6, out.width="400px", warning=FALSE, eval=FALSE}
# With compositional (relative) abundances
det <- c(0, 0.1, 0.5, 2, 5, 20)/100
prevalences <- seq(.05, 1, .05)
plot_core(pseq.rel, prevalences = prevalences, detections = det, plot.type = "lineplot") + xlab("Relative Abundance (%)")
```


### 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-exampleD, fig.width=8, fig.heigth=18, fig.show='hold', out.width = '200px', warning=FALSE,message=FALSE}
# Core with compositionals:
prevalences <- seq(.05, 1, .05)
#detections <- 10^seq(log10(1e-5), log10(.2), length = 10)
detections <- 10^seq(-5, 0, 1)

# Also define gray color palette
gray <- gray(seq(0,1,length=5))
p <- plot_core(pseq.rel, plot.type = "heatmap", colours = gray,
    prevalences = prevalences, detections = detections, min.prevalence = .5) +
    labs(x = "Detection Threshold (Relative Abundance (%))")
print(p)    


# Same with the viridis color palette
# color-blind friendly and uniform
# options: viridis, magma, plasma, inferno
# https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html
# Also discrete=TRUE versions available
library(viridis)
print(p + scale_fill_viridis())
```