library(data.table);library(pcaMethods);library(tidyverse)

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


read_counts <- read.table("rlog_reads_combined.txt", header = TRUE, row.names = 1)
lom_counts <- read.table("rlog_lom_reads.txt", header = TRUE, row.names = 1)
eck_counts <- read.table("rlog_eck_reads.txt", header = TRUE, row.names = 1)
lake_data1 <- read.csv("../RDA/lakes_annotation.csv", header = TRUE)
#lake_data2 <- read.csv("lakes_env.csv", header = TRUE)

read_counts_t <- t(read_counts)
lom_counts_t <- t(lom_counts)
eck_counts_t <- t(eck_counts)

respca <- pca(read_counts_t, method = "svd", center = TRUE, nPcs = 5)
lom_pca <- pca(lom_counts_t, method = "svd", center = TRUE, nPcs = 5)
eck_pca <- pca(eck_counts_t, method = "svd", center = TRUE, nPcs = 5)

combined_data <- scores(respca)
combined_data <- data.frame(combined_data)
combined_data$Lake <- lake_data1$Lake
combined_data$Lake <- recode(combined_data$Lake, Allt_na_lairige = "Allt na lairige")
combined_data <- combined_data %>% rownames_to_column("ID") 
write.table(combined_data, "epirad_combined_pca.txt", quote = FALSE, row.names = FALSE, sep = "\t")

lom_data <- scores(lom_pca)
lom_data <- data.frame(lom_data)
lom_data$Lake <- lake_data1$Lake[lake_data1$System == "Lomond"]
lom_data$Lake <- recode(lom_data$Lake, Allt_na_lairige = "Allt na lairige")
lom_data <- lom_data %>% rownames_to_column("ID") 
write.table(lom_data, "epirad_lom_pca.txt", quote = FALSE, row.names = FALSE, sep = "\t")

eck_data <- scores(eck_pca)
eck_data <- data.frame(eck_data)
eck_data$Lake <- lake_data1$Lake[lake_data1$System == "Eck"]
eck_data <- eck_data %>% rownames_to_column("ID") 
write.table(eck_data, "epirad_eck_pca.txt", quote = FALSE, row.names = FALSE, sep = "\t")


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

combined_data$Lake <- factor(combined_data$Lake,levels = c("Glashan", "Tarsan","Eck","Allt na lairige", "Carron","Shira","Sloy","Lomond"))
lom_data$Lake <- factor(lom_data$Lake,levels = c("Allt na lairige", "Carron","Shira","Sloy","Lomond"))
eck_data$Lake <- factor(eck_data$Lake,levels = c("Glashan", "Tarsan","Eck"))

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/epira_combined_pc12.pdf", width = 11.69, height = 8.27)
ggplot(combined_data, aes(x = PC1, y = PC2, colour = Lake)) + theme.pca() + geom_point(size = 6, alpha = 0.85) + scale_colour_manual(values = cols) + 
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') + labs(x = "PC1 5%", y = "PC2 2%")
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/epirad_lom_pc12.pdf", width = 11.69, height = 8.27)
ggplot(lom_data, aes(x = PC1, y = PC2, colour = Lake)) + theme.pca() + geom_point(size = 6, alpha = 0.85) + scale_colour_manual(values = cols) + 
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') + labs(x = "PC1 2.8%", y = "PC2 2.6%")
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/epirad_lom_pc13.pdf", width = 11.69, height = 8.27)
ggplot(lom_data, aes(x = PC1, y = PC3, colour = Lake)) + theme.pca() + geom_point(size = 6, alpha = 0.85) + scale_colour_manual(values = cols) + 
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') + labs(x = "PC1 2.8%", y = "PC3 2.4%")
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/epirad_eck_pc12.pdf", width = 11.69, height = 8.27)
ggplot(eck_data, aes(x = PC1, y = PC2, colour = Lake)) + theme.pca() + geom_point(size = 6, alpha = 0.85) + scale_colour_manual(values = cols) + 
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') + labs(x = "PC1 4%", y = "PC2 4%") 
dev.off()

pdf(file = "~/Dropbox/Marco_Crotti/Evolutionary genomics of whitefish/Translocation project/Genomic analyses/figures/epirad_eck_pc13.pdf", width = 11.69, height = 8.27)
ggplot(eck_data, aes(x = PC1, y = PC3, colour = Lake)) + theme.pca() + geom_point(size = 6, alpha = 0.85) + scale_colour_manual(values = cols) + 
  geom_vline(xintercept = 0, linetype = 'dashed') + geom_hline(yintercept = 0, linetype = 'dashed') + labs(x = "PC1 4%", y = "PC23 3%") 
dev.off()