R语言农艺性状表型数据与环境关联互作代码分析

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")

相关阅读