4_eda_88soils.qmd

---
title: "Exploratory analysis of compositional data (part II)"
format: 
  html:
    code-fold: false
    code-tools: true
bibliography: references.bib
editor: visual
---

## Exploratory analysis of a microbial data set

Here we use the "88 soils" data set [@lauber2009pyrosequencing] containing bacterial communities in 88 soils from across North and South America.

### Load packages and data data

```{r}
#| label: load_data
#| warning: false
#| message: false
#| echo: true

library(phyloseq)
library(DESeq2)
library(vsn)
library(microbiome)
library(dplyr)
library(ggplot2)
library(gridExtra)

soil <- readRDS("data/soil_processed/soil.RDS")
```

### Exploration of the phyloseq object

```{r}
#| label: explore_phyloseq
#| warning: false
#| echo: true

soil

# Read count table
otutab <- otu_table(soil)
otutab[1:10, 1:6]

# Taxonomy table
taxtab <- tax_table(soil)
head(taxtab)

# Sample data
sampdata <- sample_data(soil)
head(sampdata)
```

### Add sample variables

We add a few sample variables we will need later in this tutorial:

* sampleID
* phType (acidic if pH<6.7 and basic if pH>6.7)
* totalReads (total number of reads per sample)

```{r}
#| label: add_sample_data
sampdata <- data.frame(sampleID = rownames(sampdata),
                       totalReads = sample_sums(soil),
                       ph = sampdata$ph)

sampdata$phType <- as.factor(ifelse(sampdata$ph < 6.7, "acid", "basic")) 

head(sampdata)
table(sampdata$phType)

# Add data frame to phyloseq object
sample_data(soil) <- sampdata
```

### Sequencing depth / library size

Plot the sequencing depth (total number of reads) of each sample.

```{r}
#| label: sequencing_depth

p <- ggplot(sampdata, aes(x = sampleID, y = totalReads)) +
  theme_bw() +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_blank()) +
  xlab("Sample")
p
```


```{r}
min(sampdata$totalReads)

torm <- sampdata[which.min(sampdata$totalReads), "sampleID"]
```

Remove the sample with only one read count.

```{r}
#| label: remove_sample

soil <- subset_samples(soil, sampleID != torm)
soil

otutab <- as(otu_table(soil), "matrix")
sampdata <- sample_data(soil)
```


### Sparsity

Number of zeros and percentage of zeros in the OTU table

```{r}
#| label: sparsity

nvar <-dim(otutab)[2]
nsamp <- dim(otutab)[1]

sum(otutab == 0)
sum(otutab == 0) / (nvar * nsamp) * 100
```

`r round(sum(otutab == 0) / (nvar * nsamp) * 100)` of the data are zeros.

### Visualize microbial composition

Have a look at the microbial compositions on phylum level.

```{r}
#| label: agglomerate_phylum
# Agglomerate to phylum level
soil_phyla <- tax_glom(soil, taxrank = "Phylum")
```


Stacked bar plot of the observed "absolute" abundances:

```{r, fig.width=8, fig.height=3}
#| label: barplot_absolute
#| fig-width: 8
#| fig-height: 3

p <- plot_bar(soil_phyla, fill = "Phylum")

p + theme_bw() +
  ylab("Absolute abundance") +
  theme(axis.text.x = element_blank()) +
  scale_fill_brewer(palette = "Set2")
```

Stacked bar plot of the relative abundances:

```{r}
# Compute relative abundances
soil_phyla_rel <- transform_sample_counts(soil_phyla, function(x) x/sum(x))
```

```{r, fig.width=8, fig.height=3}
#| label: barplot_relative
#| fig-width: 8
#| fig-height: 3
p <- plot_bar(soil_phyla_rel, fill = "Phylum")

p + theme_bw() +
  ylab("Absolute abundance") +
  theme(axis.text.x = element_blank()) +
  scale_fill_brewer(palette = "Set2")
```

### Zero replacement 

Zero counts are replaced by a unit pseudo count. In doing so, ratios between non-zero counts are preserved, which is not the case if a pseudo count is added to the whole matrix (which is also common).

```{r}
#| label: replace_zeros

soil_zrepl <- soil
otu_table(soil_zrepl)[otu_table(soil_zrepl) == 0] <- 1

otutab_zrepl <- as(otu_table(soil_zrepl), "matrix")
```

The `zCompositions` R package provides more complex methods for zero replacement.

### Mean - variance relationship

```{r, fig.width=8, fig.height=5}
#| label: mean-sd-plot_raw
#| fig-width: 8
#| fig-height: 5
vsn::meanSdPlot(t(otutab), plot = FALSE)$gg +
  theme_bw() + theme(text = element_text(size = 16))
```

### Normalization

We compare three normalization methods we have already seen in the lectures:

* CLR (centered log-ratio) transformation
* Variance stabilizing transformation (vst)
* Regularized log transformation (similar to vst but more robust when the size factors vary widely).

```{r}
#| label: normalization

# clr
soil_clr <- microbiome::transform(soil_zrepl, transform = "clr")
clr_counts <- otu_table(soil_clr)

# vst (function expects samples in columns)
vst_counts <- DESeq2::varianceStabilizingTransformation(t(otutab_zrepl), 
                                                         fitType = "local")

# rlog
rlog_counts <- DESeq2::rlog(t(otutab_zrepl), fitType = "local")
```

Create mean-sd plots using the function from `vsn` package.

```{r, fig.width=12, fig.height=4}
#| label: mean-sd-plots_norm
#| fig-width: 12
#| fig-height: 4
pclr <- meanSdPlot(t(clr_counts), plot = FALSE)$gg +
  theme_bw() + theme(text = element_text(size = 16)) + ggtitle("clr")
pvst <- meanSdPlot(vst_counts, plot = FALSE)$gg +
  theme_bw() + theme(text = element_text(size = 16)) + ggtitle("vst")
prlog <- meanSdPlot(rlog_counts, plot = FALSE)$gg +
  theme_bw() + theme(text = element_text(size = 16)) + ggtitle("rlog")

gridExtra::grid.arrange(pclr, pvst, prlog, 
             ncol = 3, nrow = 1)
```

### Alpha diversity

Another common task is diversity analysis. Alpha diversity summarizes the distribution of species abundances in a given sample.

```{r, fig.width=10, fig.height=4}
#| label: alpha_diversity_soil
#| fig-width: 10
#| fig-height: 4

pAlpha <- plot_richness(soil,
                       color = "phType",
                       measures = c("Observed", "Shannon", "InvSimpson", "Chao1"),
                       title = "Alpha diveristy for 88 soil data")

pAlpha + geom_point(size = 2) + theme_bw() + theme(axis.text.x = element_blank())
```

### Low-dimensional representation

In this section, we want to plot the data in the two-dimensional space. We therefore compute distances between all samples and perform multi-dimensional scaling (MDS). Here, we use the non-normalized data (see [@mcKnight2019methods])

Distance functions provided by the phyloseq package:

```{r}
#| label: distances
dist_methods <- unlist(distanceMethodList)
dist_methods
```

#### MDS plot with Bray Curtis distance

```{r}
#| label: mds_bray

soil_mds_bray <- ordinate(soil, method = "MDS", distance = "bray")
```

First a plot without coloring metadata.

```{r, fig.width=8, fig.height=5}
#| label: MDS_plot_bray
#| fig-width: 8
#| fig-height: 5
mds_soil_bray <- plot_ordination(soil, soil_mds_bray, 
                                 title = "MDS of 88 soil data")

mds_soil_bray + 
  theme_bw() + 
  theme(text = element_text(size = 14)) + 
  geom_point(size = 3)
```

One could identify two main clusters, separated by the first axis. 

Add the new data frame to the phyloseq object.

```{r}
sample_data(soil) <- sampdata
```


```{r, fig.width=8, fig.height=5}
#| label: MDS_plot_bray_labeled
#| fig-width: 8
#| fig-height: 5
mds_soil_bray <- plot_ordination(soil, soil_mds_bray, color = "phType",
                                 title = "")

mds_soil_bray + 
  theme_bw() + 
  theme(text = element_text(size = 14)) + 
  geom_point(size = 3) +
  ggtitle("MDS of 88 soil data colored by pH type")
```

Indeed, pH type separates samples along Axis.1. 

#### MDS plot with Aitchison distance

We generate another MDS plot, but this time using the Aitchison distance, which is simply the Euclidean distance of clr-transformed counts.

```{r}
#| label: mds_aitchison

soil_mds_ait <- ordinate(soil_clr, method = "MDS", distance = "euclidean")
```

```{r, fig.width=8, fig.height=4}
#| label: MDS_plot_bray_aitchison
#| fig-width: 8
#| fig-height: 4

sample_data(soil_clr) <- sampdata

mds_soil_ait <- plot_ordination(soil_clr, soil_mds_ait, color = "phType",
                                title = "")


p1 <- mds_soil_bray + 
  theme_bw() + 
  theme(text = element_text(size = 12)) + 
  geom_point(size = 3) +
  ggtitle("Bray Curtis dissimilarity")

p2 <- mds_soil_ait + 
  theme_bw() + 
  theme(text = element_text(size = 12)) + 
  geom_point(size = 3) +
  ggtitle("Aitchison distance")

gridExtra::grid.arrange(p1, p2,  
             ncol = 2, nrow = 1)
```


#### MDS plots for several distance measures

Now, we apply all available distance functions to the data, compute the corresponding MDS embeddings, and plot them.

```{r}
# Remove some unwanted distances
sel_dist <- dist_methods[!dist_methods %in% c("unifrac", "wunifrac",
                                              "dpcoa", "ANY", "z")]
```

Now we can apply these distance metrics to our data. 

```{r}
#| label: MDS_multiple_dist
gather_ordination_data <- function(d_method, phylo_obj){
  # Calculate the MDS matrix using the distance
  ordinate_obj <- ordinate(phylo_obj, method = "MDS", distance = d_method)
  # Get the data for plotting
  plot_data <- plot_ordination(phylo_obj, ordinate_obj)$data
  # Add a column for distance
  plot_data$distance <- d_method
  return(plot_data)
}

mds_plot_data <- bind_rows(lapply(sel_dist, gather_ordination_data,
                                  phylo_obj = soil))

head(mds_plot_data)
```

We color again according to pH type.

```{r, fig.width=14, fig.height=9}
#| label: MDS_plots_multiple_dist
#| fig-width: 14
#| fig-height: 9

mds_soil <- ggplot(mds_plot_data, aes(x=Axis.1, y=Axis.2, color=phType)) +
  geom_point() +
  facet_wrap(~ distance, scales = "free") +
  labs(x="Axis.1",
       y="Axis.2",
       title="Separation of 88 soil samples according to different distances")

mds_soil + theme_bw() + theme(text = element_text(size = 12))
```