# Jan 20, 2025
# R version 4.4.0
# Script to generate nMDS plot of whole bacterial communities of 
# environmental water column, sediment and incubation samples.

# Clear the R workspace
rm(list = ls())

set.seed(9876543) # as this may influence nmds calculations

# Install required packages if not already installed
if (!require("vegan")) install.packages("vegan")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("dplyr")) install.packages("dplyr")
if (!require("readr")) install.packages("readr")
if (!require("tidyr")) install.packages("tidyr")

# Load libraries
library(vegan)
library(ggplot2)
library(dplyr)
library(readr)
library(tidyr)  


# Load the ASV table 
asv_table <- read.csv(file='asv_table.csv', row.names=1)
dim(asv_table) 
# Rows/ASVs: 27012 Columns/samples: 318

# read sample metadata (for 318 samples)
metadata <- read.csv(file= 'filtered_metadata.csv')

# rm singletons 
idx <-  which(rowSums(asv_table)>1)   
length(idx)  # 26104
asv_filt <- asv_table[idx,]

asv.SqRt <- sqrt(asv_filt) # SqRt transformation
asv.Wis <- wisconsin(t(asv.SqRt)) # transpose bc function expects samples in rows
asv.Wis <- t(asv.Wis)  # transpose back

# remove ASVs with many zeros (at least 5% of the samples must be non-zero )
zeroes = function(vec, percentage, n) {
  sum(vec==0) < n*percentage  
}

filt.5 <- asv.Wis[apply(asv.Wis,1,zeroes, percentage=0.95, n=318),]  
dim(filt.5)  # 5266  ASVs, 318 samples

# compute Bray-Curtis dissimilarity matrix
# vegdist will cluster rows, so transpose asvs!
bray_curtis <- vegdist(t(filt.5), method = "bray")   

# Perform NMDS. 
# set autotransform=FALSE, otherwise a data maximum 50 triggers sqrt and values >9 triggers wisconsin double standardization!!
# (should not be triggered but to be on the safe side)
nmds_result <- metaMDS(bray_curtis, k = 2, trymax = 100, autotransform=FALSE)

# Extract NMDS scores
nmds_scores <- as.data.frame(scores(nmds_result))
# Map correct sampleID from transposed_data or merged_data 
nmds_scores$sampleID <- colnames(filt.5)  # Use sampleID from filtered asv data

# PLOT
plot_data <- nmds_scores %>% 
  left_join(metadata, by = "sampleID")

# Combine location and type into a single variable for shape
plot_data <- plot_data %>% 
  mutate(treatment = paste(location, type, sep = " "))

# Identify T0 samples
t0_samples <- plot_data %>%
  filter(grepl("^T0", description))  # Select rows where description starts with "T0"

# Create the NMDS plot with an additional layer for T0 samples
ggplot(plot_data, aes(x = NMDS1, y = NMDS2, color = season, shape = treatment)) +
  geom_point(size = 4) +  # Main points
  geom_point(data = t0_samples, aes(x = NMDS1, y = NMDS2), shape = 4, size = 1, color = "black", stroke=1.5) +  # Cross for T0
  scale_shape_manual(
    values = c(
      "north sea incubation" = 15,    # Filled square
      "north sea watercolumn" = 17, # Filled triangle
      "north sea sediment" = 16,     # Filled circle
      "wadden sea incubation" = 0,   # Open square
      "wadden sea watercolumn" = 2, # Open triangle
      "wadden sea sediments" = 1      # Open circle
    )
  ) +
  scale_color_manual(
    values = c(
      "winter" = "blue",
      "spring" = "green",
      "summer" = "red",
      "autumn" = "brown"
    )
  ) +
  theme_minimal() +
  labs(
    title = "NMDS of Microbial Community with T0 Samples Highlighted",
    x = "NMDS1",
    y = "NMDS2",
    color = "Season",
    shape = "Treatment"
  ) +
  theme(
    legend.position = "right",
    text = element_text(size = 12)
  )

# Display the stress level of the NMDS result
nmds_result$stress

# save plot
dev.print(device=pdf, file='Plots/nmds_5percFiltered_BC_stress0.16_20250120.pdf')


############################################
sessionInfo()

# R version 4.4.0 (2024-04-24 ucrt)
#Platform: x86_64-w64-mingw32/x64
#Running under: Windows 11 x64 (build 22621)

#Matrix products: default

#locale:
#  [1] LC_COLLATE=English_Europe.utf8  LC_CTYPE=English_Europe.utf8   
#[3] LC_MONETARY=English_Europe.utf8 LC_NUMERIC=C                   
#[5] LC_TIME=English_Europe.utf8    

# time zone: Europe/Amsterdam
#tzcode source: internal

#attached base packages:
#  [1] stats     graphics  grDevices utils     datasets  methods   base     

#other attached packages:
#  [1] tidyr_1.3.1    readr_2.1.5    dplyr_1.1.4    ggplot2_3.5.1  vegan_2.6-6.1 
#[6] lattice_0.22-6 permute_0.9-7 

#loaded via a namespace (and not attached):
#  [1] bit_4.0.5            Matrix_1.7-0         gtable_0.3.5        
#[4] crayon_1.5.2         compiler_4.4.0       tidyselect_1.2.1    
#[7] Rcpp_1.0.12          colorBlindness_0.1.9 gridGraphics_0.5-1  
#[10] parallel_4.4.0       cluster_2.1.6        splines_4.4.0       
#[13] scales_1.3.0         R6_2.5.1             plyr_1.8.9          
#[16] labeling_0.4.3       generics_0.1.3       MASS_7.3-60.2       
#[19] tibble_3.2.1         munsell_0.5.1        tzdb_0.4.0          
#[22] pillar_1.9.0         rlang_1.1.4          utf8_1.2.4          
#[25] bit64_4.0.5          pkgload_1.4.0        cli_3.6.2           
#[28] withr_3.0.0          magrittr_2.0.3       mgcv_1.9-1          
#[31] grid_4.4.0           vroom_1.6.5          rstudioapi_0.16.0   
#[34] hms_1.1.3            cowplot_1.1.3        lifecycle_1.0.4     
#[37] nlme_3.1-164         vctrs_0.6.5          glue_1.7.0          
#[40] farver_2.1.2         fansi_1.0.6          colorspace_2.1-0    
#[43] purrr_1.0.2          tools_4.4.0          pkgconfig_2.0.3