Landscaping.Rmd

---
title: "Microbiome Landscapes"
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 - density}
  %\usepackage[utf8]{inputenc}
  %\VignetteEncoding{UTF-8}  
-->


[Microbiome Landscaping](https://academic.oup.com/femsre/article/doi/10.1093/femsre/fuw045/2979411/Intestinal-microbiome-landscaping-insight-in#58802539) refers to the analysis and illustration of population frequencies. Typically, these are wrappers based on standard ordination methods (for more examples, see [ordination examples](Ordination.html))


## Two-dimensional microbiome landscape

Load example data:

```{r landscaping, message=FALSE, warning=FALSE, eval=TRUE}
library(microbiome)
library(phyloseq)
library(ggplot2)

data(dietswap)
pseq <- dietswap

# Convert to compositional data
pseq.rel <- microbiome::transform(pseq, "compositional")

# Pick core taxa
pseq.core <- core(pseq.rel, detection = 5/100, prevalence = 5/100)
pseq.core <- subset_samples(pseq.core, sex == "female" &
	                               bmi_group == "overweight")
```


## Landscape figure

Visualize the microbiome landscape (sample similarities on two-dimensional projection). When using these tools, kindly cite Shetty et al. FEMS Microbiology Reviews, 41(2):182–199, 2017 [doi:10.1093/femsre/fuw045](https://doi.org/10.1093/femsre/fuw045).


### PCA

```{r landscape_pca, message=FALSE, warning=FALSE, fig.width=8, fig.height=6, fig.show="hold", out.width="400px"}
# PCA with euclidean distance and CLR transformation
p <- plot_landscape(pseq, method = "PCA", transformation = "clr") +
       labs(title = paste("PCA / CLR"))
print(p)
```


### PCoA / MDS

```{r landscape_pcoa, message=FALSE, warning=FALSE, fig.width=8, fig.height=6, fig.show="hold", out.width="400px", eval=FALSE}
# PCoA for compositional data with Bray-Curtis distances
p <- plot_landscape(microbiome::transform(pseq.core, "compositional"),
                      method = "PCoA", distance = "bray") +
       labs(title = paste("PCoA / Compositional / Bray-Curtis"))
print(p)
```


### t-SNE

```{r landscape_tsne, message=FALSE, warning=FALSE, fig.width=8, fig.height=6, fig.show="hold", out.width="400px"}
p <- plot_landscape(pseq, "t-SNE",
       distance = "euclidean", transformation = "hellinger") +
       labs(title = paste("t-SNE / Hellinger / Euclidean"))       
print(p)
```


### NMDS

```{r landscape3, message=FALSE, warning=FALSE, fig.width=8, fig.height=6, fig.show="hold", out.width="400px"}
# Landscape plot directly from phyloseq object
p <- plot_landscape(pseq.core, "NMDS", "bray", col = "nationality") +
       labs(title = paste("NMDS / Bray-Curtis"))       
```





For direct access to the ordination coordinates, use the following:

```{r landscape4, message=FALSE, warning=FALSE, fig.width=8, fig.height=6, fig.show="hold", out.width="400px"}
# Project the samples with the given method and dissimilarity measure. 
# Ordinate the data; note that some ordinations are sensitive to random seed
# "quiet" is used to suppress intermediate outputs
set.seed(423542)
x <- pseq.core
quiet(x.ord <- ordinate(x, "NMDS", "bray"))
# Pick the projected data (first two columns + metadata)
proj <- phyloseq::plot_ordination(x, x.ord, justDF=TRUE)
# Rename the projection axes
names(proj)[1:2] <- paste("Comp", 1:2, sep=".")

# Same with a generic data.frame
# (note that random seed will affect the exact ordination)
p <- plot_landscape(proj[, 1:2], col = proj$nationality, legend = T)
print(p)

# Visualize sample names:
ax1 <- names(proj)[[1]]
ax2 <- names(proj)[[2]]
p <- ggplot(aes_string(x = ax1, y = ax2, label = "sample"), data = proj) +
       geom_text(size = 2)
print(p)
```


## Abundance histograms (one-dimensional landscapes)

Population densities for Dialister:

```{r hist, fig.width=6, fig.width=8, fig.height=6, fig.show="hold", out.width="400px"}
# Load libraries
library(microbiome)
library(phyloseq)
pseq <- dietswap

# Visualize population densities for specific taxa
plot_density(pseq, "Dialister") + ggtitle("Absolute abundance")

# Same with log10 compositional abundances
x <- microbiome::transform(pseq, "compositional")
tax <- "Dialister"
plot_density(x, tax, log10 = TRUE) +
  ggtitle("Relative abundance") +
  xlab("Relative abundance (%)")
```