## plotting relatedness (plink)

# library ----
library(tidyverse);library(reshape2);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.position = "none")
}

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

# Lomond ----

# set working directory 
setwd("~/Desktop/ddRAD_epiRAD_powan/demultiplexed_reads/EpiRAD/Intermediate_files3/06.Population_genetics/lomond/relatedness/")

# read data 
my_matrix_lom <- read.table("lom.LD.maf.rel", header = FALSE)
my_id_lom <- read.table("lom.LD.maf.rel.id", header = FALSE)

colnames(my_matrix_lom) <- my_id_lom$V1
rownames(my_matrix_lom) <- my_id_lom$V2
lom_matrix <- as.matrix(my_matrix_lom)

# manipulate 

lai <- lom_matrix[1:8,1:8]
lai2 <- melt(lai)
lai3 <- filter(lai2, value < 0.5 & value != 0)
lai3$lake <- "Allt na lairige"

shi <- lom_matrix[9:21,9:21]
shi2 <- melt(shi)
shi3 <- filter(shi2, value < 0.5 & value != 0)
shi3$lake <- "Shira"

car <- lom_matrix[22:37,22:37]
car2 <- melt(car)
car3 <- filter(car2, value < 0.5 & value != 0)
car3$lake <- "Carron"

sloy <- lom_matrix[38:66,38:66]
sloy2 <- melt(sloy)
slo3 <- filter(sloy2, value < 0.5 & value != 0)
slo3$lake <- "Sloy"

lom <- lom_matrix[67:110,67:110]
lom2 <- melt(lom)
lom3 <- filter(lom2, value < 0.5 & value != 0)
lom3$lake <- "Lomond"

# combine all values in a dataframe
lom_relatedness <- rbind(lai3,shi3,car3,slo3,lom3)

# plot 
pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/lom_relat_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(lom_relatedness, aes(x = lake, y = value, colour = lake)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.1, size = 3, alpha=0.85) + 
  scale_color_manual(values = cols) + theme.plot() + scale_x_discrete(limits = c("Lomond", "Allt na lairige", "Shira", "Carron", "Sloy")) + 
  labs(x = "Lake", y = "Relatedness")  +
  geom_signif(comparisons = list(c("Lomond","Carron")),map_signif_level = TRUE, colour = "black", y_position = 0.25, textsize = 10) +
  geom_signif(comparisons = list(c("Lomond","Sloy")),map_signif_level = TRUE, colour = "black", y_position = 0.3, textsize = 10) +
  geom_signif(comparisons = list(c("Lomond","Shira")),map_signif_level = TRUE, colour = "black", y_position = 0.2, textsize = 5) +
  geom_signif(comparisons = list(c("Lomond","Allt na lairige")),map_signif_level = TRUE, colour = "black", y_position = 0.15, textsize = 5) 
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/lom_relat_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(lom_relatedness, aes(x = lake, y = value, fill = lake)) + geom_violin() + geom_boxplot(width=0.1,outlier.shape = NA) +
  scale_fill_manual(values = cols) + theme.plot() + scale_x_discrete(limits = c("Lomond", "Allt na lairige", "Shira", "Carron", "Sloy")) + 
  labs(x = "Lake", y = "Relatedness")  +
  geom_signif(comparisons = list(c("Lomond","Carron")),map_signif_level = TRUE, colour = "black", y_position = 0.3, textsize = 10) +
  geom_signif(comparisons = list(c("Lomond","Sloy")),map_signif_level = TRUE, colour = "black", y_position = 0.35, textsize = 10) +
  geom_signif(comparisons = list(c("Lomond","Shira")),map_signif_level = TRUE, colour = "black", y_position = 0.25, textsize = 10) +
  geom_signif(comparisons = list(c("Lomond","Allt na lairige")),map_signif_level = TRUE, colour = "black", y_position = 0.20, textsize = 10)
dev.off()

# Kolmogorov-Smirnov Tests

ks.test(lom_relatedness$value[lom_relatedness$lake == "Lomond"], lom_relatedness$value[lom_relatedness$lake == "Carron"])  # 0.68614, p-value < 2.2e-16
ks.test(lom_relatedness$value[lom_relatedness$lake == "Lomond"], lom_relatedness$value[lom_relatedness$lake == "Sloy"])  # D = 0.63377, p-value < 2.2e-16
ks.test(lom_relatedness$value[lom_relatedness$lake == "Lomond"], lom_relatedness$value[lom_relatedness$lake == "Shira"])  # D = 0.14251, p-value = 0.1051
ks.test(lom_relatedness$value[lom_relatedness$lake == "Lomond"], lom_relatedness$value[lom_relatedness$lake == "Allt na lairige"])  # D = 0.18001, p-value = 0.3416



# Eck ----

# set working directory 
setwd("~/Desktop/ddRAD_epiRAD_powan/demultiplexed_reads/EpiRAD/Intermediate_files3/06.Population_genetics/eck/relatedness/")

# read data 

my_matrix_eck <- read.table("eck.LD.maf.rel", header = FALSE)
my_id_eck <- read.table("eck.LD.maf.rel.id", header = FALSE)

colnames(my_matrix_eck) <- my_id_eck$V1
rownames(my_matrix_eck) <- my_id_eck$V1
eck_matrix <- as.matrix(my_matrix_eck)

# manipulate 

gla <- eck_matrix[1:17,1:17]
gla2 <- melt(gla)
gla3 <- filter(gla2, value < 0.5 & value != 0)
gla3$lake <- "Glashan"

tar <- eck_matrix[18:33,18:33]
tar2 <- melt(tar)
tar3 <- filter(tar2, value < 0.5 & value != 0)
tar3$lake <- "Tarsan"

eck <- eck_matrix[34:77,34:77]
eck2 <- melt(eck)
eck3 <- filter(eck2, value < 0.5 & value != 0)
eck3$lake <- "Eck"

# combine all values in a dataframe
eck_relatedness <- rbind(gla3,tar3,eck3)

# plot
pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/eck_relat_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(eck_relatedness, aes(x = lake, y = value, colour = lake)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.1, size = 3, alpha=0.85) + 
  scale_color_manual(values = cols) + theme.plot() + scale_x_discrete(limits = c("Eck", "Glashan", "Tarsan")) + 
  labs(x = "Lake", y = "Relatedness") +
  geom_signif(comparisons = list(c("Eck","Tarsan")),map_signif_level = TRUE, colour = "black", y_position = 0.2, textsize = 10) +
  geom_signif(comparisons = list(c("Eck","Glashan")),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/eck_relat_maf_stacks.pdf", width = 11.69, height = 8.27)
ggplot(eck_relatedness, aes(x = lake, y = value, fill = lake)) + geom_violin() + geom_boxplot(width=0.1,outlier.shape = NA) + 
  scale_fill_manual(values = cols) + theme.plot() + scale_x_discrete(limits = c("Eck", "Glashan", "Tarsan")) + 
  labs(x = "Lake", y = "Relatedness") +
  geom_signif(comparisons = list(c("Eck","Tarsan")),map_signif_level = TRUE, colour = "black", y_position = 0.2, textsize = 10) +
  geom_signif(comparisons = list(c("Eck","Glashan")),map_signif_level = TRUE, colour = "black", y_position = 0.15, textsize = 10) 
dev.off()

# Kolmogorov-Smirnov Tests

ks.test(eck_relatedness$value[eck_relatedness$lake == "Eck"], eck_relatedness$value[eck_relatedness$lake == "Glashan"])  # D = 0.10011, p-value = 0.1843
ks.test(eck_relatedness$value[eck_relatedness$lake == "Eck"], eck_relatedness$value[eck_relatedness$lake == "Tarsan"])  # D = 0.15298, p-value = 0.01369