#使用方法：Rscript hap_fig.R genename
library(agricolae)
library(gridExtra)
library(ggplot2)
args = commandArgs(T)
genename=args[1]  #基因名称
md <- read.table(paste0('./04.R_data/',genename,"_R_data.xls"),header=T) #读取文件
colnames=names(md) #将表头写入colnames
plot_list = list()
#批量处理文件
for(i in 3:length(colnames)){
  real= colnames[i]
  #Hap= colnames[1]
  data_i <- md[,c(1,i)] #将作图数据写入data_i
  data_i <- na.omit(data_i) #？
  #将数据表头改为"Hap trait"
  colnames(data_i)=c("Hap","trait")
  #计算多重比较
  mod = aov(data_i$trait ~ data_i$Hap, data=data_i)
  re = LSD.test(mod,"data_i$Hap",alpha = 0.05)
  mean<- re$means[,1]
  mar<- re$groups
  rownamemar <- row.names(mar)
  newmar<-data.frame(rownamemar,mar)
  sort<- newmar[order(newmar$rownamemar),]
  pop<-sort$rownamemar
  marker<-sort$groups
  sd<-re$means[,2]
  plotdata<-data.frame(pop,marker,mean,sd)
  #开始作图
  pi <- ggplot(data_i,aes(x=Hap,y=trait))+
    geom_violin(aes(fill = Hap),position = position_dodge(width = 1), scale = 'width',
                alpha=0.8,width=0.8,show.legend = F)+
    theme(panel.grid = element_blank(), 
          panel.background = element_rect(color = "black",fill = "transparent"),
          legend.title=element_blank(),
          axis.title.x=element_blank())+
    geom_boxplot(position = position_dodge(width = 1), 
                 outlier.colour="black", 
                 outlier.shape=16,outlier.size = 0.5, width = 0.2, 
                 show.legend = FALSE,fill="lightgray")+
    geom_text(data=plotdata,aes(x=factor(pop),y=mean+sd*5,label=marker),hjust=0)+
    theme_classic()+
    ylab(real)+xlab("")
  plot_list[[i]] = pi
}
nrow_num = ceiling((length(colnames)-1)/4) 
height_num = nrow_num*2
#同时保存png与pdf
pdf(file=paste0("./05.hap_fig/",genename,".pdf"),width = 10, height = height_num)
#a<-dev.cur()   #记录pdf
#png(file=paste0("./05.hap_fig/",genename,".png"),width = 10, height = height_num)
#dev.control("enable")
new_list=plot_list[3:length(colnames)] 
grid.arrange(grobs = new_list, ncol = 4,nrow=nrow_num)
#dev.off()
dev.off()