## Morphological analysis of linear traits

# Packages ----
library(geomorph); library(dplyr); library(ggplot2);library(car);library(dunn.test);library(Morpho);library(dabestr)
library(ggsignif)
# Functions ----
# Consistent plot aesthetics for PCA
theme.pca <- function() {
  theme_bw() +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          panel.background = element_rect(colour="black",fill="white", size=1),
          axis.text = element_text(size=16, color = "black"),
          axis.ticks = element_line(size = 0.5, colour = "black"),
          axis.ticks.length = unit(3, "mm"),
          axis.title.y = element_text(size = 30),
          axis.title.x = element_text(size = 30),
          axis.text.x = element_text(size=20),
          axis.text.y = element_text(size=20),
          legend.title = element_text(size = 20),
          legend.text = element_text(size = 20))
}

theme.boxplot <- function() {
  theme_bw() +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          panel.background = element_rect(colour="black",fill="white", size=1),
          axis.text = element_text(size=16, color = "black"),
          axis.ticks = element_line(size = 0.5, colour = "black"),
          axis.ticks.length = unit(3, "mm"),
          axis.title.y = element_text(size = 15),
          axis.title.x = element_text(size = 15),
          axis.text.x = element_text(size=10),
          axis.text.y = element_text(size=10),
          legend.position = "none")
}

cols <- c("Allt na lairige"="#a4dede","Carron"="#008396","Eck"="gray60","Glashan"="#9b2915","Lomond"="gray60","Shira"="#287c71",
          "Sloy"="#6c966f","Tarsan"="#fe7c73")

# Set working directory ----
setwd("~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Morphometrics_ecology")

# read data ----
act <- read.csv("allometry_corrected_traits_20200421.csv",header=TRUE)
lake_names <- dplyr::recode(act$Lake, Lairige = "Allt na lairige")
act$Lake <- lake_names
act$Lake <- factor(act$Lake, levels = c("Eck","Glashan", "Tarsan","Allt na lairige", "Shira","Carron","Sloy","Lomond"))


# split data into system ----
lom <- filter(act, System == "Lomond")
lom$Lake <- factor(lom$Lake, levels = c("Lomond","Allt na lairige","Shira","Carron","Sloy"))


eck <- subset(act, System == "Eck")
eck$Lake <- factor(eck$Lake, levels = c("Eck","Glashan","Tarsan"))


# Dunn test analysis ----

# Lom
model_ml_lom <- dunn.test(lom$ML, lom$Lake, method = "bh")
model_ed_lom <- dunn.test(lom$ED, lom$Lake, method = "bh")
model_sl_lom <- dunn.test(lom$SL, lom$Lake, method = "bh")
model_hl_lom <- dunn.test(lom$HL, lom$Lake, method = "bh")
model_hd_lom <- dunn.test(lom$HD, lom$Lake, method = "bh")
model_bda_lom <- dunn.test(lom$BDA, lom$Lake, method = "bh")
model_bdp_lom <- dunn.test(lom$BDP, lom$Lake, method = "bh")
model_cpd_lom <- dunn.test(lom$CPD, lom$Lake, method = "bh")
model_fl_lom <- dunn.test(lom$FL, lom$Lake, method = "bh")

# Eck
model_ml_eck <- dunn.test(eck$ML, eck$Lake, method = "bh")
model_ed_eck <- dunn.test(eck$ED, eck$Lake, method = "bh")
model_sl_eck <- dunn.test(eck$SL, eck$Lake, method = "bh")
model_hl_eck <- dunn.test(eck$HL, eck$Lake, method = "bh")
model_hd_eck <- dunn.test(eck$HD, eck$Lake, method = "bh")
model_bda_eck <- dunn.test(eck$BDA, eck$Lake, method = "bh")
model_bdp_eck <- dunn.test(eck$BDP, eck$Lake, method = "bh")
model_cpd_eck <- dunn.test(eck$CPD, eck$Lake, method = "bh")
model_fl_eck <- dunn.test(eck$FL, eck$Lake, method = "bh")

# between systems
(ml_lom_eck <- kruskal.test(ML ~ System, data = act))
(ed_lom_eck <- kruskal.test(ED ~ System, data = act))
(sl_lom_eck <- kruskal.test(SL ~ System, data = act))
(hl_lom_eck <- kruskal.test(HL ~ System, data = act))
(hd_lom_eck <- kruskal.test(HD ~ System, data = act))
(bda_lom_eck <- kruskal.test(BDA ~ System, data = act))
(bdp_lom_eck <- kruskal.test(BDP ~ System, data = act))
(cpd_lom_eck <- kruskal.test(CPD ~ System, data = act))
(fl_lom_eck <- kruskal.test(FL ~ System, data = act))



# plot 

ml_plot <- ggplot(act, aes(x = Lake, y = ML, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Maxilla length") +
  geom_signif(comparisons = list(c("Eck","Tarsan")),map_signif_level = TRUE, colour = "black", y_position = -1.35, textsize = 5) +
  geom_signif(comparisons = list(c("Eck","Glashan")),map_signif_level = TRUE, colour = "black", y_position = -1.3, textsize = 5) 


ed_plot <- ggplot(act, aes(x = Lake, y = ED, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Eye diameter")

sl_plot <- ggplot(act, aes(x = Lake, y = SL, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Snout length")

hl_plot <- ggplot(act, aes(x = Lake, y = HL, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Head length")

hd_plot <- ggplot(act, aes(x = Lake, y = HD, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Head depth")

bda_plot <- ggplot(act, aes(x = Lake, y = BDA, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Body depth anterior")

bdp_plot <- ggplot(act, aes(x = Lake, y = BDP, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Body depth posterior")

cpd_plot <- ggplot(act, aes(x = Lake, y = CPD, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Caudal peduncule depth")

fl_plot <- ggplot(act, aes(x = Lake, y = FL, colour = Lake, shape = System)) + geom_boxplot(outlier.shape = NA) + 
  geom_jitter(size = 3, alpha = 0.85, width = 0.1) + theme.boxplot() + 
  scale_color_manual(values = cols) + labs(y = "Fin length")

panel <- gridExtra::grid.arrange(ml_plot, ed_plot, sl_plot, hl_plot, hd_plot, bda_plot, bdp_plot, cpd_plot, fl_plot, ncol = 3, nrow = 3)

ggsave("./figures/linear_traits.pdf", panel, device = "pdf", width = 18, height = 12, units = "in", limitsize = FALSE)

ggplot(act, aes(x = BL, y = ED, colour = System, shape = System)) + theme_bw() + geom_point(size = 5) + geom_smooth(method = "lm")