nMDS.R

library(vegan)
library(rioja)
library(readxl)
library(writexl)
library(funrar)
library(tidyverse)
library(analogue)
library(janitor)
library(svglite)
library(viridis)

#read dataframe
df <- read_excel("input/ATP23-012_Analiza.xlsx", sheet = "Analiza-merged")

#metadata
meta <- read_excel("input/ATP23-012_Analiza.xlsx", sheet = "Metadata")
ages <- read_excel("input/Dating_biomarkers.xlsx")
metadata <- merge(meta, ages) %>%
  select(DNA_ID, Age, Site)

#taxonomy dataframe
taxo <- df %>%
  separate(Taxon, c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")) %>%
  select(1:8)

#only numbers about abundance
matrix <- df %>% select(1, 3:51) %>%
  column_to_rownames(var="#ASV ID")

#relative abundance of matrix
matrix_t <- t(matrix)
rel_abundance <- make_relative(matrix_t)
rel_abundance <- t(rel_abundance)
rel_abundance <- data.frame(rel_abundance) %>%
  mutate(Total = rowSums(.)) %>%
  rownames_to_column("#ASV ID") %>%
  inner_join(taxo)

rel_abundance <- data.frame(t(rel_abundance)) %>%
  row_to_names(row_number = 1) %>%
  rownames_to_column("DNA_ID") %>%
  merge(metadata) 
rel_abundance[2:3724] <- sapply(rel_abundance[2:3724],as.numeric)
rel_abundance <- rel_abundance %>%
  arrange(Age) %>%
  select(where(~ any(. != 0)))

#make community matrix - extract columns with abundance information
com = rel_abundance[,2:2548]

#turn abundance data frame into a matrix
m_com = as.matrix(com)

#calculate nMDS
set.seed(123)
nmds = metaMDS(m_com, distance = "bray")
nmds

#extract NMDS scores (x and y coordinates)
data.scores = as.data.frame(scores(nmds)$sites)
#add columns to data frame 
data.scores$Site = rel_abundance$Site
data.scores$Age = rel_abundance$Age

nmds_plot = ggplot(data.scores, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(size = 6, aes(colour = Age, shape = Site))+ 
  theme(axis.text.y = element_text(colour = "black", size = 12, face = "bold"), 
        axis.text.x = element_text(colour = "black", face = "bold", size = 12), 
        legend.text = element_text(size = 12, face ="bold", colour ="black"), 
        legend.position = "right", axis.title.y = element_text(face = "bold", size = 14), 
        axis.title.x = element_text(face = "bold", size = 14, colour = "black"), 
        legend.title = element_text(size = 14, colour = "black", face = "bold"), 
        panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA, linewidth = 1.2),
        legend.key=element_blank()) + 
  labs(x = "NMDS1", colour = "Age", y = "NMDS2", shape = "Site") +
  scale_color_viridis(discrete = FALSE) +
  geom_text(aes(label = Age), size = 4, vjust = 1.5)
nmds_plot
stressplot(nmds)
ggsave("nMDS_CHARME.png", plot = nmds_plot , bg = "white", width = 12, height = 8, device = "png", path = "output")