library(adegenet) library(ape) rm(list=ls()) ### set the working directory setwd("/your/working/directory") ### read in the data in PLINK format. #snp_data <- read.PLINK("Para.1.raw") #snp_data <- read.PLINK("Para.2.raw") #snp_data <- read.PLINK("Para.3.raw") snp_data ### assign pops to the data geog_grps <- c(1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,7,7,7,7,7,7) snp_geogr_grps <- as.factor(geog_grps) ### select number of clusters (answer the prompt with K for the lowest BIC) clust_snp_data <- find.clusters(snp_data, n.pca=15, max.n.clust=10) ### assign clusters snp_clusters <- as.factor(clust_snp_data$grp) snp_clusters ### choosing PCs for DAPC no_pcs <- dapc(snp_data, n.pc = 25, n.da = 6) a_score <- optim.a.score(no_pcs) ### final DAPC - # PC from above, n.da = number of groups-1 final.dapc.para <- dapc(snp_data, clust_snp_data$grp, n.pca = 1, n.da = 6) scatter(final.dapc.para, scree.pca = TRUE) ### visualize clusters as assigments # 2 clusters scatter(final.dapc.para, scree.pca = TRUE, leg=TRUE, txt.leg=paste("group", c('A', 'B')), col=c("red", "blue")) # 3 clusters scatter(final.dapc.para, scree.pca = TRUE, leg=TRUE, txt.leg=paste("group", c('A', 'B', 'C')), col=c("red", "blue", "forest green")) ### visualize clusters as populations scatter(final.dapc.para, grp = snp_geogr_grps, scree.pca = TRUE, leg=TRUE, txt.leg=paste("group", c('1', '2', '3', '4', '5', '6', '7')), col=c("red", "blue", "forest green", "orange", "purple", "yellow", "pink")) ### make an assignment plot # 2 clusters compoplot(final.dapc.para, only.grp=NULL, subset=NULL, new.pred=NULL, col=c("red", "blue"), lab=NULL, legend=TRUE, txt.leg=NULL, ncol=4, posi=NULL, cex.axis=.75, cleg=.5, bg=transp("white")) # 3 clusters compoplot(final.dapc.para, only.grp=NULL, subset=NULL, new.pred=NULL, col=c("red", "blue", "forest green"), lab=NULL, legend=TRUE, txt.leg=NULL, ncol=4, posi=NULL, cex.axis=.75, cleg=.5, bg=transp("white")) # 4 clusters compoplot(final.dapc.para, only.grp=NULL, subset=NULL, new.pred=NULL, col=c("red", "blue", "forest green", "yellow"), lab=NULL, legend=TRUE, txt.leg=NULL, ncol=4, posi=NULL, cex.axis=.75, cleg=.5, bg=transp("white"))