ordination.plots.Rmd

---
title: "NMDS Plots Visualizing Community Divergence"
date: '2022.04.30'
author: 'Nathan Malamud'
---

```{r setup}
# markdown setup
library(knitr)
knitr::opts_chunk$set(echo = F) # Don't print code
knitr::opts_chunk$set(warning = F) # Don't print warnings
knitr::opts_chunk$set(message = F) # Don't print messages
knitr::opts_chunk$set(echo = TRUE)

# Other important libraries
library(phyloseq)
library(vegan)
library(dplyr)
library(reshape2)
library(EcolUtils)
library(spaa)
library(ggplot2)
library(ggpubr)
library(tidyverse)

# Clear all objects from workspace
rm(list = ls())

asvtab <- readRDS('./Data/asvtab.rds')
metadata <- readRDS('./Data/metadata.rds')
guilds <- readRDS('./Data/guilds.rds')
itsphyseq <- readRDS('./Data/itsphyseq.rds')

# subset itsphyseq and metadata across soil and litter
itsphyseq.soil = subset_samples(itsphyseq, Soil_Litter=="Soil")
itsphyseq.litter = subset_samples(itsphyseq, Soil_Litter=="Litter")

metadata.soil = subset(metadata, Soil_Litter=="Soil")
metadata.litter = subset(metadata, Soil_Litter=="Litter")

# subset asvs across functional guilds
guilds.amf <- subset(guilds, Guild=="Arbuscular Mycorrhizal")
asvs.amf <- guilds.amf$ASV_ID

guilds$groups[grepl("Plant Pathogen", guilds$Guild)] = "Plant Pathogen"
guilds.pathogen <- subset(guilds, groups=="Plant Pathogen")
asvs.pathogen <- guilds.pathogen$ASV_ID

guilds$groups[grepl("Saprotroph", guilds$Guild)] = "Saprotroph"
guilds.sap <- subset(guilds, groups=="Saprotroph")
asvs.sap <- guilds.sap$ASV_ID

# Function for removing plot legend
library(gridExtra)
g_legend<-function(a.gplot){
    tmp <- ggplot_gtable(ggplot_build(a.gplot))
    leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
    legend <- tmp$grobs[[leg]]
    legend
}
```

```{r soil and litter}
set.seed(111)
nmds <- ordinate(
  physeq = itsphyseq,
  method = "NMDS",
  formula = ~Soil_Litter,
  distance = "bray"
)

sl=plot_ordination(
  physeq = itsphyseq,
  ordination = nmds,
  color = "Soil_Litter",
  title = "All Fungi across Soil and Litter"
) + 
  aes(group=Soil_Litter)+
  geom_point(aes(color = Soil_Litter), size = 2) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),strip.text = element_text(colour = 'white'),text = element_text(size=12),legend.text=element_text(size=12),legend.title = element_blank())

plot(sl)

ggsave("Figures/thesis.fig3.tiff", plot=sl, width=5, height = 5, units="in", dpi=300, compression="lzw")
```

```{r tree plotting function}

plot.across.trees <- function (physeq.obj, plot.title = "") {
  # This is a function that returns a plot of ASVs across tree species
  # Uses NMDS to generate ordination plots.
  #
  # Parameters:
  #   physeq_obj, a physeq object
  #   plot_title, a string of text
  #
  # Returns: an plot object.
  #
  
  set.seed(111)
  nmds <- ordinate(
    physeq = physeq.obj,
    method = "NMDS",
    formula = ~Tree_species,
    distance = "bray"
  )
  
  # Create dataframe of scores + covariates
  # Only want scores to return values for "sites"
  nmds.scores <- as.data.frame(scores(nmds, display="sites"))
  nmds.scores$Tree_species <- physeq.obj@sam_data$Tree_species
  physeq.obj@sam_data$Tree_species <- as.factor(physeq.obj@sam_data$Tree_species)
  
  # calculating means and sd's to get centroid + error bars
  means1 <- aggregate(nmds.scores$NMDS1,
                      by = list(nmds.scores$Tree_species),
                      FUN = mean)

  means2 <- aggregate(nmds.scores$NMDS2,
                      by = list(nmds.scores$Tree_species),
                      FUN = mean)

  sd1 <- aggregate(nmds.scores$NMDS1,
                      by = list(nmds.scores$Tree_species),
                      FUN = mean)

  sd2 <- aggregate(nmds.scores$NMDS2,
                      by = list(nmds.scores$Tree_species),
                      FUN = sd)

  centroids <- data.frame("means1" = means1$x, "means2" = means2$x,
                          "sd1" = sd1$x, "sd2" = sd2$x,
                          "Tree_species"=means1$Group.1)
  
  nmds.scores$Tree_species <- as.factor(nmds.scores$Tree_species)
  centroids$Tree_species <- as.factor(centroids$Tree_species)

  # plotting!
  pl = ggplot(data = nmds.scores, aes(x = NMDS1, y = NMDS2, color = Tree_species)) +
    geom_point(alpha = 0.3) + 
    geom_errorbar(data = centroids, aes(x = means1, y = means2, ymin = means2 - sd2, ymax = means2 + sd2)) +
    geom_errorbar(data = centroids, aes(x = means1, y = means2, xmin = means1 - sd1, xmax = means1 + sd1)) +
    geom_point(data = centroids, aes(x = means1, y = means2, color = Tree_species), size = 5) +
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),strip.text = element_text(colour = 'white'),text = element_text(size=12),legend.text=element_text(size=12),legend.title = element_blank())

  return(pl)
}

```

```{r amf fungi, echo=F}
itsphyseq.soil.amf <- subset_taxa(itsphyseq.soil, ASV_ID %in% asvs.amf)
itsphyseq.litter.amf <- subset_taxa(itsphyseq.litter, ASV_ID %in% asvs.amf)

# We need to prune out AMF samples that do not show up in soil or litter
amfsoilPS <- prune_samples(sample_sums(itsphyseq.soil.amf)>0, itsphyseq.soil.amf)
amflitterPS <- prune_samples(sample_sums(itsphyseq.litter.amf)>0, itsphyseq.litter.amf)

amf.soil.plt = plot.across.trees(amfsoilPS, "AMF in Soil")
amf.litter.plt = plot.across.trees(amflitterPS, "AMF in Litter")

plot(amf.soil.plt)
plot(amf.litter.plt)
```

```{r saprotrophic fungi}
itsphyseq.soil.sap <- subset_taxa(itsphyseq.soil, ASV_ID %in% asvs.sap)
itsphyseq.litter.sap <- subset_taxa(itsphyseq.litter, ASV_ID %in% asvs.sap)

sap.soil.plt = plot.across.trees(itsphyseq.soil.sap, "Decomposers in Soil")
sap.litter.plt = plot.across.trees(itsphyseq.litter.sap, "Decomposers in Litter")

plot(sap.soil.plt)
plot(sap.litter.plt)
```

```{r pathogenic fungi}
itsphyseq.soil.pathogen <- subset_taxa(itsphyseq.soil, ASV_ID %in% asvs.pathogen)
itsphyseq.litter.pathogen <- subset_taxa(itsphyseq.litter, ASV_ID %in% asvs.pathogen)

# Remove outlier due to low read count
pathogensoilPS <- subset_samples(itsphyseq.soil.pathogen, Sample_ID != "HITR_182HS")

pathogen.soil.plt = plot.across.trees(pathogensoilPS, "Pathogens in Soil")
pathogen.litter.plt = plot.across.trees(itsphyseq.litter.pathogen, "Pathogens in Litter")

plot(pathogen.soil.plt)
plot(pathogen.litter.plt)
```

```{r save all plots, include=False, results='hide', echo=F}
# Soil on the left, litter on the right

tree.legend = g_legend(amf.soil.plt)

amf.charts <- ggarrange(amf.soil.plt + theme(legend.position="none"), NULL,
                        amf.litter.plt + theme(legend.position="none"),
                        tree.legend,
                        nrow=1, ncol=4,
                        widths=c(1, 0.05, 1, 0.5),
                        labels=c("A", "", "B", ""))

sap.charts <- ggarrange(sap.soil.plt + theme(legend.position="none"), NULL,
                        sap.litter.plt + theme(legend.position="none"),
                        tree.legend,
                        nrow=1, ncol=4,
                        widths=c(1, 0.05, 1, 0.5),
                        labels=c("C", "", "D", ""))

pathogen.charts <- ggarrange(pathogen.soil.plt + theme(legend.position="none"), NULL,
                        pathogen.litter.plt + theme(legend.position="none"),
                        tree.legend,
                        nrow=1, ncol=4,
                        widths=c(1, 0.05, 1, 0.5),
                        labels=c("E", "", "F", ""))

ggsave(filename = "Figures/thesis.fig4AB.amf.tiff", plot=amf.charts, width=10, height = 5, units="in", dpi=300, compression="lzw")
ggsave(filename = "Figures/thesis.fig4CD.saps.tiff", plot=sap.charts, width=10, height = 5, units="in", dpi=300, compression="lzw")
ggsave(filename = "Figures/thesis.fig4EF.pathogens.tiff", plot=pathogen.charts, width=10, height = 5, units="in", dpi=300, compression="lzw")
```