library(tidyverse)
library(data.table)
library(vegan)
library(qvalue)
library(psych)
library(adegenet)
library(ComplexHeatmap)

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")


# Load pca data ----

epi_combined_pca <- read.table("../datasets/epirad_combined_pca.txt", header = TRUE)

gen_combined_pca <- read.csv("~/Desktop/ddRAD_epiRAD_powan/demultiplexed_reads/EpiRAD/Intermediate_files3/06.Population_genetics/combined/PCA/gen_pca_combined.csv", header = TRUE)

morpho_combined_pca <- read.csv("~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Morphometrics_ecology/morphology_combined_pca.csv", header = TRUE)


morpho_combined_pca2 <- morpho_combined_pca %>% filter(ID %in% gen_combined_pca$sample.id)

epi_combined_pca2 <- epi_combined_pca %>% filter(ID %in% morpho_combined_pca2$ID)

gen_combined_pca2 <- gen_combined_pca %>% filter(sample.id %in% epi_combined_pca2$ID)

# modify the dataframes
morpho_combined_pca2 <- morpho_combined_pca2[c(21:60,1:20,69:97,61:68),]
#
#gen_combined_pca3 <- gen_combined_pca2 %>% filter(sample.id %in% morpho_combined_pca2$ID)
#epi_combined_pca3 <- epi_combined_pca2 %>% filter(ID %in% morpho_combined_pca2$ID)

epi_combined_pca2$System <- morpho_combined_pca2$System
gen_combined_pca2$System <- morpho_combined_pca2$System

epi_combined_pca2$Lake <- recode(epi_combined_pca2$Lake, Allt_na_lairige = "Allt na lairige")
gen_combined_pca2$Lake <- recode(gen_combined_pca2$Lake, Allt_na_lairige = "Allt na lairige")

gen_combined_pca2$Lake <- factor(gen_combined_pca2$Lake, levels = c("Glashan", "Tarsan","Eck","Allt na lairige", "Shira","Carron","Sloy","Lomond"))
epi_combined_pca2$Lake <- factor(epi_combined_pca2$Lake, levels = c("Glashan", "Tarsan","Eck","Allt na lairige", "Shira","Carron","Sloy","Lomond"))
morpho_combined_pca2$Lake <- factor(morpho_combined_pca2$Lake, levels = c("Glashan", "Tarsan","Eck","Allt na lairige", "Shira","Carron","Sloy","Lomond"))



morpho_pca <- ggplot(morpho_combined_pca2, aes(x = PC1, y = PC2,col = Lake, shape = System)) + geom_point(size=6, alpha=0.85) + scale_color_manual(values = cols) +
  theme.pca() +
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') 

epi_pca <- ggplot(epi_combined_pca2, aes(x = PC1, y = PC2,col = Lake, shape = System)) + geom_point(size=6, alpha=0.85) + scale_color_manual(values = cols) +
  theme.pca() +
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') 

gen_pca <- ggplot(gen_combined_pca2, aes(x = EV1, y = EV2,col = Lake, shape = System)) + geom_point(size=6, alpha=0.85) + scale_color_manual(values = cols) +
  theme.pca() +
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') 



info_data <- data.frame(gen_combined_pca2[,1:5],epi_combined_pca2[,2:5], gen_combined_pca2[,7:8])
info_data$Lake <- factor(info_data$Lake,levels = c("Glashan", "Tarsan","Eck","Allt na lairige", "Shira","Carron","Sloy","Lomond"))

# RDA full dataset
# Genomic and epigenomic effects
rda.morpho.all <- rda(morpho_combined_pca2[,2:4] ~ PC1 + PC2 + PC3 + EV1 + EV2 + EV3, data = info_data, scale = TRUE)
RsquareAdj(rda.morpho.all)
signif.all <- anova.cca(rda.morpho.all, parallel = 4, permutations = 999) # default is permutation=999
signif.all

# genomic effects
rda.morpho.gen <- rda(morpho_combined_pca2[,2:4] ~ EV1 + EV2 + EV3, data = info_data, scale = TRUE)
signif.gen <- anova.cca(rda.morpho.gen, parallel = 4, permutations = 999) # default is permutation=999
signif.gen

# epigenomic effects
rda.morpho.epi<- rda(morpho_combined_pca2[,2:4] ~   PC1 + PC2 + PC3, data = info_data, scale = TRUE)
signif.epi <- anova.cca(rda.morpho.epi, parallel = 4, permutations = 999) # default is permutation=999
signif.epi

# genomic effect while accounting for epigenomic
rda.morpho.gen.epi <- rda(morpho_combined_pca2[,2:4] ~   EV1 + EV2 + EV3 + Condition(PC1) + Condition(PC2) + Condition(PC3), data = info_data, scale = TRUE)
RsquareAdj(rda.morpho.gen.epi)
signif.gen.epi <- anova.cca(rda.morpho.gen.epi, parallel = 4, permutations = 999) # default is permutation=999
signif.gen.epi

# epigenomic effect while accounting for genomic
rda.morpho.epi.gen <- rda(morpho_combined_pca2[,2:4] ~   PC1 + PC2 + PC3 + Condition(EV1) + Condition(EV2) + Condition(EV3), data = info_data, scale = TRUE)
RsquareAdj(rda.morpho.epi.gen)
signif.epi.gen <- anova.cca(rda.morpho.epi.gen, parallel = 4, permutations = 999) # default is permutation=999
signif.epi.gen

varpart(morpho_combined_pca2[,2:4], ~ EV1 + EV2 + EV3, ~ PC1+PC2+PC3, data = info_data, scale = TRUE)

## Divide in ECK and LOMOND datasets

eck_data <- filter(info_data, System == "Eck")
lom_data <- filter(info_data, System == "Lomond")
eck_morpho <- filter(morpho_combined_pca2, System == "Eck")
lom_morpho <- filter(morpho_combined_pca2, System == "Lomond")

# eck
plot(eck_data$PC1, eck_data$PC1, col = eck_data$Lake)

rda.eck.morpho.all <- rda(eck_morpho[,2:4] ~ PC1 + PC2 + PC3 + EV1 + EV2 + EV3, data = eck_data, scale = TRUE)
RsquareAdj(rda.eck.morpho.all)
signif.eck.morpho.all <- anova.cca(rda.eck.morpho.all, parallel = 4, permutations = 999) # default is permutation=999
signif.eck.morpho.all

rda.eck.morpho.gen<- rda(eck_morpho[,2:4] ~ (EV1 + EV2 + EV3), data = eck_data, scale = TRUE)
(signif.eck.morpho.gen <- anova.cca(rda.eck.morpho.gen, parallel = 4, permutations = 999)) # default is permutation=999

rda.eck.morpho.epi <- rda(eck_morpho[,2:4] ~ (PC1 + PC2 + PC3), data = eck_data, scale = TRUE)
(signif.eck.morpho.epi <- anova.cca(rda.eck.morpho.epi, parallel = 4, permutations = 999)) # default is permutation=999

rda.eck.morpho.gen.epi <- rda(eck_morpho[,2:4] ~ EV1 + EV2 + EV3 + Condition(PC1) + Condition(PC2) + Condition(PC3), data = eck_data, scale = TRUE)
RsquareAdj(rda.eck.morpho.gen.epi)
signif.eck.morpho.gen.epi <- anova.cca(rda.eck.morpho.gen.epi, parallel = 4, permutations = 999) # default is permutation=999
signif.eck.morpho.gen.epi

rda.eck.morpho.epi.gen <- rda(eck_morpho[,2:4] ~ PC1 + PC2 + PC3 + Condition(EV1) + Condition(EV2) + Condition(EV3), data = eck_data, scale = TRUE)
RsquareAdj(rda.eck.morpho.epi.gen)
signif.eck.morpho.epi.gen <- anova.cca(rda.eck.morpho.epi.gen, parallel = 4, permutations = 999) # default is permutation=999
signif.eck.morpho.epi.gen

varpart(eck_morpho[,2:4], ~ EV1 + EV2 + EV3, ~ PC1+PC2+PC3, data = eck_data, scale = TRUE)

# lomond
plot(lom_data$EV1, lom_data$EV2, col = lom_data$Lake)

# divide lomond in old and young translocated populations
lom_data1 <- filter(lom_data, Lake != "Sloy" & Lake != "Carron")
lom_data2 <- filter(lom_data, Lake == "Sloy" | Lake == "Carron" | Lake == "Lomond")
lom_morpho1 <- filter(lom_morpho, Lake != "Sloy" & Lake != "Carron")
lom_morpho2 <- filter(lom_morpho, Lake == "Sloy" | Lake == "Carron" | Lake == "Lomond")

# lomond and young translocated populations
rda.lom.morpho.all <- rda(lom_morpho1[,2:4] ~ PC1 + PC2 + PC3 + EV1 + EV2 + EV3 , data = lom_data1, scale = TRUE)
RsquareAdj(rda.lom.morpho.all)
signif.lom.morpho.all <- anova.cca(rda.lom.morpho.all, parallel = 4, permutations = 999) # default is permutation=999
signif.lom.morpho.all

rda.lom.morpho.gen <- rda(lom_morpho1[,2:4] ~ EV1 + EV2 + EV3, data = lom_data1, scale = TRUE)
(signif.lom.morpho.gen <- anova.cca(rda.lom.morpho.gen, parallel = 4, permutations = 999)) # default is permutation=999

rda.lom.morpho.epi <- rda(lom_morpho1[,2:4] ~ PC1 + PC2 + PC3, data = lom_data1, scale = TRUE)
(signif.lom.morpho.epi <- anova.cca(rda.lom.morpho.epi, parallel = 4, permutations = 999)) # default is permutation=999

rda.lom.morpho.gen.epi <- rda(lom_morpho1[,2:4] ~ EV1 + EV2 + EV3 + Condition(PC1) + Condition(PC2) + Condition(PC3), data = lom_data1, scale = TRUE)
RsquareAdj(rda.lom.morpho.gen.epi)
signif.lom.morpho.gen.epi <- anova.cca(rda.lom.morpho.gen.epi, parallel = 4, permutations = 999) # default is permutation=999
signif.lom.morpho.gen.epi

rda.lom.morpho.epi.gen <- rda(lom_morpho1[,2:4] ~ PC1 + PC2 + PC3 + Condition(EV1) + Condition(EV2) + Condition(EV3), data = lom_data1, scale = TRUE)
RsquareAdj(rda.lom.morpho.epi.gen)
signif.lom.morpho.epi.gen <- anova.cca(rda.lom.morpho.epi.gen, parallel = 4, permutations = 999) # default is permutation=999
signif.lom.morpho.epi.gen

varpart(lom_morpho1[,2:4], ~ EV1 + EV2 + EV3, ~ PC1+PC2+PC3, data = lom_data1, scale = TRUE)

# lomond and old translocated populations

rda.lom.morpho.all <- rda(lom_morpho2[,2:4] ~ PC1 + PC2 + PC3 + EV1 + EV2 + EV3, data = lom_data2, scale = TRUE)
RsquareAdj(rda.lom.morpho.all)
signif.lom.morpho.all <- anova.cca(rda.lom.morpho.all, parallel = 4, permutations = 999) # default is permutation=999
signif.lom.morpho.all

rda.lom.morpho.gen <- rda(lom_morpho2[,2:4] ~ (EV1 + EV2 + EV3), data = lom_data2, scale = TRUE)
signif.lom.morpho.gen <- anova.cca(rda.lom.morpho.gen, parallel = 4, permutations = 999) # default is permutation=999

rda.lom.morpho.epi <- rda(lom_morpho2[,2:4] ~ (PC1 + PC2 + PC3), data = lom_data2, scale = TRUE)
signif.lom.morpho.epi <- anova.cca(rda.lom.morpho.epi, parallel = 4, permutations = 999) # default is permutation=999

rda.lom.morpho.gen.epi <- rda(lom_morpho2[,2:4] ~ EV1 + EV2 + EV3 + Condition(PC1) + Condition(PC2) + Condition(PC3), data = lom_data2, scale = TRUE)
RsquareAdj(rda.lom.morpho.gen.epi)
signif.lom.morpho.gen.epi <- anova.cca(rda.lom.morpho.gen.epi, parallel = 4, permutations = 999) # default is permutation=999
signif.lom.morpho.gen.epi

rda.lom.morpho.epi.gen <- rda(lom_morpho2[,2:4] ~ PC1 + PC2 + PC3 + Condition(EV1) + Condition(EV2) + Condition(EV3), data = lom_data2, scale = TRUE)
RsquareAdj(rda.lom.morpho.epi.gen)
signif.lom.morpho.epi.gen <- anova.cca(rda.lom.morpho.epi.gen, parallel = 4, permutations = 999) # default is permutation=999
signif.lom.morpho.epi.gen

varpart(lom_morpho2[,2:4], ~ EV1 + EV2 + EV3, ~ PC1+PC2+PC3, data = lom_data2, scale = TRUE)