# script for inbreeding coefficient from plink

# load package ----
library(ggplot2); library(tidyverse);library(ggsignif)

# set theme ----

theme.plot <- function() {
  theme_bw() +
    theme(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_blank(),
          legend.text = element_blank(),
          legend.position = "none")
}

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

# Lomond ----

# setwd 
setwd("~/Desktop/ddRAD_epiRAD_powan/demultiplexed_reads/EpiRAD/Intermediate_files3/06.Population_genetics/lomond/inbreeding/")

# read data 
my_data <- read.csv("inbreeding.LD.maf.csv", header= TRUE)

my_data$Lake <- recode(my_data$Lake, Allt_na_lairige = "Allt na lairige")
my_data$Lake <- factor(my_data$Lake, levels = c( "Lomond", "Allt na lairige", "Shira","Carron", "Sloy"))  # reorder the levels of the lakes

# plot 
pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/lom_inbreeding_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(my_data, aes(x = Lake, y = F, colour = Lake)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.1, size = 6, alpha=0.85) + 
  scale_color_manual(values = cols) + theme.plot()  + ylim(-0.26,0.2) +
  geom_signif(comparisons = list(c("Lomond","Carron")),map_signif_level = TRUE, colour = "black", y_position = -0.2, textsize = 10, tip_length = -0.03) +
  geom_signif(comparisons = list(c("Lomond","Sloy")),map_signif_level = TRUE, colour = "black", y_position = -0.25, textsize = 10, tip_length = -0.03) +
  geom_signif(comparisons = list(c("Lomond","Shira")),map_signif_level = TRUE, colour = "black", y_position = 0.1, textsize = 10) +
  geom_signif(comparisons = list(c("Lomond","Allt na lairige")),map_signif_level = TRUE, colour = "black", y_position = 0.15, textsize = 10) 
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/lom_inbreeding_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(my_data, aes(x = Lake, y = F, fill = Lake)) + geom_violin() + geom_boxplot(width=0.1,outlier.shape = NA) + 
  scale_fill_manual(values = cols) + theme.plot()  + ylim(-0.26,0.2) +
  geom_signif(comparisons = list(c("Lomond","Carron")),map_signif_level = TRUE, colour = "black", y_position = -0.2, textsize = 10, tip_length = -0.03) +
  geom_signif(comparisons = list(c("Lomond","Sloy")),map_signif_level = TRUE, colour = "black", y_position = -0.25, textsize = 10, tip_length = -0.03) +
  geom_signif(comparisons = list(c("Lomond","Shira")),map_signif_level = TRUE, colour = "black", y_position = 0.1, textsize = 10) +
  geom_signif(comparisons = list(c("Lomond","Allt na lairige")),map_signif_level = TRUE, colour = "black", y_position = 0.15, textsize = 10)
dev.off()

# Kolmogorov-Smirnov Tests

ks.test(my_data$F[my_data$Lake == "Lomond"], my_data$F[my_data$Lake == "Carron"])  # D = 0.69318, p-value = 2.535e-05
ks.test(my_data$F[my_data$Lake == "Lomond"], my_data$F[my_data$Lake == "Sloy"])  # D = 0.64577, p-value = 9.326e-07
ks.test(my_data$F[my_data$Lake == "Lomond"], my_data$F[my_data$Lake == "Shira"])  # D = 0.52273, p-value = 0.008305
ks.test(my_data$F[my_data$Lake == "Lomond"], my_data$F[my_data$Lake == "Allt na lairige"])  # D = 0.46591, p-value = 0.1058

mean(my_data$F[my_data$Lake == "Lomond"])
mean(my_data$F[my_data$Lake == "Allt na lairige"])
median(my_data$F[my_data$Lake == "Lomond"])
median(my_data$F[my_data$Lake == "Allt na lairige"])



# Eck ----

# setwd 
setwd("~/Desktop/ddRAD_epiRAD_powan/demultiplexed_reads/EpiRAD/Intermediate_files3/06.Population_genetics/eck/inbreeding/")

# read data
my_data <- read.csv("inbreeding.LD.maf.csv", header= TRUE)

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/eck_inbreeding_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(my_data, aes(x = Lake, y = F, colour = Lake)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.1, size = 6, alpha=0.85) + 
  scale_color_manual(values = cols) + theme.plot() +
  geom_signif(comparisons = list(c("Eck","Tarsan")),map_signif_level = TRUE, colour = "black", y_position = 0.15, textsize = 5) +
  geom_signif(comparisons = list(c("Eck","Glashan")),map_signif_level = TRUE, colour = "black", y_position = 0.1, textsize = 5) 
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/eck_inbreeding_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(my_data, aes(x = Lake, y = F, fill = Lake)) + geom_violin() + geom_boxplot(width=0.1,outlier.shape = NA) + 
  scale_fill_manual(values = cols) + theme.plot() +
  geom_signif(comparisons = list(c("Eck","Tarsan")),map_signif_level = TRUE, colour = "black", y_position = 0.15, textsize = 5) +
  geom_signif(comparisons = list(c("Eck","Glashan")),map_signif_level = TRUE, colour = "black", y_position = 0.1, textsize = 5) 
dev.off()

# Kolmogorov-Smirnov Tests

ks.test(my_data$F[my_data$Lake == "Eck"], my_data$F[my_data$Lake == "Glashan"])  # D = 0.18182, p-value = 0.7338
ks.test(my_data$F[my_data$Lake == "Eck"], my_data$F[my_data$Lake == "Tarsan"])  # D = 0.28977, p-value = 0.229