#R包的安装（1）    install.packages('devtools')devtools::install('royfrancis/pophelper')

#R包的安装（2）install.packages(c('ggplot2','gridExtra','label','switching','tidyr','remotes'),                 repos = "http://cloud.r-project.org")#检查是否成功安装，如果成功安装可报pophelper的版本信息library(pophelper)

#首先设置工作目录，我的structure结果文件被存储在‘F:群体遗传学习2.群体结构分析2.5structure分析’中setwd('F:\群体遗传学习\2.群体结构分析\2.5structure分析')#读取文件，小伙伴们在做structure时可知，其输入文件如下，因此为避免重复

files <- list.files(path = '.', #文件路径                    pattern = '_f$',  #以_f结尾的文件                    full.names = F)  #是否展示路径信息files  #查看files的内容

#Q矩阵阅读slist <- readQ(files = files, #文件               indlabfromfile = T,  #文件中每行的命名               filetype = 'structure')   #文件类型auto', 'structure','tess2','baps','basic' or 'clumpp'

#对齐xlist <- alignK(slist)

#合并，共生成11个列表形式的文件，对应一开始读入的11个k值文件xlist2 <- mergeQ(qlist)#如果报错如下

#需要利用以下函数解决，可直接复制，但正常输出的结果都可用上述mergeQ进行合并，应为alignK继而mergeQ是同一个函数下的功能且应用顺序紧挨着mergy <- function(x) {return(list(Reduce(`+`, x)/length(x)))}fun1 <- function(x) as.matrix(unlist(attributes(x)))a <- lapply(xlist, fun1)b <- as.data.frame(t(as.data.frame(lapply(a, function(x,y) x[y, ], "k"), stringAsFactors = FALSE)), stringsAsFactors = FALSE)
fun2 <- function(x) if (all(!is.na(as.numeric(as.character(x))))) {return(as.numeric(as.character(x)))}else {return(x)}b <- as.data.frame(sapply(b, fun2), stringAsFactors = FALSE)ord <- order(b[, "k"])qlist <- xlist[ord]labels <- summariseQ(tabulateQ(qlist, sorttable = FALSE))$kx <- unlist(lapply(splitQ(qlist), mergy), recursive = FALSE)names(x) <- labelsxlist2 <- as.qlist(x)xlist2

#随后进行图形绘制P1 <- plotQ(xlist2,  #Q矩阵#图片合并,因为我们输入了11个文件所以不能选用sep(单个输出)      imgoutput = "join",           imgtype = "pdf",#输出格式"png","jpeg","tiff" or "pdf".      width = 20,height = 1.2, #宽度和高度      showindlab = T,  #是否对每个K值图形绘制边框      useindlab = T, #展示单株标签      indlabhjust = 0.5,  #单株标签水平居中      indlabangle = 0,     #单株标签角度      indlabsize = 1.2,  #单株标签文本大小      sortind = 'all',   #排序，基于聚类排序#K值图不会共享单株标签名,且只在 imgoutput = "join"时有效      sharedindlab = F,        panelspacer = 0.02,   #不同k值图的间隙      indlabspacer = 0,  #不同k值图中单株空隙      indlabcol = "black",  #k值文本颜色      showtitle = T,titlelab = "",titlecol = "black",titlehjust = 0.5,    #标题      showsp = T,sppos = "left",   #k值标签展示,也有‘right’      splab = paste0("K=",1:11),  #k值标签      splabangle = 180,   #标签角度      dpi = 600,  #输出分辨率      outputfilename = paste0("structure1"),  #输出文件的名称      units = "cm",   #输出文件的大小      exportplot = T,  #图片保存      exportpath = getwd())  #保存路径#届时会在自己设置的工作目录下出现名为structure1.pdf的目标文件

#此外我们还可以进行分组显示#引入分组信息pop <- read.table('group.txt', header = T, sep = 't', stringsAsFactors = F,row.names = 1)pop

##绘图P2 <- plotQ(xlist2,imgoutput = "join",imgtype = "png",  #输出格式的变化height = 1.2,width = 20,showindlab = T,useindlab = T,sharedindlab = T,showsp = T,sppos = "left",splab = paste0("K=",1:11),splabangle = 180,      showlegend = F,ordergrp = T,grplab = pop,   #添加分组      subsetgrp = LETTERS[1:9],  #添加组名      dpi = 600,      outputfilename = paste0("structure_barplot_pop1"),units = "cm",exportpath = getwd())

#此时已经可以在结构图下面看到分组信息
#最后推断最优的分群数目#利用k值及deltaK靓列数确定最优k值，deltaK最大的最优dt <- evannoMethodStructure(summariseQ(tabulateQ(slist,writetable = F,sorttable = T)),writetable = T,exportpath = getwd())

#写出最优Q矩阵,第6个是最优的write.table(x = xlist2[[6]],file ='best_k6_Qscore',quote = F,sep = 't')