local({r <- getOption("repos")  
r["CRAN"] <- "http://mirrors.tuna.tsinghua.edu.cn/CRAN/"   
options(repos=r)}) 
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
#install.packages("GGEBiplots")
#install.packages("GGEBiplotGUI")
library(GGEBiplotGUI)
library("GGEBiplots")
library("reshape2")
library(vegan)
library(corrplot)
library("psych")
library(Cairo)
setwd("D:/potato/")
data1<-read.table("张北data.txt",header = F,row.names = 1,comment.char = "",sep = "\t",encoding = "UTF-8")
data.t<-as.data.frame(t(data1))
data.l<-melt(data.t,id.vars = c('年度', '性状'),variable.name = "基因型")
mydata1<-dcast(data.l, `年度` + `基因型` ~ `性状`)
#`单株产量` +   `淀粉含量（%）` + `干物质含量（%）` +  `平均出苗率（%）`+ `生育日数(天)`+ `株高`
data2<-read.table("康保data.txt",header = F,row.names = 1,comment.char = "",sep = "\t",encoding = "UTF-8")
data.t<-as.data.frame(t(data2))
data.l<-melt(data.t,id.vars = c('年度', '性状'),variable.name = "基因型")
mydata2<-dcast(data.l, `年度` + `基因型` ~ `性状`)
#`单株产量` +   `淀粉含量（%）` + `干物质含量（%）` +  `平均出苗率（%）`+ `生育日数(天)` +   `株高`
#inter<-intersect(colnames(mydata1),colnames(mydata2))
#去掉一些数据
inter<-c("年度"   ,           "基因型"    ,        "单株产量"    ,     
 "淀粉含量(%)"    ,   "干物质含量(%)"  ,   "平均出苗率(%)"  ,  
"生育日数(天)"   ,  "株高" )
mydata1.c<-mydata1[,inter]
mydata2.c<-mydata2[,inter]
mydata1.c$`环境`="张北"
mydata2.c$`环境`="康保"
mydata=rbind(mydata1.c,mydata2.c)
########################计算平均值
#traits<-c("单株产量" , "淀粉含量(%)" , "干物质含量(%)" , "平均出苗率(%)", "生育日数(天)" , "小区产量" ,"折合单位产量kg/hm", "株高")

traits=c(   "单株产量"  ,        "淀粉含量(%)"  ,    
            "干物质含量(%)"  ,   "平均出苗率(%)" ,   
            "生育日数(天)"    ,       
            "株高"      )
for (i in traits){
  mydata[,i]=as.double(mydata[,i])
}
SD=round(apply(mydata[,traits],2,sd),2)
MEAN=round(apply(mydata[,traits],2,mean),2)
RANGE=round(apply(mydata[,traits],2,range),2)
CV=round(SD/MEAN*100,2)
H=round(diversity(mydata[,traits], index = "shannon", MARGIN = 2),2)
trait.mean=cbind(`平均值`=MEAN,`标准差`=SD,`变异系数（%）`=CV,`变幅`=paste(range.min=RANGE[1,],range.max=RANGE[2,],sep="-"),`多样性指数 H’`=H)
write.csv(trait.mean,file="01.马铃薯表型性状统计.csv")

#################################相关性分析
occor = corr.test(mydata[,traits],use="pairwise",method="spearman",adjust="BH",alpha=.05)
occor.r = occor$r  # 取相关性矩阵R值
occor.p = occor$p  # 取相关性矩阵p值
#设置色板
palette_2 <- colorRampPalette(c("yellow","red"))(100)
#绘图
#png(file="traits_cor.png",w=10*300,h=10*300,res = 300)
pdf(file="traits_cor.pdf",w=10,h=10,family="GB1")
corrplot(occor.r, method = "square", type = "lower", order="AOE",
         col = palette_2, tl.pos = "tp", tl.col = "blue", cl.pos = "r",
         title = "Traits corr", mar = c(2,4,4,1))
corrplot(occor.r, method = "number", type = "upper",order="AOE",
         col = palette_2,add = TRUE, diag = FALSE, tl.pos = "n", cl.pos = "n", 
         number.cex=1,mar = c(2,4,4,1))
dev.off()
write.csv(occor.p,file="02.相关性P值.csv")
write.csv(occor.r,file="03.相关性R值.csv")
########################分布范围柱状图
#########################求平均
genotype.mean=aggregate(mydata[,traits],by=list(`基因型`=mydata$`基因型`),mean)



for(t in traits){
  
  #t="平均出苗率(%)"
  d<-  genotype.mean[,t]
  ff=gsub("\\(\\S+\\)", "", t)
  ff=gsub("kg/hm", "", ff)
  
  
  pdf(paste0(ff,"-表型性状值分布.pdf"),w=5,h=5,family="GB1")
  par(mar=c(6,4.1,2,4))
  xhist <- hist(d,breaks=seq(from=floor(min(d)),to=ceiling(max(d)),length.out=10),
                plot=F,xlab =t ,ylim=c(0,),col="white",xaxt="n", yaxs="i",fill=NULL,ylab = "样本数",freq=T,main="")
  
  xhist <- hist(d,breaks=seq(from=floor(min(d)),to=ceiling(max(d)),length.out=10),
                plot=T,xlab ="" ,ylim=c(0,max(xhist$counts)*1.1),col="white",xaxt="n", yaxs="i",fill=NULL,ylab = "样本数",freq=T,main="")
  
  mtext(t,side=1,line=4.5)
  axis(side=1,xhist$mids,labels=F)
  
  b=(xhist$mids[2]-xhist$mids[1])/2
  s=round(xhist$mids-b,2)
  s[1]=0
  e=round(xhist$mids+b,2)
  lab=paste(s,e,sep="-")
  
  text(x=xhist$mids, y=-3, labels=lab,cex=1,  xpd=TRUE, srt=45,adj=1)
  
  dd=density(d)
  
  
  
  par(new=T)
  plot(dd, ylim=c(0,max(dd$y)*1.1), col ='red',lwd=2,yaxs="i",xlab = "",ylab="",xaxt="n",yaxt="n",main = "") 
  axis(4, las=1,at=seq(from = 0, to = max(xhist$density), length.out = 6), labels=round(seq(0, max(xhist$density),length.out =6)/sum(xhist$density),2), col = 'red', col.axis = 'red')
  mtext("频率",side=4,line=3,col="red")
  box(bty="l")
  dev.off()
}
##############GGEbioplots
ggdata=mydata[,c("环境","基因型","单株产量")]
ggdata1=aggregate(`单株产量` ~ 基因型 + 环境, data = mydata, mean)
gg.f=dcast(ggdata1,基因型  ~环境)
row.names(gg.f)=gg.f[,1]
gg.f=gg.f[,-1]*10
#GGEBiplot(Data = gg.f)
library(ggplot2)
GGE1<-GGEModel(gg.f,centering = "tester",
               SVP = "symmetrical",
               scaling = "none")
g1=GGEPlot(GGE1,type=6,sizeGen=4,sizeEnv = 2,largeSize=3)
g1=g1+ labs(title = "")
pdf(paste0("适应性分析bioplot.pdf"),w=5,h=5,family="GB1")
print(g1)
dev.off()

g2=GGEPlot(GGE1,type=9,sizeGen=2,sizeEnv = 3,largeSize=3)
g2=g2+ labs(title = "")
pdf(paste0("丰产性-稳定性分析bioplot.pdf"),w=5,h=5,family="GB1")
print(g2)
dev.off()

#################################表型主成分分析
p=genotype.mean
row.names(p)=p[,1]
p=p[,-1]
pca = prcomp(t(p), scale=TRUE)
#pca = prcomp(t(decostand(p, "hellinger")), scale=TRUE)
ss=summary(pca)
pc=pca$x
xx=ss$importance
row.names(xx)<-c("特征值E","贡献率CR","累计贡献率CCR")
pca.res=rbind(pc,xx)
write.csv(pca.res,file="04.PCA分析结果.csv")

#表型数据聚类分析
dist_mat <- dist(genotype.mean[,2:5], method = "euclidean")
clustering <- hclust(dist_mat, method = "complete")
plot(clustering,labels=genotype.mean$基因型,xlab="sample",ylab="Height" ,main="")
#根据进化树分组
group.data=genotype.mean
group.data$group=as.factor(cutree(clustering,h=15))

names=c("基因型"   ,     "单株产量" ,     "淀粉含量"  , "干物质含量",
        "平均出苗率" ,"生育日数"  ,"株高"    ,      "group" )
colnames(group.data)=names
library("ggpubr")
phenotype=names[2:(length(names)-1)]
test.res=data.frame()
for(i in phenotype){
  f=as.formula(paste0("`",i,"`~group"))
  print(f)
  res=compare_means(f, data = group.data,method="t.test")
  print(res)
  test.res=rbind(test.res,res)
}
write.csv(test.res,file="聚类树分组差异统计检验.csv")
write.csv(group.data,file="聚类树分组信息表.csv")
g1.res=aggregate(group.data[,2:7],by=list(group.data$group),FUN=mean)
g2.res=aggregate(group.data[,2:7],by=list(group.data$group),FUN=sd)
write.csv(g1.res,file="聚类树分组统计平均值.csv")
write.csv(g2.res,file="聚类树分组统计标准差.csv")

#表型数据在不同环境下变异
for (i in traits){
  mydata1[,i]=as.double(mydata1[,i])
}
SD=round(apply(mydata1[,traits],2,sd),2)
MEAN=round(apply(mydata1[,traits],2,mean),2)
RANGE=round(apply(mydata1[,traits],2,range),2)
CV=round(SD/MEAN*100,2)
H=round(diversity(mydata1[,traits], index = "shannon", MARGIN = 2),2)
trait.mean=cbind(`平均值`=MEAN,`标准差`=SD,`变异系数（%）`=CV,`变幅`=paste(range.min=RANGE[1,],range.max=RANGE[2,],sep="-"),`多样性指数 H’`=H)
write.csv(trait.mean,file="01.马铃薯表型性状统计(张北).csv")

#表型数据在不同环境下变异
for (i in traits){
  mydata2[,i]=as.double(mydata2[,i])
}
SD=round(apply(mydata2[,traits],2,sd),2)
MEAN=round(apply(mydata2[,traits],2,mean),2)
RANGE=round(apply(mydata2[,traits],2,range),2)
CV=round(SD/MEAN*100,2)
H=round(diversity(mydata2[,traits], index = "shannon", MARGIN = 2),2)
trait.mean=cbind(`平均值`=MEAN,`标准差`=SD,`变异系数（%）`=CV,`变幅`=paste(range.min=RANGE[1,],range.max=RANGE[2,],sep="-"),`多样性指数 H’`=H)
write.csv(trait.mean,file="01.马铃薯表型性状统计(康保).csv")

######关联分析
mydata$年度=as.factor(mydata$年度)
mydata$环境=as.factor(mydata$环境)
F.res=data.frame()
for (i in traits){
  print(i)
  fl=c(as.formula(paste0("`",i,"`~基因型")),
       as.formula(paste0("`",i,"`~环境")),
       as.formula(paste0("`",i,"`~年度")),
       as.formula(paste0("`",i,"`~基因型*环境")),
       as.formula(paste0("`",i,"`~基因型*年度"))
       #as.formula(paste0("`",i,"`~基因型*年度*环境"))
  )
  
  for(f in fl){ 
    blp=lm(f,data=mydata)
    aa=summary(blp)
    #
    #print(aa)
    print(aa$df)
    print(f)
    fstatistic = aa$fstatistic
    
    p_value = pf(as.numeric(fstatistic[1]), as.numeric(fstatistic[2]), as.numeric(fstatistic[3]), lower.tail = FALSE)
    print(p_value)
    aa=data.frame(`公式`=as.character(f)[3],df=aa$df[1],F=fstatistic[1],P=p_value,`性状`=i)
    F.res=rbind(F.res,aa)
    
    
    # question1 <- readline("Would you like to proceed untill the loop ends? (Y/N)") 
    # if(regexpr(question1, 'y', ignore.case = TRUE) == 1){ 
    #   continue = TRUE 
    #   next 
    # } else{
    #   break
    # }
  }
  
}
write.csv(F.res,file="05.表型与基因型环境年度互作效应.csv")

local({r <- getOption("repos")  
r["CRAN"] <- "http://mirrors.tuna.tsinghua.edu.cn/CRAN/"   
options(repos=r)}) 
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
#install.packages("GGEBiplots")
#install.packages("GGEBiplotGUI")
library(GGEBiplotGUI)
library("GGEBiplots")
library("reshape2")
library(vegan)
library(corrplot)
library("psych")
library(Cairo)
setwd("D:/potato/第二次分析")
mydata<-read.table("data.txt",header = T,row.names = 1,comment.char = "",sep = "\t",encoding = "UTF-8")


########################计算平均值

traits=colnames(mydata)
for (i in traits){
  mydata[,i]=as.double(mydata[,i])
}
SD=round(apply(mydata[,traits],2,sd),2)
MEAN=round(apply(mydata[,traits],2,mean),2)
RANGE=round(apply(mydata[,traits],2,range),2)
CV=round(SD/MEAN*100,2)
H=round(diversity(mydata[,traits], index = "shannon", MARGIN = 2),2)
J=round(H/log(specnumber(t(mydata[,traits]))),2)

trait.mean=cbind(`平均值`=MEAN,`标准差`=SD,`变异系数（%）`=CV,`变幅`=paste(range.min=RANGE[1,],range.max=RANGE[2,],sep="~"),`H’`=H,J=J)
write.csv(trait.mean,file="01.马铃薯表型性状统计.csv")


########################分布范围柱状图
#########################求平均


for(t in traits){
  
  ##t="顶小叶宽"
  #t="出苗率"
  d<-  mydata[,t]
  pdf(paste0(t,"-表型性状值分布.pdf"),w=5,h=5,family="GB1")
  par(mar=c(6,4.1,2,4))
  
  MIN=min(d)
  MAX=max(d)
  step=(MAX-MIN)/(length(d)^(1/2)+1)
  
  xhist <- hist(d,breaks=seq(from=floor(min(d)),to=ceiling(max(d))+step,by=step),
                plot=F,xlab =t ,ylim=c(0,),col="white",xaxt="n", yaxs="i",fill=NULL,ylab = "样本数",freq=T,main="")
  
  xhist <- hist(d,breaks=seq(from=floor(min(d)),to=ceiling(max(d))+step,by=step),
                plot=T,xlab ="" ,ylim=c(0,max(xhist$counts)*1.1),col="green",xaxt="n", yaxs="i",fill=NULL,ylab = "样本数",freq=T,main="")
  
  mtext(t,side=1,line=4.5)
  axis(side=1,xhist$mids,labels=F)
  
  
  s=round(xhist$mids-step,2)
  if (s[1]<0) {s[1]=0}
  e=round(xhist$mids+step,2)
  lab=paste(s,e,sep="-")
  
  text(x=xhist$mids, y=-3, labels=lab,cex=1,  xpd=TRUE, srt=45,adj=1)
  axis(4, las=1,at=seq(from = 0, to = max(xhist$counts), length.out = 4), labels=paste0(round(seq(0, max(xhist$counts),length.out =4)/sum(xhist$counts),4)*100,"%"), col = 'red', col.axis = 'red')
  
  dd=density(d)
  par(new=T)
  plot(dd, ylim=c(0,max(dd$y)*1.1), col ='red',lwd=2,yaxs="i",xlab = "",ylab="",xaxt="n",yaxt="n",main = "") 
  #mtext("频率",side=4,line=3,col="red")
  box(bty="l")
  l.at="topright"
  if(t=="出苗率" || t=="叶缘形状"){
    l.at="topleft"
  }
  legend(l.at,legend = c("直方图"), 
         fill = c('green'),  
         bty='n')
  legend(l.at,legend = c("正态图"), 
         col = c("red"), lty=1, 
         bty='n',inset = c(0,0.08), lwd=3)
  dev.off()
}
################相关性分析###################
#################################相关性分析
occor = corr.test(mydata,use="pairwise",method="spearman",adjust="BH",alpha=.05)
occor.r = occor$r  # 取相关性矩阵R值
occor.p = occor$p  # 取相关性矩阵p值
#设置色板
palette_2 <- colorRampPalette(c("yellow","red"))(100)
#绘图
#png(file="traits_cor.png",w=10*300,h=10*300,res = 300)
pdf(file="traits_cor.pdf",w=10,h=10,family="GB1")
corrplot(occor.r, method = "square", type = "lower", order="AOE",
         col = palette_2, tl.pos = "tp", tl.col = "blue", cl.pos = "r",
         title = "Traits corr", mar = c(2,8,4,1))
corrplot(occor.r, method = "number", type = "upper",order="AOE",
         col = palette_2,add = TRUE, diag = FALSE, tl.pos = "n", cl.pos = "n", 
         number.cex=1,mar = c(2,4,4,1))
dev.off()
write.csv(occor.p,file="02.相关性P值.csv")
write.csv(occor.r,file="03.相关性R值.csv")

###############因子分析###############
mydata.l=melt(mydata)
bt=bartlett.test(value~variable,data=mydata.l)
kmo=KMO(mydata)
sink("04.因子分析.txt")
bt
kmo
sink()


#################################表型主成分分析
#pca = prcomp(mydata, scale=TRUE)
pca = princomp(mydata, scores=TRUE, cor=TRUE)
#pca = prcomp(t(mydata), scale=TRUE)
fa1 = factanal(mydata, factor=2, rotation="varimax", scores="regression")
fa1
#write.csv(abs(fa1$loadings), "loadings.csv")
pdf(file="04.PCA分析特征根碎石图.pdf",w=8,h=4,family="GB1")
screeplot(pca, type="line", main="Scree Plot") 
dev.off()

sink(file="04.PCA分析summary.txt")
summary(pca)
sink()
pc=loadings(pca)[1:21,]
pc=as.data.frame(pc)
write.csv(pc,file="04.PCA分析结果.csv")
library(RColorBrewer)
mycolor=c(brewer.pal(9, "Set1"),brewer.pal(8, "Set2"),brewer.pal(9, "Paired"))
pdf(file="04.PCA分析表型分布.pdf",w=7,h=6,family="GB1")
par(mar=c(5,4,4,10))
plot(x=pc[,1],y=pc[,2],cex=1.5,col=mycolor,pch=16,xlab="factor 1",ylab = "factor 2")
abline(h=0,lty=1)
abline(v=0,lty=1)
legend("topright",legend =rownames(pc) ,pch=16,col=mycolor,xpd =T,inset = c(-0.3,0),ncol=1,bty='n',cex=1)
dev.off()

#表型数据聚类分析
dist_mat <- dist(mydata, method = "euclidean")
clustering <- hclust(dist_mat, method = "complete")
pdf(file="05.聚类图.pdf",w=15,h=4,family="GB1")
plot(clustering,labels=rownames(mydata),cex=0.55,xlab="sample",ylab="Height" ,main="")
abline(h=40,col="red")
dev.off()
#根据聚类树分组
group.data=mydata
group.data$group=as.factor(cutree(clustering,h=40))


###################分组PCA图################
#install.packages("ggfortify")
library(ggfortify)
pdf(file="04.PCA分析表型样本双序图.pdf",w=5,h=4,family="GB1")
#par(mar=c(5,4,4,10))
 #pcs=pca$scores
 #plot(x=pcs[,1],y=pcs[,2],cex=1.5,col=mycolor[group.data$group],pch=1,xlab="factor 1",ylab = "factor 2")
 #par(new=T)
 #plot(x=pc[,1],y=pc[,2],cex=1.5,col="blue",pch=16,xlab="factor 1",ylab = "factor 2")
#legend("topright",legend =rownames(pc) ,pch=16,col=mycolor,xpd =T,inset = c(-0.3,0),ncol=1,bty='n',cex=1)

#legend("topright",legend =1:6 ,pch=16,col=mycolor[1:6],xpd =T,inset = c(-0.3,0),ncol=1,bty='n',cex=1)
#biplot(pca)
g=autoplot(pca, data = group.data, colour = 'group',
         loadings = TRUE, loadings.colour = 'blue',
         loadings.label = TRUE, loadings.label.size = 3,size=2)#,frame = TRUE, frame.type = 'norm')
g=g+scale_color_manual(name="group",values =mycolor)+
  scale_fill_manual(name="group",values =mycolor)+
  theme_bw()+ theme(  
    panel.grid=element_blank(), 
    axis.text.x=element_text(colour="black"),
    axis.text.y=element_text(colour="black"),
    panel.border=element_rect(colour = "black"),
    legend.key = element_blank(),
    plot.title = element_text(hjust = 0.5))
g
dev.off()




library("ggpubr")
phenotype=traits
test.res=data.frame()
for(i in phenotype){
  f=as.formula(paste0("`",i,"`~group"))
  print(f)
  fit=aov(f,data=group.data)
  ss=summary(fit)
  res=as.data.frame(ss[[1]])
  res$`表型`=i
  rownames(res)=c("组间","组内")
  print(res)
  test.res=rbind(test.res,res)
}
colnames(test.res)=c("Df","平方和SS","均方MS","F","显著性","性状trait")
write.csv(test.res,file="05.聚类树分组差异统计检验.csv")
write.csv(group.data,file="05.聚类树分组信息表.csv")
g1.res=aggregate(group.data[,traits],by=list(group.data$group),FUN=mean)
Average=apply(group.data[,traits],2,mean)
Average=c("Group.1"="总体平均值 Average",Average)
g.mean=rbind(g1.res,t(as.data.frame(Average)))
write.csv(t(g.mean),file="06.聚类树分组统计平均值.csv")