load(".Rdata")
library("reshape2")
library("ape")
library("diveRsity")
library("dplyr")
library("concaveman")
library(adegenet)
library(ggplot2)
library(hierfstat)
library("PopGenReport")
library(extrafont)
library(extrafontdb)
loadfonts()
library(rworldmap)
library("ggforce")
library("poppr")
library("scatterpie")
library("ggmap") 
library("pegas")
library("related")
library("ggthemes")
library("pophelper")
library("tess3r")
library("ShannonGen")
library("pals")
## Wprowadzenie danych -----------------------
setwd("D:/R/Dzik/")
Dzik <- read.structure("Last_N_pops.stru", n.ind = 169, n.loc = 12, onerowperind = T,
                       col.lab = 1, col.pop = 3, col.others = 2, row.marknames = 1)
Dzik@other$latlong <- read.table("coords.txt", sep = "\t")
colnames(Dzik@other$latlong) <- c("lat","long")
N_map <- ggmap::get_stamenmap(bbox = c(left = (min(Dzik@other$latlong$long)-0.02), 
                                       bottom = (min(Dzik@other$latlong$lat)-0.02),
                                       right = (max(Dzik@other$latlong$long)+0.02), 
                                       top = (max(Dzik@other$latlong$lat)+0.02)), 
                              force = TRUE,zoom = 11)

data <- data.frame(ID=indNames(Dzik),pop=Dzik$pop,Dzik@other$latlong)
head(data)
pdf("Populations.pdf",width = 10,height = 6)
ggmap::ggmap(N_map) +
  scale_fill_manual(values = as.character(glasbey(length(levels(data$pop))))) +
  geom_point(data = data, shape = 21, size = 3, stroke = 0.2, 
             aes(long, lat, fill = factor(pop))) +
  labs(fill = "pop") +
  coord_equal()
dev.off()

strata(Dzik) <- data.frame(pop = pop(Dzik)) #ustawianie strata

TableOfVariance <- data.frame(N = sapply(seppop(Dzik),nInd))
diversities <- basic.stats(Dzik)
TableOfVariance <- cbind(TableOfVariance,Ho = apply(diversities$Ho,2,mean))
TableOfVariance <- cbind(TableOfVariance,t(apply(diversities$Ho,2,boot.div)))
TableOfVariance <- cbind(TableOfVariance,Ho = apply(diversities$Hs,2,mean))
TableOfVariance <- cbind(TableOfVariance,t(apply(diversities$Hs,2,boot.div)))
ARSes <- hierfstat::allelic.richness(Dzik)$Ar
TableOfVariance <- cbind(TableOfVariance,Rs = apply(ARSes,2,mean))
TableOfVariance <- cbind(TableOfVariance,t(apply(ARSes,2,boot.div)))
Shannons <- ShannonGen(Dzik)[[1]]
TableOfVariance <- cbind(TableOfVariance,H = apply(Shannons,2,mean))
TableOfVariance <- cbind(TableOfVariance,t(apply(Shannons,2,boot.div)))
TableOfVariance <- cbind(TableOfVariance,Fis = apply(basic.stats(Dzik)$Fis,2,mean))
TableOfVariance  <- cbind(TableOfVariance,boot.ppfis(genind2hierfstat(Dzik),9999,dig=3)[[2]])

write.table(TableOfVariance,file = "BasicStats.xls",sep = "\t",quote = FALSE,col.names = NA)


### Null Alleles by popgenreport
Dzik.PGR.null <- popgenreport(Dzik, mk.null.all =  TRUE,path.pgr = "D:/R/Dzik/Nullalleles", 
                              foldername = "Total",mk.pdf = TRUE)
for (pop in levels(Dzik$pop)){
popgenreport(Dzik[Dzik$pop==pop], mk.counts = FALSE, mk.null.all =  TRUE,path.pgr = "D:/R/Dzik/Nullalleles", 
             foldername = pop,mk.pdf = TRUE)
}
# Hardy Weinberg expectations

cat(file="HWeq.xls")
for (pop in popNames(Dzik)){
  popul <- Dzik[pop(Dzik)==pop]
  cat(paste(pop,nInd(popul),"\n"),file="HWeq.xls", sep = "\t",append = TRUE)
  write.table(pegas::hw.test(popul, B = 9999),file="HWeq.xls", 
              quote = FALSE, sep = "\t",append = TRUE,col.names = NA)
}
boot.div <- function(data,n.boot = 99999){
  boots <- sapply(1:n.boot,function(i) mean(sample(data,length(data),replace = TRUE),na.rm = TRUE))
  return(c(quantile(boots,0.025),  quantile(boots,0.975)))
}
sHannon <- ShannonGen(Dzik)[[1]]
apply(sHannon,2,boot.div)

write.table(poppr.amova(Dzik, ~pop)$componentsofcovariance,file = "AMOVA.txt",sep = "\t",quote = FALSE)
Dzik_genpop <- genind2genpop(Dzik)
Dzik_dist <- as.matrix(dist(Dzik_genpop))
genotypy <- genind2df(Dzik)
genotypy$pop <- as.numeric(Dzik$pop)
genotypy[is.na(genotypy)] <- "0"
Dzik.matFst <- pairwise.WCfst(genind2hierfstat(Dzik))
Dzik.matFst_melted <- reshape2::melt(Dzik.matFst, na.rm = TRUE)
Dzik.barplot.Fst <- Dzik.matFst_melted[Dzik.matFst_melted$Var1!=Dzik.matFst_melted$Var2,]
pdf("Fig_2_Fst.pdf", width = 12, height = 9)
boxplot(value~Var1, data = Dzik.barplot.Fst,
        las=1,  ylab="Fst",horizontal=TRUE,outline=FALSE)
dev.off()

boots <- hierfstat::boot.ppfst(genind2hierfstat(Dzik),nboot=999)

pdf("FST heat.pdf", width = 5,height = 5)
ggplot(data = Dzik.matFst_melted, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "white", high = "black",  
                       midpoint = 0.012, limit = c(0,0.08), space = "Lab") +
  theme_minimal()+ 
  geom_label(aes(Var2, Var1, label = value), size = 4, 
             fill = "#EEEEEE", label.size = 0, label.padding = unit(0.1, "lines")) + #to jest do porz?dkowania wykresu
  theme(
    axis.text.x = element_text(vjust = 1, size = 12, hjust = 0.5),
    axis.text.y = element_text(vjust = 0.5, size = 12, hjust = 0),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.major = element_blank(),
    panel.border = element_blank(),
    panel.background = element_blank(),
    axis.ticks = element_blank())
dev.off()

pops <- levels(pop(Dzik))
fstTable <- pairwise.WCfst(genind2hierfstat(Dzik))
nperm <- 10000
for (pop1 in 1:(length(pops)-1)){
  for (pop2 in (pop1+1):length(pops)){
    newDzik <- repool(Dzik[pop(Dzik)==pops[pop1]],Dzik[pop(Dzik)==pops[pop2]])
    newDzik <- genind2hierfstat(newDzik)
    permuted <- c() 
    for (perm in 1:nperm){
      newDzik$pop <- sample(newDzik$pop)
      permuted <- append(permuted,pairwise.WCfst(newDzik)[1,2])
    }
    fstTable[pop2,pop1] <- round((sum(permuted>=fstTable[pop1,pop2])+1)/(nperm+1),4)
    cat(paste(pops[pop1],"/",pops[pop2],"- Done!!!\n"))
  }
}

# PCA ---------------------------------------------------------------------
x.Dzik <- tab(Dzik, freq=T, NA.method="mean")

par(mar=c(4.5, 7.5, 1, 1),bg="white") #to s? jakie? parametry rysowania wykres?w

xval<-xvalDapc(x.Dzik, grp = pop(Dzik), n.pca.max = 100, training.set=0.8, result="groupMean", 
               scale=TRUE,n.rep=1000,n.pca=seq(15,30,by=1),xval.plot = T,
               parallel = "snow",ncpus=4)
xval$`Number of PCs Achieving Highest Mean Success`
xval$`Number of PCs Achieving Lowest MSE`

dzik_dapc <- dapc(Dzik, grp = pop(Dzik), scale = FALSE, n.pca = 24)
5

pdf("Fig_5a.pdf",width = 12,height = 8)
scatter(dzik_dapc,col = colorblind_pal()(7),cex =3,leg=TRUE,clab=0,pch=20,solid = 0.9,
        cstar=0, scree.da = FALSE)
dev.off()

pdf("compoplot2.pdf",width = 12,height = 4)
compoplot(dzik_dapc, col = colorblind_pal()(7),show.lab=TRUE,lab = pop(Dzik),
          leg=FALSE,cex.names=0.5)
dev.off()
plot()

#======== Structure visualisation =======
StructureDir <- "d:/R/Dzik/Structure/"
projects <- list.dirs(paste0(StructureDir),full.names =FALSE,recursive = FALSE)
geo <- read.table(file=paste0(StructureDir,"project_data"),header = TRUE,sep = "\t",
                      stringsAsFactors = FALSE, row.names = 1)[c(2,6,7)]
colnames(geo) <- c("Pop","lat","long")
for (i in projects){
  structFiles <- list.files(path=paste0(StructureDir,i,"/Results/"),full.names = TRUE)
  Structs <- pophelper::readQ(files = structFiles, filetype = "structure", 
                              indlabfromfile = TRUE)
  qAligned <- alignK(Structs)
  qTable <- tabulateQ(qAligned)
  qSummary <- summariseQ(qTable)
  write.table(qSummary,file=paste0("Summary_",i,".xls"),sep = "\t",quote = FALSE)
  evan <- evannoMethodStructure(qSummary,returnplot = TRUE, returndata = TRUE, 
                                exportplot = TRUE, outputfilename = paste0("Evanno_",i), imgtype = "pdf",
                                basesize = 6,xaxisbreaks = unique(qTable[,"k"]),
                                exportpath = getwd())
  K = evan$data[which(evan$data[,"deltaK"]==max(evan$data[,"deltaK"],na.rm = TRUE)),"k"]
  qMerged <- mergeQ(qAligned)
  for (j in 1:length(qSummary[,1])){
    attr(qMerged[[j]],"ind") <- qSummary[j,"ind"]
    attr(qMerged[[j]],"k") <- qSummary[j,"k"]
  }
  plotQ(qMerged[K],imgtype = "pdf",outputfilename = paste0("StrBarplot_",i,"_K",K), 
        returnplot = TRUE, width = 30, height =6, grplabsize=3.5,
        showindlab=FALSE, useindlab = TRUE, grplab = geo[1], basesize = 15,
        exportpath = getwd())
}


# Zrzucanie danych ze Structure do pliku
Structured <- data.frame(ID=rownames(geo),pop=geo[,1])
Structured[,3:4] <- data.frame(qMerged$`2`) 
Structured <- Structured[Structured[,2]!="TON",]
Structured$X <- coords[,2]
Structured$Y <- coords[,1]
Structured <- Structured[,c(1,2,5,6,3,4)]
write.table(Structured,file = "Structure_K2.txt",sep = "\t",quote = FALSE,col.names = NA)
write.table(Structured,file = "Structure_K2.csv",sep = ",",quote = FALSE,col.names = NA)


Structured <- data.frame(ID=rownames(geo),pop=geo[,1])
Structured[,3:5] <- data.frame(qMerged$`3`) 
Structured <- Structured[Structured[,2]!="TON",]
Structured$X <- coords[,2]
Structured$Y <- coords[,1]
Structured <- Structured[,c(1,2,6,7,3,4,5)]
write.table(Structured,file = "Structure_K3.txt",sep = "\t",quote = FALSE,col.names = NA)
write.table(Structured,file = "Structure_K3.csv",sep = ",",quote = FALSE,col.names = NA)

# Geneland
x <- read.table("Geneland_K2/proba.pop.membership.indiv.txt")


#Mapka z wynikiem ze Structured - K2 albo K3
write.table(Structured,file = "StructureK2.txt",sep = "\t",quote = TRUE)
write.table(Structured,file = "StructureK2.csv",sep = ",",quote = TRUE)