setwd("~/Desktop/ddRAD_epiRAD_powan/demultiplexed_reads/EpiRAD/Intermediate_files3/08.EpiRAD/RDA/")

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

# epigenetic RDA ----
combined_rlog_data <- read.table("rlog_reads_combined.txt",header=TRUE,row.names = 1)
combined_rlog_data <- data.frame(t(combined_rlog_data))
lake_data <- read.csv("lakes_annotation.csv",header=TRUE, row.names = 1)

identical(rownames(combined_rlog_data), rownames(lake_data))

lakes.rda1 <- rda(combined_rlog_data ~ Lake_type + System, data = lake_data, scale = TRUE)
lakes.rda1

vif.cca(lakes.rda1)
RsquareAdj(lakes.rda1)
summary(eigenvals(lakes.rda1, model = "constrained"))
screeplot(lakes.rda1)

signif.full <- anova.cca(lakes.rda1, parallel = 4, permutations = 999) # default is permutation=999
signif.full

signif.axis <- anova.cca(lakes.rda1, by="axis", parallel = 4)
signif.axis

lake <- lake_data$Lake
bg <- c("#1B9E77","#1B9E77","#D95F02","#1B9E77","#D95F02","#1B9E77","#1B9E77","#1B9E77")
syst <- lake_data$System
shape <- c(21,24)

plot(lakes.rda1, type="n", scaling=3)
points(lakes.rda1, display="species", pch=20, cex=0.7, col="gray32", scaling=3)           # the SNPs
points(lakes.rda1, display="sites", pch=shape[syst], cex=3, col="gray32", scaling=3, bg=bg[lake]) # the wolves
text(lakes.rda1, scaling=3, display="bp", col="#0868ac", cex=2)                           # the predictors
legend("bottomright", legend=c("refuge","source"), bty="n", fill=c("#1B9E77","#D95F02"), cex=3)

lakes.rda1_scores <- vegan::scores(lakes.rda1, choices = 1:3, display = "species") 
lakes.rda1_scores <- as.data.frame(lakes.rda1_scores)
lakes.rda1_scores$zscore1 <- ((lakes.rda1_scores$RDA1 - (mean(lakes.rda1_scores$RDA1)))/sd(lakes.rda1_scores$RDA1))
lakes.rda1_scores$zscore2 <- ((lakes.rda1_scores$RDA2 - (mean(lakes.rda1_scores$RDA2)))/sd(lakes.rda1_scores$RDA2))

head(lakes.rda1_scores)
write.csv(lakes.rda1_scores, "RDA_loadings_lakes_RDAzscore_methylation.csv")

outliers <- lakes.rda1_scores %>% 
  rownames_to_column('Locus') %>% 
  filter(zscore2 >= 2.5 | zscore2 <= -2.5) 

write.csv(outliers, "epirad_RDA_laketype_outliers.csv", row.names = FALSE)