## Geometric morphometric analysis

# 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))
}

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

# library ----
library(geomorph);library(tidyverse)

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

# read data ----
fish <- readland.tps("unbent2.TPS",specID="imageID") #read species name from photo's name
str(fish)
classifier <- read.csv("classifier2.csv",header = T) # read classifier table

# data manupulation ----
#Remove landmarks used for linear measurements and other unnecessary ones
omit3 <- c(3,4,5,6,7,10,25,23,24,19,15)
fish2 <- fish[-omit3,,]

### Generalised Procrustes Analysis
fish.gpa <- gpagen(fish2)
plotAllSpecimens(fish.gpa$coords)

# split into separate systems
eck.coords <- fish.gpa$coords[,,c(1:118,142:213,239:271,397:456)]
lom.coords <- fish.gpa$coords[,,c(119:141,214:238,272:396,457:508)]

###SUBSET classifier table
ECK.table <- subset(classifier, System == "Eck")
LOM.table <- subset(classifier, System == "Lomond")

# create goemorph data frames
gdf.eck <- geomorph.data.frame(Shape = eck.coords, Lake = ECK.table$Lake, Size = fish.gpa$Csize[c(1:118,142:213,239:271,397:456)], Population = ECK.table$Population)
gdf.lom <- geomorph.data.frame(Shape = lom.coords, Lake = LOM.table$Lake, Size = fish.gpa$Csize[c(119:141,214:238,272:396,457:508)], Population = LOM.table$Population)
gdf.lakes <- geomorph.data.frame(Shape = fish.gpa$coords, Lake = classifier$Lake, Size = fish.gpa$Csize, 
                                 Habitat = classifier$Habitat, System = classifier$System, Population = classifier$Population)

# check allometry ----
eck.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.eck,iter=999,alpha = 0.01)
eck.allo; plot(eck.allo, method="PredLine", warpgrids=FALSE)

lom.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.lom,iter=999,alpha = 0.05)
lom.allo; plot(lom.allo, method="PredLine", warpgrids=FALSE)

# check allometry within eck system separately ----
# eck and tar
eck.tar.coords <- fish.gpa$coords[,,c(1:118,172:213,397:456)]
ECK.TAR.table <- subset(classifier, Lake == "Eck" | Lake == "Tarsan")
gdf.eck.tar <- geomorph.data.frame(Shape = eck.tar.coords, Lake = ECK.TAR.table$Lake, Size = fish.gpa$Csize[c(1:118,172:213,397:456)])

eck.tar.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.eck.tar,iter=999,alpha = 0.01)
eck.tar.allo; plot(eck.tar.allo, method="PredLine", warpgrids=FALSE)

# check allometry within lomond system separately ----
### test lomond lakes separately
# lom  and car
lom.car.coords <- fish.gpa$coords[,,c(214:238,276:364,477:508)]
LOM.CAR.table <- subset(classifier, Lake == "Lomond" | Lake == "Carron")
gdf.lom.car <- geomorph.data.frame(Shape = lom.car.coords, Lake = LOM.CAR.table$Lake, Size = fish.gpa$Csize[c(214:238,276:364,477:508)])

lom.car.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.lom.car,iter=999,alpha = 0.05)
lom.car.allo; plot(lom.car.allo, method="PredLine", warpgrids=FALSE)

# lom and sloy
lom.sloy.coords <- fish.gpa$coords[,,c(214:221,276:381,457:497)]
LOM.SLOY.table <- subset(classifier, Lake == "Lomond" | Lake == "Sloy")
gdf.lom.sloy <- geomorph.data.frame(Shape = lom.sloy.coords, Lake = LOM.SLOY.table$Lake, Size = fish.gpa$Csize[c(214:221,276:381,457:497)])

lom.sloy.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.lom.sloy,iter=999,alpha = 0.05)
lom.sloy.allo; plot(lom.sloy.allo, method="PredLine", warpgrids=FALSE)

# lom and shira
lom.shi.coords <- fish.gpa$coords[,,c(128:141,214:221,276:364,382:396,477:497)]
LOM.SHI.table <- subset(classifier, Lake == "Lomond" | Lake == "Shira")
gdf.lom.shi <- geomorph.data.frame(Shape = lom.shi.coords, Lake = LOM.SHI.table$Lake, Size = fish.gpa$Csize[c(128:141,214:221,276:364,382:396,477:497)])

lom.shi.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.lom.shi,iter=999,alpha = 0.05)
lom.shi.allo; plot(lom.shi.allo, method="PredLine", warpgrids=FALSE)

# lom and allt na lairige
lom.lai.coords <- fish.gpa$coords[,,c(119:127,214:221,272:364,477:497)]
LOM.LAI.table <- subset(classifier, Lake == "Lomond" | Lake == "Lairige")
gdf.lom.lai <- geomorph.data.frame(Shape = lom.lai.coords, Lake = LOM.LAI.table$Lake, Size = fish.gpa$Csize[c(119:127,214:221,272:364,477:497)])

lom.lai.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.lom.lai,iter=999,alpha = 0.05)
lom.lai.allo; plot(lom.lai.allo, method="PredLine", warpgrids=FALSE)

# lom car sloy

lom.car.sloy.coords <- fish.gpa$coords[,,c(214:238,276:381,457:508)] 
LOM.CAR.SLOY.table <- subset(classifier, Lake == "Lomond" | Lake == "Carron" | Lake == "Sloy")
gdf.lom.car.sloy <- geomorph.data.frame(Shape = lom.car.sloy.coords, Lake = LOM.CAR.SLOY.table$Lake, Size = fish.gpa$Csize[c(214:238,276:381,457:508)])

lom.car.sloy.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.lom.car.sloy,iter=999,alpha = 0.05)
lom.car.sloy.allo; plot(lom.car.sloy.allo, method="PredLine", warpgrids=FALSE)

# lom car sloy lai

lom.car.sloy.lai.coords <- fish.gpa$coords[,,c(119:127,214:238,272:381,457:508)] 
LOM.CAR.SLOY.LAI.table <- subset(classifier, Lake == "Lomond" | Lake == "Carron" | Lake == "Sloy" | Lake == "Lairige")
gdf.lom.car.sloy.lai <- geomorph.data.frame(Shape = lom.car.sloy.lai.coords, Lake = LOM.CAR.SLOY.LAI.table$Lake, Size = fish.gpa$Csize[c(119:127,214:238,272:381,457:508)])

lom.car.sloy.lai.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.lom.car.sloy.lai,iter=999,alpha = 0.05)
lom.car.sloy.lai.allo; plot(lom.car.sloy.lai.allo, method="PredLine", warpgrids=FALSE)

# car shi sloy

car.shi.coords <- fish.gpa$coords[,,c(128:141,222:238,365:396,457:476,498:508)] 
CAR.SHI.table <- subset(classifier, Lake == "Carron" | Lake == "Shira" | Lake == "Sloy")
gdf.car.shi <- geomorph.data.frame(Shape = car.shi.coords, Lake = CAR.SHI.table$Lake, Size = fish.gpa$Csize[c(128:141,222:238,365:396,457:476,498:508)])

car.shi.allo <- procD.allometry(Shape ~ log(Size), ~Lake, data = gdf.car.shi,iter=999,alpha = 0.05)
car.shi.allo; plot(car.shi.allo, method="PredLine", warpgrids=FALSE)


# pairwise comparison ----
###Pairwise comparison
syst.pairwise <- advanced.procD.lm(Shape~log(Size)+System, ~log(Size)*System,iter=999,groups= ~System,slope= ~log(Size),data=gdf.lakes)
syst.pairwise

eck.pairwise <- advanced.procD.lm(Shape~log(Size)+Lake,~log(Size)*Lake,groups=~Lake,slope = ~log(Size),iter=999,data=gdf.eck)
eck.pairwise

# comparison uncluding Shira and Lake*Size interaction term
lom1.pairwise <- advanced.procD.lm(Shape~log(Size)+Lake,~log(Size)*Lake,groups=~Lake,slope = ~log(Size),data=gdf.lom,iter=999)
lom1.pairwise

# comparison excluding shira, comparing mean shape
lom2.pairwise <- advanced.procD.lm(Shape~log(Size),~log(Size)+Lake,groups=~Lake,data=gdf.lom.car.sloy.lai,iter=999)
lom2.pairwise

# PCA ----
syst.pca <- plotTangentSpace(fish.gpa$coords, groups = classifier$System, legend = TRUE, warpgrids = TRUE, axis1 = 1, axis2 = 2)
eck.pca <- plotTangentSpace(eck.coords, groups = ECK.table$Lake, legend = TRUE,warpgrids = TRUE, axis1 = 1, axis2 = 2)
lom.pca <- plotTangentSpace(lom.coords, groups = LOM.table$Lake, legend = TRUE,warpgrids = TRUE, axis1 = 1, axis2 = 2)
lakes.pca <- plotTangentSpace(fish.gpa$coords, groups = classifier$Lake, legend = TRUE, axis1 = 1, axis2 = 2)

# ggplot PCA plotting ----
lakes_new_df <- data.frame(classifier$Lake,lakes.pca$pc.scores)


df_pca_lakes <- as.data.frame(lakes.pca$pc.scores[,c(1:4)]);df_pca_lakes[,5] <- classifier$Lake; df_pca_lakes[,6] <- classifier$System
colnames(df_pca_lakes)[5] <- "Lake";colnames(df_pca_lakes)[6] <- "System"
lake_names <- dplyr::recode(df_pca_lakes$Lake, Lairige = "Allt na lairige")
df_pca_lakes$Lake <- lake_names
df_pca_lakes2 <- df_pca_lakes %>% rownames_to_column("ID")
write.csv(df_pca_lakes2, "morphology_combined_pca.csv", row.names = FALSE)
df_pca_lakes$Lake <- factor(df_pca_lakes$Lake, levels = c("Glashan", "Tarsan","Eck","Allt na lairige", "Shira","Carron","Sloy","Lomond"))
df_pca_lakes <- df_pca_lakes %>% 
  mutate(Population = case_when(
    grepl("Eck", Lake) ~ "source",
    grepl("Glashan", Lake) ~ "refuge",
    grepl("Tarsan", Lake) ~ "refuge",
    grepl("Allt na lairige", Lake) ~ "refuge",
    grepl("Shira", Lake) ~ "refuge",
    grepl("Carron", Lake) ~ "refuge",
    grepl("Sloy", Lake) ~ "refuge",
    grepl("Lomond", Lake) ~ "source"))


# create a datatframe of mean pc scores
lakes_pc_mean <- plyr::ddply(df_pca_lakes,c("Lake","System","Population"), summarise,
                       PC1.mean = mean(PC1),
                       PC1.se = sd(PC1)/sqrt(length(PC1)),
                       PC2.mean = mean(PC2),
                       PC2.se = sd(PC2)/sqrt(length(PC2)),
                       PC3.mean = mean(PC3),
                       PC3.se = sd(PC3)/sqrt(length(PC3)),
                       PC4.mean = mean(PC4),
                       PC4.se = sd(PC4)/sqrt(length(PC4)))

# Plot
pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Morphometrics_ecology/figures/pca_geomorph_12_lake.pdf", width = 11.69, height = 8.27)
ggplot() + geom_point(data = lakes_pc_mean, aes(x = PC1.mean, y = PC2.mean,color = Lake, shape = System),size = 8, alpha = 0.85) + 
  geom_errorbarh(data = lakes_pc_mean, aes(xmax = PC1.mean + PC1.se, xmin = PC1.mean - PC1.se, height = 0, y = PC2.mean)) +
  geom_errorbar(data = lakes_pc_mean, aes(ymax = PC2.mean + PC2.se, ymin = PC2.mean - PC2.se, width=0, x = PC1.mean)) + theme.pca() +
  labs(x = "PC1 (22%)",y = "PC2 (20%)") + scale_color_manual(values = cols)
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Morphometrics_ecology/figures/pca_geomorph_13_lake.pdf", width = 11.69, height = 8.27)
ggplot() + geom_point(data = lakes_pc_mean, aes(x = PC1.mean, y = PC3.mean,color = Lake, shape = System),size = 8, alpha = 0.85) + 
  geom_errorbarh(data = lakes_pc_mean, aes(xmax = PC1.mean + PC1.se, xmin = PC1.mean - PC1.se, height = 0, y = PC3.mean)) +
  geom_errorbar(data = lakes_pc_mean, aes(ymax = PC3.mean + PC3.se, ymin = PC3.mean - PC3.se, width=0, x = PC1.mean)) + theme.pca() +
  labs(x = "PC1 (22%)",y = "PC3 (12%)") + scale_color_manual(values = cols)
dev.off()

ggplot() + geom_point(data = lakes_pc_mean, aes(x = PC3.mean, y = PC4.mean,color = Lake, shape = System),size = 8, alpha = 0.85) + 
  geom_errorbarh(data = lakes_pc_mean, aes(xmax = PC3.mean + PC3.se, xmin = PC3.mean - PC3.se, height = 0, y = PC4.mean)) +
  geom_errorbar(data = lakes_pc_mean, aes(ymax = PC4.mean + PC4.se, ymin = PC4.mean - PC4.se, width=0, x = PC3.mean)) + theme.pca() +
  labs(x = "PC3 (22%)",y = "PC4 (8%)") + scale_color_manual(values = cols)

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Morphometrics_ecology/figures/pca_geomorph_12_laketype.pdf", width = 11.69, height = 8.27)
ggplot() + geom_point(data = lakes_pc_mean, aes(x = PC1.mean, y = PC2.mean,color = Population, shape = System),size = 8, alpha = 0.85) + 
  geom_errorbarh(data = lakes_pc_mean, aes(xmax = PC1.mean + PC1.se, xmin = PC1.mean - PC1.se, height = 0, y = PC2.mean)) +
  geom_errorbar(data = lakes_pc_mean, aes(ymax = PC2.mean + PC2.se, ymin = PC2.mean - PC2.se, width=0, x = PC1.mean)) + theme.pca() +
  labs(x = "PC1 (22%)",y = "PC2 (20%)") + scale_colour_brewer(palette = "Dark2")
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Morphometrics_ecology/figures/pca_geomorph_13_laketype.pdf", width = 11.69, height = 8.27)
ggplot() + geom_point(data = lakes_pc_mean, aes(x = PC1.mean, y = PC3.mean,color = Population, shape = System),size = 8, alpha = 0.85) + 
  geom_errorbarh(data = lakes_pc_mean, aes(xmax = PC1.mean + PC1.se, xmin = PC1.mean - PC1.se, height = 0, y = PC3.mean)) +
  geom_errorbar(data = lakes_pc_mean, aes(ymax = PC3.mean + PC3.se, ymin = PC3.mean - PC3.se, width=0, x = PC1.mean)) + theme.pca() +
  labs(x = "PC1 (22%)",y = "PC3 (12%)") + scale_colour_brewer(palette = "Dark2")
dev.off()


# trajectory analysis ----

TA <- trajectory.analysis(Shape ~ System * Population, data = gdf.lakes)
summary(TA, angle.type = "deg")
plot(TA,pt.scale=1.5)
TA$means


# reference shape graphs ----
## Make nice graphs
vector <- c(1,2,2,5,5,9,9,10,10,11,11,12,12,13,13,14,14,7,7,4,4,1,5,6,6,7)
matrix <- matrix(vector,13,byrow = TRUE)
matrix

ref <- mshape(fish.gpa$coords)
define.links(ref)

plotRefToTarget(M1=ref,M2=lakes.pca$pc.shapes$PC1min,links=matrix,method = "points",mag = 2)
plotRefToTarget(M1=ref,M2=lakes.pca$pc.shapes$PC1max,links=matrix,method = "points",mag = 2)
plotRefToTarget(M1=ref,M2=lakes.pca$pc.shapes$PC4max,links=matrix,method = "points",mag = 2)
plotRefToTarget(M1=ref,M2=lakes.pca$pc.shapes$PC4min,links=matrix,method = "points",mag = 2)



# Extracting linear measurements from file ----
# Needed for linear measurements analysis
fish2.gpa <- gpagen(fish)
plot(fish2.gpa)

# Make a dataframe defining three interlandmark distances, calculate distances and write as csv

lmks <- data.frame(ML = c(1,2), HL = c(1,12), ED = c(4,5), SL = c(1,7), HD = c(9,10), 
                   BDA = c(22,16), BDP = c(21,17), CPD = c(20,18), FL = c(14,15), OL = c(25,11), 
                   BL = c(1,19), row.names = c("start", "end")) 
A <- fish2.gpa$coords
lineardists.WLC <- interlmkdist(A, lmks)

write.csv(lineardists.WLC, "LinearMeasurements_Whitefish_raw.csv")