# Jan 20, 2025
# R version 4.4.0
# Script to generate nMDS plot of MOB 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)

# remove singleton ASVs (only one count)
asv_filt <- asv_table[rowSums(asv_table)>1,]   
dim(asv_filt)  # 26104 ASVs remain

# load metadata 
metadata <- read.csv(file='filtered_metadata.csv')

# List of MOB ASVs (excl. Methylomonadaceae-NA)
mob_asvs <- c(
  "asv.47", "asv.48", "asv.55", "asv.81", "asv.106", "asv.125",
  "asv.146", "asv.160", "asv.188", "asv.201", "asv.226", "asv.244",
  "asv.280", "asv.283", "asv.294", "asv.322", "asv.383", "asv.769",
  "asv.896", "asv.921", "asv.924", "asv.979", "asv.1054", "asv.1285",
  "asv.1296", "asv.1320", "asv.1344", "asv.1351", "asv.1356", "asv.1394",
  "asv.1419", "asv.1520", "asv.1720", "asv.1826", "asv.1833", "asv.1840",
  "asv.1873", "asv.1917", "asv.1980", "asv.2010", "asv.2143", "asv.2835",
  "asv.2873", "asv.2915", "asv.2951", "asv.3704", "asv.3791", "asv.4036",
  "asv.4222", "asv.5187", "asv.5351", "asv.5542", "asv.6147", "asv.6651",
  "asv.6982", "asv.7206", "asv.7243", "asv.7844", "asv.8207", "asv.8501",
  "asv.10009", "asv.10023", "asv.10489", "asv.11175", "asv.11761", "asv.12005",
  "asv.12108", "asv.12111", "asv.12216", "asv.12384", "asv.12483", "asv.12923",
  "asv.13051", "asv.13386", "asv.13554", "asv.13822", "asv.14104", "asv.14153",
  "asv.14480", "asv.14618", "asv.15609", "asv.15631", "asv.16445", "asv.16933",
  "asv.17396", "asv.17641", "asv.18453", "asv.18501", "asv.19722", "asv.20575",
  "asv.20917", "asv.20930", "asv.20939", "asv.21344", "asv.21535", 
  "asv.21784", "asv.21847", "asv.22392", "asv.22395", "asv.22497", "asv.22965",
  "asv.23000", "asv.23008", "asv.23652", "asv.23656", "asv.24385", "asv.24397",
  "asv.24450", "asv.24518", "asv.24545", "asv.25236", "asv.25342", "asv.25447",
  "asv.25620", "asv.26323", "asv.26354", "asv.26368", "asv.27693", "asv.27814",
  "asv.29157"
)

####################
# Use double Wisconsin normalization on sqrt-counts (all non-singleton counts)
asv.SqRt <- sqrt(asv_filt) # SqRt transformation (omit ASV_ID col)
asv.Wis <- wisconsin(t(asv.SqRt)) # transpose bc function expects samples in rows
asv.Wis <- t(asv.Wis)  # transpose back
# add back ASV IDs as rownames
rownames(asv.Wis) <- rownames(asv_filt)  

# reduce asvs to MOB asvs (120), remove samples that have no MOB counts (all zero)
asv.Wis.mob <- asv.Wis[mob_asvs,]
idx <- which(colSums(asv.Wis.mob)> 0) 
length(idx)  # 291 samples remain
asv.Wis.mob <- asv.Wis.mob[, idx]

# same as with all asvs, require non-zero count in at least 5% of the samples
# to keep an ASV
zeroes = function(vec, percentage, n) {
  sum(vec==0) < n*percentage  
}

# only 17 ASVs remain
filt5 <- asv.Wis.mob[apply(asv.Wis.mob,1,zeroes, percentage=0.95, n=291),]   # keep ASVs present in at least 5% of the samples, 276 zeros are ok
dim(filt5)  

# remove again samples that have no counts in these ASVs
toRM <- which(colSums(filt5)==0)  # 8 samples
filt5 <- filt5[, -toRM]
dim(filt5) # 17 283


# calculate Bray-Curtis distances
mob_bc_wis_filt <- vegdist(t(filt5), method = "bray")
# Perform NMDS on the cleaned MOB ASV data
mob_nmds_wis_filt <- metaMDS(mob_bc_wis_filt, k = 2, trymax = 100)

# Extract NMDS scores for MOBs 
mob_nmds_scores <- as.data.frame(scores(mob_nmds_wis_filt))
mob_nmds_scores$sampleID <- colnames(filt5)


####################################
# Merge with metadata for plotting
mob_plot_data <- mob_nmds_scores %>%
  left_join(metadata, by = "sampleID") %>%
  mutate(treatment = paste(location, type, sep = " "))

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

# Plot the MOB NMDS
ggplot(mob_plot_data, aes(x = NMDS1, y = NMDS2, color = season, shape = treatment)) +
  geom_point(size = 4) +
  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 MOB Subset",
    x = "NMDS1",
    y = "NMDS2",
    color = "Season",
    shape = "Treatment"
  ) +
  theme(
    legend.position = "right",
    text = element_text(size = 12)
  )

mob_nmds_wis_filt$stress  # 0.09259482

# save plot to file
dev.print(device=pdf, file='Plots/nmds_MOBs_Filt_BC_stress0.093_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