#set the woking directory to the location of your files
setwd("D:/Users/elias/Documents/R/Internship_2")

#install and open all necessary packages
install.packages("remotes")
install.packages("ggfortify")
remotes::install_github("gavinsimpson/ggvegan", force=TRUE)
library(vegan)
library(readxl)
library(ggplot2)
library(ggfortify)




#Part 1a: PCA on coral dataset based on all timepoints

#read the data and check if the data is correctly imported
data_coral <- read_xlsx("Edited_Data/datacoral.xlsx")
summary(data_coral)

#perform the PCA and plot the graph
PCA <- prcomp(data_coral[, c(5:31)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.1 | a[,1] > 0.1) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.1 | a[,2] > 0.1) 
a4 <- rownames(a3)
plot1 <- autoplot(PCA, data= data_coral, colour= 'Zone', shape='Time', size= 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a1[,1]/5), y = (a1[,2]/5),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a3[,1]/5), y = (a3[,2]/5),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery", "Peak-bleaching","Pre-bleaching", "Short-term recovery")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot1

rm(a, a1, a2, a3, a4, PCA, data_coral)


#Part 1b: PCA on coral dataset based on January 2016 and October 2019 timepoints

#read the data and check if the data is correctly imported
data_coral_jan2019 <- read_xlsx("Edited_Data/datacoral_JAN2019.xlsx")
summary(data_coral_jan2019)

#perform the PCA and plot the graph
PCA <- prcomp(data_coral_jan2019[, c(5:31)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.2 | a[,1] > 0.2) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.2 | a[,2] > 0.2) 
a4 <- rownames(a3)
plot2 <- autoplot(PCA, data= data_coral_jan2019, colour= 'Zone', shape='Time', size= 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a1[,1]/5), y = (a1[,2]/5),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a3[,1]/5), y = (a3[,2]/5),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery","Pre-bleaching")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot2

rm(a, a1, a2, a3, a4, PCA, data_coral_jan2019)

#Part 2a: PCA on benthic dataset based on all timepoints

#read the data and check if the data is correctly imported
data_benthos <- read_xlsx("Edited_Data/databenthos.xlsx")
summary(data_benthos)

#perform the PCA and plot the graph
PCA <- prcomp(data_benthos[, c(5:10)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.1 | a[,1] > 0.1) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.1 | a[,2] > 0.1) 
a4 <- rownames(a3)
plot3 <- autoplot(PCA, data= data_benthos, colour= 'Zone', shape='Time', size= 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a1[,1]/3), y = (a1[,2]/3),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a3[,1]/3), y = (a3[,2]/3),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery", "Peak-bleaching","Pre-bleaching", "Short-term recovery")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot3

rm(a, a1, a2, a3, a4, PCA, data_benthos)

#Part 2b: PCA on benthic dataset based on January 2016 and October 2019 timepoints

#read the data and check if the data is correctly imported
data_benthos_jan2019 <- read_xlsx("Edited_Data/databenthos_JAN2019.xlsx")
summary(data_benthos_jan2019)

#perform the PCA and plot the graph
PCA <- prcomp(data_benthos_jan2019[, c(5:10)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.1 | a[,1] > 0.1) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.1 | a[,2] > 0.1) 
a4 <- rownames(a3)
plot4 <- autoplot(PCA, data= data_benthos_jan2019, colour= 'Zone', shape='Time', size= 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a1[,1]/3), y = (a1[,2]/3),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a3[,1]/3), y = (a3[,2]/3),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery","Pre-bleaching")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot4

rm(a, a1, a2, a3, a4, PCA, data_benthos_jan2019)


#Part 3a: PCA on life history stratagies based all timepoints

#read the data and check if the data is correctly imported
data_lhs <- read_xlsx("Edited_Data/datalhs_updated.xlsx")
summary(data_lhs)

#perform the PCA and plot the graph
PCA <- prcomp(data_lhs[, c(5:7)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.2 | a[,1] > 0.2) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.2 | a[,2] > 0.2) 
a4 <- rownames(a3)
plot5 <- autoplot(PCA, data= data_lhs, colour= 'Zone', shape='Time', size = 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a1[,1]/3), y = (a1[,2]/3),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a3[,1]/3), y = (a3[,2]/3),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery", "Peak-bleaching","Pre-bleaching", "Short-term recovery")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot5

rm(a, a1, a2, a3, a4, PCA, data_lhs)

#Part 3b: PCA on life history strategies based on January 2016 and October 2019 timepoints

#read the data and check if the data is correctly imported
data_lhs_jan2019 <- read_xlsx("Edited_Data/datalhs_JAN2019_updated.xlsx")
summary(data_lhs_jan2019)

#perform the PCA and plot the graph
PCA <- prcomp(data_lhs_jan2019[, c(5:7)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.2 | a[,1] > 0.2) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.2 | a[,2] > 0.2) 
a4 <- rownames(a3)
plot6 <- autoplot(PCA, data= data_lhs_jan2019, colour= 'Zone', shape='Time', size = 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a1[,1]/3), y = (a1[,2]/3),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/3),yend = (PC2/3)),
               color = "red") +
  annotate("text", x = (a3[,1]/3), y = (a3[,2]/3),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery", "Pre-bleaching")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot6

rm(a, a1, a2, a3, a4, PCA, data_lhs_jan2019)

#Part 4a: PCA on growth forms based on all timepoints

#read the data and check if the data is correctly imported
data_growthform <- read_xlsx("Edited_Data/datagrowthform.xlsx")
summary(data_growthform)

#perform the PCA and plot the graph
PCA <- prcomp(data_growthform[, c(4:8)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.2 | a[,1] > 0.2) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.2 | a[,2] > 0.2) 
a4 <- rownames(a3)
plot7 <- autoplot(PCA, data= data_growthform, colour= 'Zone', shape='Time', size = 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a1[,1]/5), y = (a1[,2]/5),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a3[,1]/5), y = (a3[,2]/5),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery", "Peak-bleaching","Pre-bleaching", "Short-term recovery")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot7

rm(a, a1, a2, a3, a4, PCA, data_growthform)

#Part 4b: PCA on growth forms based on January 2016 and October 2019 timepoints

#read the data and check if the data is correctly imported
data_growthform_jan2019 <- read_xlsx("Edited_Data/datagrowthform_JAN2019.xlsx")
summary(data_growthform_jan2019)

#perform the PCA and plot the graph
PCA <- prcomp(data_growthform_jan2019[, c(4:8)], scale=FALSE)
a <- PCA$rotation
a1 <- subset(a, a[,1] < -0.2 | a[,1] > 0.2) 
a2 <- rownames(a1)
a3 <- subset(a, a[,2] < -0.2 | a[,2] > 0.2) 
a4 <- rownames(a3)
plot8 <- autoplot(PCA, data= data_growthform_jan2019, colour= 'Zone', shape='Time', size = 5) +
  geom_point(size=0.75, colour='black') +
  geom_segment(data = a1, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a1[,1]/5), y = (a1[,2]/5),
           label = a2) +
  geom_segment(data = a3, aes(x = 0, y = 0, xend = (PC1/5),yend = (PC2/5)),
               color = "red") +
  annotate("text", x = (a3[,1]/5), y = (a3[,2]/5),
           label = a4) +
  scale_shape(name="Time", labels= c("Long-term recovery", "Pre-bleaching")) +
  theme(legend.direction = "horizontal", legend.box = "horizontal", legend.position="bottom")+
  theme(legend.background = element_rect(fill="lightblue", size=0.5, linetype="solid", colour ="darkblue"))
plot8

rm(a, a1, a2, a3, a4, PCA, data_growthform_jan2019)


#all plots individually
plot1
plot2
plot3
plot4
plot5
plot6
plot7
plot8