### Geneland #### library(Geneland) library(PBSmapping) library(RgoogleMaps) coords <- read.table(file = "Coords.txt",sep = "\t") colnames(coords) <- c("Y","X") genetable <- read.csv(file = "Last_N_pops.stru",header = TRUE,sep="\t") attr(coords, "projection") <- "LL" attr(coords, "zone") <- NA colnames(coords) <- c("Y","X") coord.utm <- PBSmapping::convUL(coords) # converting Long/Lat to UTM write.table(cbind(coords, pop = Dzik@pop), file = "pops.csv", sep = ",",quote = FALSE,col.names = NA) gene <- genetable[,9:32] folder <- "Geneland" MCMC(coordinates = coord.utm, geno.dip.codom = gene, varnpop = TRUE, npopmax = 2, npopmin = 2, spatial = TRUE, freq.model = "Correlated", nit = 10000000, thinning = 1000,path.mcmc = folder) PostProcessChain(coordinates = coord.utm, path.mcmc = paste0(folder,"/"), nxdom = 100, nydom = 100, burnin = 1000) Plotnpop(path.mcmc = paste0(folder,"/"), burnin = 1000) GenlandRes <- cbind(coords,Clust = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3]) write.table(GenlandRes,file = paste0("Genland_opt.csv"), sep = ",",quote = FALSE,col.names = NA) Cluster = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3] Geneland_opt = cbind(coords,Cluster) table(Cluster) pdf("Geneland_K_opt.pdf", width = 10, height = 6) ggmap::ggmap(N_map)+ scale_fill_manual(values = as.character(glasbey(10))) + geom_point(data = data.frame(Geneland_opt),shape = 21, size = 4, aes(X, Y, fill = factor(Cluster))) + labs(fill = "Cluster")+ coord_equal() dev.off() ####### ####K2 folder <- "Geneland_K2" MCMC(coordinates = coord.utm, geno.dip.codom = gene, varnpop = TRUE, npopmax = 2, npopmin = 2, spatial = TRUE, freq.model = "Correlated", nit = 10000000, thinning = 1000,path.mcmc = folder) PostProcessChain(coordinates = coord.utm, path.mcmc = paste0(folder,"/"), nxdom = 100, nydom = 100, burnin = 1000) Plotnpop(path.mcmc = paste0(folder,"/"), burnin = 1000) GenlandRes <- cbind(coords,Clust = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3]) write.table(GenlandRes,file = paste0("Genland_K2.csv"), sep = ",",quote = FALSE,col.names = NA) Cluster = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3] Geneland_opt = cbind(coords,Cluster) pdf("Geneland_K2.pdf", width = 10, height = 6) ggmap::ggmap(N_map)+ scale_fill_manual(values = as.character(glasbey(2))) + geom_point(data = data.frame(Geneland_opt),shape = 21, size = 4, aes(X, Y, fill = factor(Cluster))) + labs(fill = "Cluster")+ coord_equal() dev.off() ###### ####K3 folder <- "Geneland_K3" MCMC(coordinates = coord.utm, geno.dip.codom = gene, varnpop = TRUE, npopmax = 3, npopmin = 3, spatial = TRUE, freq.model = "Correlated", nit = 10000000, thinning = 1000,path.mcmc = folder) PostProcessChain(coordinates = coord.utm, path.mcmc = paste0(folder,"/"), nxdom = 100, nydom = 100, burnin = 1000) Plotnpop(path.mcmc = paste0(folder,"/"), burnin = 1000) GenlandRes <- cbind(coords,Clust = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3]) write.table(GenlandRes,file = paste0("Genland_K3.csv"), sep = ",",quote = FALSE,col.names = NA) Cluster = read.csv(file = paste0(folder,"/modal.pop.indiv.txt"),sep = " ",header = FALSE)[,3] Geneland_opt = cbind(coords,Cluster) pdf("Geneland_K3.pdf", width = 10, height = 6) ggmap::ggmap(N_map)+ scale_fill_manual(values = as.character(glasbey(3))) + geom_point(data = data.frame(Geneland_opt),shape = 21, size = 4, aes(X, Y, fill = factor(Cluster))) + labs(fill = "Cluster")+ coord_equal() dev.off()