#####################################################
# Marco FW Gauger
# 2021/02/23
#
# This is the script to obtain fig 5 of the manuscript submitted to PeerJ Gauger et al 2021 Diel influences on bottlenose dolphin acoustic detection in a coastal lagoon in the southwestern Gulf of California
#
#
# personel computer Intel Core I7-7500 U, @2.6 GHz, RAM 8 GB, Windows 10 
# start 09:27
# end 09:27
#####################################################
ifelse("geosphere" %in% rownames(installed.packages()) == FALSE, 
       install.packages("geosphere", dependencies = T),
       ("geosphere is already installed"))
ifelse("dplyr" %in% rownames(installed.packages()) == FALSE, 
       install.packages("dplyr", dependencies = T),
       ("dplyr is already installed"))
ifelse("factoextra" %in% rownames(installed.packages()) == FALSE, 
       install.packages("factoextra", dependencies = T),
       ("factoextra is already installed"))
ifelse("FactoMineR" %in% rownames(installed.packages()) == FALSE, 
       install.packages("FactoMineR", dependencies = T),
       ("FactoMineR is already installed"))
ifelse("ggsci" %in% rownames(installed.packages()) == FALSE, 
       install.packages("ggsci", dependencies = T),
       ("ggsci is already installed"))
ifelse("cowplot" %in% rownames(installed.packages()) == FALSE, 
       install.packages("cowplot", dependencies = T),
       ("cowplot is already installed"))
ifelse("vegan" %in% rownames(installed.packages()) == FALSE, 
       install.packages("vegan", dependencies = T),
       ("vegan is already installed"))
ifelse("RColorBrewer" %in% rownames(installed.packages()) == FALSE, 
       install.packages("RColorBrewer", dependencies = T),
       ("RColorBrewer is already installed"))
ifelse("gplots" %in% rownames(installed.packages()) == FALSE, 
       install.packages("gplots", dependencies = T),
       ("gplots is already installed"))
ifelse("viridis" %in% rownames(installed.packages()) == FALSE, 
       install.packages("viridis", dependencies = T),
       ("viridis is already installed"))
ifelse("parameters" %in% rownames(installed.packages()) == FALSE, 
       install.packages(c("parameters","bayestestR")),
       ("parameters is already installed"))



require(gplots)
require(RColorBrewer)
require(vegan)
require(cowplot)
require(ggsci)
require(FactoMineR)
require(factoextra)
require(geosphere)
require(dplyr)



analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
names(analysed_data_hour)
analysed_data_hour$month_campaign[analysed_data_hour$month_campaign==c("11_13")] <- "11_12"
j <- aggregate(analysed_data_hour$dp10m_ELAP1, list(analysed_data_hour$month_campaign), sum)
names(j) <- c("campaign","dph")
j$d10pm <- 
  aggregate(analysed_data_hour$dp10m_ELAP1, list(analysed_data_hour$month_campaign), sum)$x
j$moon <- 
  aggregate(analysed_data_hour$lunar_phase , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$SST <- 
  aggregate(analysed_data_hour$SST , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$SST_BLAP <- 
  aggregate(analysed_data_hour$SST_BLAP , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$CHL <- 
  aggregate(analysed_data_hour$CHL_BLAP , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$depth <- 
  aggregate(analysed_data_hour$depth , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$sediment <- 
  aggregate(analysed_data_hour$sediment , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$sediment <- as.numeric(j$sediment)
j$distance_mangrove <- 
  aggregate(analysed_data_hour$distance_mangrove , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$distance_coast <- 
  aggregate(analysed_data_hour$distance_coast , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$distance_BLAP <- 
  aggregate(analysed_data_hour$distance_BLAP , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$delta_SST <- j$SST - j$SST_BLAP
j$daylength <- geosphere::daylength(24, as.Date(aggregate(analysed_data_hour$datehour , list(analysed_data_hour$month_campaign),  dplyr::first)$x, format="%d.%m.%Y"))
j$date <- 
  aggregate(analysed_data_hour$date , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$station <- 
  aggregate(analysed_data_hour$Station , list(analysed_data_hour$month_campaign),  dplyr::first)$x
j$date[j$campaign == "07_04"] <- "13.07.2017"

normalize <- function(x) {
return ((x - min(x)) / (max(x) - min(x)))
}

pca_j <- j
for(k in 2:length(j)){pca_j[,k] <- normalize(j[,k])}
names(pca_j)


parameters::check_kmo(pca_j[,c(5,7,14)]) # SST, CHL & daylength

rownames(pca_j) <- paste0("St.",j$station,"_", j$date)

res.pca_j2 <- PCA(pca_j[,c(5,7,14)], ncp = 2, graph = FALSE)
# Compute hierarchical clustering on principal components
res.hcpc_j2 <- HCPC(res.pca_j2, graph = FALSE)

map1 <- fviz_dend(res.hcpc_j2, 
                  cex = 0.7,                     
                  # Label size
                  palette = c("#1B9E77","grey25","magenta"),               
                  # Color palette see ?ggpubr::ggpar
                  rect = TRUE, rect_fill = TRUE, 
                  # Add rectangle around groups
                  rect_border = c("#1B9E77","grey25","magenta"),           
                  # Rectangle color
                  labels_track_height = 0.95    ,
                  main = ""
                  # Augment the room for labels
)
map1 + scale_fill_manual(values=c("#1B9E77","grey25","magenta"))
(map2 <- fviz_cluster(res.hcpc_j2,
                     repel = TRUE,            # Avoid label overlapping
                     show.clust.cent = TRUE, # Show cluster centers
                     palette =c("#1B9E77","grey25","magenta"),         # Color palette see ?ggpubr                     ggtheme = theme_minimal(),
                     main = "",
                     labelsize = 9)+
  theme_classic() )
map2 + guides(color = FALSE )+ 
  scale_color_manual(values=c("#1B9E77","grey25","magenta"))


map2 <- map2 + guides(color = FALSE )



plot_grid(map1+ylim(-.4,1.25), map2, labels = c('A', 'B'), label_size = 9)

ggsave("Supplemental Figure S3.png",
       plot_grid(map1+ylim(-.4,1.25), map2, 
                 labels = c('A', 'B'), label_size = 9), 
       dpi=300, units="cm", width=15, height=15)



dta = res.hcpc_j2$data.clust
str(dta)
dta$clust <-as.numeric(dta$clust)
dta$cluster <- dta$clust
dta$cluster[dta$clust == 3] <- "May-Jul"
dta$cluster[dta$clust == 1] <- "Aug-Oct"
dta$cluster[dta$clust == 2] <- "Oct-Mar"

ad = adonis(clust ~ SST+daylength+CHL, data= dta, method='eu')
summary(ad)
dis <- vegdist(dta[,1:3])

mod <- betadisper(dis, dta$cluster)
anova(mod)
pmod <- permutest(mod, permutations = 999, pairwise = TRUE)

(mod.HSD <- TukeyHSD(mod))
plot(mod.HSD, yaxt="n")
axis(side=2, las=2, at=1:3, 
     labels=c("May-\nJul","Aug-\nOct","Oct-\nMar"))


library(dendextend)
res.hcpc_j2$call$t$tree %>% 
  as.dendrogram() %>% 
  color_branches(k = 3, groupLabels = unique(res.hcpc_j2$data.clust$clust)) %>% 
  plot()

dend <- as.dendrogram((res.hcpc_j2$call$t$tree))
# order it the closest we can to the order of the observations:
dend <- rotate(dend, 1:21)

# Color the branches based on the clusters:
dend <- color_branches(dend, k=3) #, groupLabels=iris_species)
library(colorspace) 
# Manually match the labels, as much as possible, to the real classification of the flowers:
labels_colors(dend) <-c(rep(1,20),2)

# We hang the dendrogram a bit:
dend <- hang.dendrogram(dend,hang_height=0.1)
# reduce the size of the labels:
# dend <- assign_values_to_leaves_nodePar(dend, 0.5, "lab.cex")
dend <- set(dend, "labels_cex", 0.5)



#######
some_col_func <- function(n) 
  rev(colorspace::heat_hcl(n, c = c(80, 30),l = c(30, 90), power = c(1/5, 1.5)))


M_heatmap <- cbind(pca_j[,c(5,7,14)],aggregate(analysed_data_hour$dp10m_ELAP1, list(analysed_data_hour$month_campaign ), mean)$x*5)
names(M_heatmap) <- c("SST", "Chl","daylength","dp10m")
summary(M_heatmap)


M_heatmap$cluster <- res.hcpc_j2$data.clust$clust
M_heatmap <- M_heatmap[ order(M_heatmap$cluster), ]


Rowv  <- M_heatmap[1:3] %>% 
  scale %>% 
  dist %>% 
  hclust %>% 
  as.dendrogram %>% 
  set("branches_k_color", k = 3) %>% 
  set("branches_lwd", 1.2) %>%
  ladderize
Rowv <- color_branches(Rowv, k = 4, col = c(1,1,1,1))


#x11();par(oma=c(2,1,0,1))
#jpeg(file="heatmap_hcpc_scaled.jpeg",width=15, height=10, units="cm", res=150)

heatmap.2(as.matrix(scale(M_heatmap[1:4])), Rowv=Rowv, Colv=NA, tracecol=1,
          dendrogram = c("row"), col = grey.colors(500), 
          density.info = "none", rowsep = c(6,13), 
          colsep = c(3), srtCol=45, key.title="scaled", cexRow=.5, cexCol=.85,
          RowSideColors = c(    # grouping row-variables into different
            rep(rainbow(100)[92], 6),   # categories, Measurement 1-3: green
            rep(rainbow(100)[61], 8),    # Measurement 4-8: blue
            rep(rainbow(100)[32], 7)))
dev.off()

normalize <- function(x) {
return ((x - min(x)) / (max(x) - min(x)))
}
M_heatmap$SST <-        normalize(M_heatmap[1])
M_heatmap$Chl <-        normalize(M_heatmap[2])
M_heatmap$daylength <-  normalize(M_heatmap[3])
M_heatmap$dp10m <-      normalize(M_heatmap[4])
head(M_heatmap)
library(RColorBrewer)
library(viridis)
Colors = grey.colors(10)
Colors=rev(viridis_pal( begin=.15, end=0.8, option ="inferno")(12))

#x11()
png(file="Figure 6b.png",width=15, height=10, units="cm", res=300)

heatmap.2(as.matrix(M_heatmap[1:4]), Rowv=Rowv, Colv=NA, 
          dendrogram = c("row"), col = Colors, density.info = "none", 
          srtCol=45,
          key.ylab =1.5, keysize = 2, key.title="scaled", 
          #lmat = lmat2, lwid = c(1.5,4),lhei = c(1.5,4,1),
          cexRow=.85, cexCol=.85, xaxt="n", tracecol = "black",
          RowSideColors = c(    # grouping row-variables into different
            rep("grey35", 6),   # categories, Measurement 1-3: green
            rep("magenta", 8),    # Measurement 4-8: blue
            rep("#1B9E77", 7)), margin=c(4,8))
par(lend = 1, xpd=T)           # square line ends for the color legend
legend(0.75,1.25,legend = c( "May-Jul", "Aug-Oct", "Oct-Mar"), 
       col = c("#1B9E77","grey25","magenta"),
       bty ="n", lty= 1, lwd = 10 )           # line width)
par(xpd=F)
dev.off()

######