library(tidyverse)

# 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 = 15),
          axis.title.x = element_text(size = 15),
          axis.text.x = element_text(size=10),
          axis.text.y = element_text(size=10),
          legend.title = element_text(size = 10),
          legend.text = element_text(size = 10))
}

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

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

source("~/Desktop/1-s2.0-S1470160X15001569-mmc1/si_div.R")
# importing stable isotope values (d13C and d15N) for 97 individuals fish (different sampling size between lake)

# muscle data
individuals_si <- read.csv("muscle_isotope_2020.csv")
individuals_si$Lake <- recode(individuals_si$Lake, Allt_na_lairige = "Allt na lairige")

summary(individuals_si[,"Lake"]) # number of individuals per lake

# stomach content data
gut_si <- read.csv("gut_isotope_2020.csv")
gut_si$Lake <- recode(gut_si$Lake, Allt_na_lairige = "Allt na lairige")
summary(gut_si[,"Lake"]) # number of individuals per lake


# combined muscle and stomach content for overlap analysis

individuals_si$tissue <- "muscle"
gut_si$tissue <- "stomach_content"

combined_si <- rbind(individuals_si, gut_si)

# computing mean Stable Isotope values for each lake
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per lake did not mirror actual lake biomass
individuals_si <- data.frame(group = individuals_si[,"Lake"], individuals_si)
mean_si_lake <- meanSI_group(individuals_si)

gut_si <- data.frame(group = gut_si[,"Lake"], gut_si)
mean_gut_si_lake <- meanSI_group(gut_si)

combined_si <- data.frame(group = combined_si[, "Lake"], combined_si)
mean_combined_si_lake <- meanSI_group(combined_si)

# computing coefficent of variation within each lake to assess intraspecific variability
cbind( CV_d13C=mean_si_lake[,"sd_d13C"]/mean_si_lake[,"d13C"], CV_d15N=mean_si_lake[,"sd_d15N"]/mean_si_lake[,"d15N"] )

cbind( CV_d13C=mean_gut_si_lake[,"sd_d13C"]/mean_gut_si_lake[,"d13C"], CV_d15N=mean_gut_si_lake[,"sd_d15N"]/mean_gut_si_lake[,"d15N"] )


# Example of how computing indices at the lake level based on individuals isotopic values
data_ind_scl <- scaleSI_range01(individuals_si)
summary(data_ind_scl)

data_gut_scl <- scaleSI_range01(gut_si)

data_combined_scl <- scaleSI_range01(combined_si)

find_hull <- function(df) df[chull(df$d13C, df$d15N), ]

# plot together muscle and gut isotopic niches
df_muscle <- data_ind_scl
df_muscle$Lake <- factor(df_muscle$Lake, levels = c("Glashan", "Tarsan", "Allt na lairige", "Shira", "Carron", "Sloy", "Lomond"))

hulls_muscle <- plyr::ddply(df_muscle, "Lake", find_hull)
muscle_plot <- ggplot(df_muscle, aes(x = d13C, y = d15N, fill = Lake, colour = Lake)) + theme.pca() +
  geom_point(size = 3) + geom_polygon(data = hulls_muscle, alpha = 0.3) + scale_fill_manual(values = cols) + 
  scale_colour_manual(values = cols) + labs(x =  expression(delta^{13}*C), y = expression(delta^{15}*N)) + xlim(0,1) + ylim(0,1)

df_gut <- data_gut_scl
hulls_gut <-  plyr::ddply(df_gut, "Lake", find_hull)
gut_plot <- ggplot(df_gut, aes(x = d13C, y = d15N, fill = Lake, colour = Lake)) + theme.pca() +
  geom_point(size = 3) + geom_polygon(data = hulls_gut, alpha = 0.3) + scale_fill_manual(values = cols) + 
  scale_colour_manual(values = cols) + labs(x =  expression(delta^{13}*C), y = expression(delta^{15}*N)) + xlim(0,1) + ylim(0,1)


df_combined <- data_combined_scl
df_combined$Lake <- factor(df_combined$Lake, levels = c("Glashan", "Tarsan", "Allt na lairige", "Shira", "Carron", "Sloy", "Lomond"))
hulls_combined <- plyr::ddply(df_combined, c("Lake","tissue"), find_hull)



# computing isotopic diversity for each lake based on individuals values (no weighting by fish mass) ----
IDLAI_scl_ind<-IDiversity(cons=data_ind_scl[which( data_ind_scl[,"group"]=="Allt na lairige"),c("d13C","d15N")], col="#a4dede", nm_plot="Figure ID_LAI", scaled=T)
IDSHI_scl_ind<-IDiversity(cons=data_ind_scl[which( data_ind_scl[,"group"]=="Shira"),c("d13C","d15N")], col="#287c71", nm_plot="Figure ID_SHI", scaled=T)
IDCAR_scl_ind<-IDiversity(cons=data_ind_scl[which( data_ind_scl[,"group"]=="Carron"),c("d13C","d15N")], col="#008396", nm_plot="Figure ID_CAR", scaled=T)
IDSLO_scl_ind<-IDiversity(cons=data_ind_scl[which( data_ind_scl[,"group"]=="Sloy"),c("d13C","d15N")], col="#6c966f", nm_plot="Figure ID_SLO", scaled=T)
IDLOM_scl_ind<-IDiversity(cons=data_ind_scl[which( data_ind_scl[,"group"]=="Lomond"),c("d13C","d15N")], col="#999999", nm_plot="Figure ID_LOM", scaled=T)
IDGLA_scl_ind<-IDiversity(cons=data_ind_scl[which( data_ind_scl[,"group"]=="Glashan"),c("d13C","d15N")], col="#9b2915", nm_plot="Figure ID_GLA", scaled=T)
IDTAR_scl_ind<-IDiversity(cons=data_ind_scl[which( data_ind_scl[,"group"]=="Tarsan"),c("d13C","d15N")], col="#fe7c73", nm_plot="Figure ID_TAR", scaled=T)

IDLAI_scl_gut<-IDiversity(cons=data_gut_scl[which( data_gut_scl[,"group"]=="Allt na lairige"),c("d13C","d15N")], col="#a4dede", nm_plot="Figure ID_LAI", scaled=T)
IDSHI_scl_gut<-IDiversity(cons=data_gut_scl[which( data_gut_scl[,"group"]=="Shira"),c("d13C","d15N")], col="#287c71", nm_plot="Figure ID_SHI", scaled=T)
IDCAR_scl_gut<-IDiversity(cons=data_gut_scl[which( data_gut_scl[,"group"]=="Carron"),c("d13C","d15N")], col="#008396", nm_plot="Figure ID_CAR", scaled=T)
IDSLO_scl_gut<-IDiversity(cons=data_gut_scl[which( data_gut_scl[,"group"]=="Sloy"),c("d13C","d15N")], col="#6c966f", nm_plot="Figure ID_SLO", scaled=T)
IDLOM_scl_gut<-IDiversity(cons=data_gut_scl[which( data_gut_scl[,"group"]=="Lomond"),c("d13C","d15N")], col="#999999", nm_plot="Figure ID_LOM", scaled=T)
IDGLA_scl_gut<-IDiversity(cons=data_gut_scl[which( data_gut_scl[,"group"]=="Glashan"),c("d13C","d15N")], col="#9b2915", nm_plot="Figure ID_GLA", scaled=T)
IDTAR_scl_gut<-IDiversity(cons=data_gut_scl[which( data_gut_scl[,"group"]=="Tarsan"),c("d13C","d15N")], col="#fe7c73", nm_plot="Figure ID_TAR", scaled=T)


iso_diversity <- cbind(LAI = IDLAI_scl_ind , SHI = IDSHI_scl_ind, CAR = IDCAR_scl_ind, SLO = IDSLO_scl_ind, LOM = IDLOM_scl_ind, GLA = IDGLA_scl_ind, TAR =IDTAR_scl_ind)
write.table(iso_diversity, "~/Desktop/1-s2.0-S1470160X15001569-mmc1/isotopic_diversity.txt", sep = "\t", quote = FALSE)

iso_gut_diversity <- cbind(LAI = IDLAI_scl_gut , SHI = IDSHI_scl_gut, CAR = IDCAR_scl_gut, SLO = IDSLO_scl_gut, LOM = IDLOM_scl_gut, GLA = IDGLA_scl_gut, TAR =IDTAR_scl_gut)
write.table(iso_gut_diversity, "~/Desktop/1-s2.0-S1470160X15001569-mmc1/isotopic_gut_diversity.txt", sep = "\t", quote = FALSE)



###  stable isotope ~ morphology ----

# get isotope data 
setwd("~/Desktop/1-s2.0-S1470160X15001569-mmc1")

iso_muscle <- read.table("isotopic_diversity.txt", header = TRUE)
iso_stomach <- read.table("isotopic_gut_diversity.txt", header = TRUE)

iso_muscle2 <- iso_muscle %>% 
  column_to_rownames("Index") 
iso_muscle2 <- data.frame(t(iso_muscle2))
iso_muscle2 <- iso_muscle2 %>% rownames_to_column("Lake")

iso_stomach2 <- iso_stomach %>% 
  column_to_rownames("Index") 
iso_stomach2 <- data.frame(t(iso_stomach2))


summary(lm(iso_stomach2$IPos_d15N ~ iso_muscle2$IPos_d15N))


df_d15N <- data.frame(iso_muscle2$Lake, iso_muscle2$IPos_d15N, iso_stomach2$IPos_d15N)
colnames(df_d15N) <- c("lake", "muscle","gut")

df_d15N$lake <- recode(df_d15N$lake, LAI = "Allt na lairige", SHI = "Shira", CAR = "Carron", SLO = "Sloy", LOM = "Lomond",
                       GLA = "Glashan", TAR = "Tarsan")
df_d15N$lake <- factor(df_d15N$lake, levels = c("Glashan","Tarsan","Allt na lairige", "Shira", "Carron", "Sloy", "Lomond"))

ggplot(df_d15N, aes(x = muscle, y = gut, colour = lake)) + theme.pca() + 
  geom_smooth(method = "lm", colour = "black", fill = "#FFF68F", alpha = 0.3) + 
  geom_point(size = 10) + 
  ylim(0,1) + xlim(0,0.8) + labs(x = expression(paste("muscle ", delta^{15}*N)), y = expression(paste("stomach content ", delta^{15}*N))) + 
  scale_color_manual(values = cols)