#CROPS545_Lecture23_CV #Zhiwu Zhang #April 3, 2016 rm(list=ls()) #Import GAPIT #source("http://www.bioconductor.org/biocLite.R") #biocLite("multtest") #install.packages("gplots") #install.packages("scatterplot3d")#The downloaded link at: http://cran.r-project.org/package=scatterplot3d library('MASS') # required for ginv library(multtest) library(gplots) library(compiler) #required for cmpfun library("scatterplot3d") source("http://www.zzlab.net/GAPIT/emma.txt") #source("http://www.zzlab.net/GAPIT/gapit_functions.txt") #source("~/Dropbox/GAPIT/Functions/emma.txt") source("~/Dropbox/GAPIT/Functions/gapit_functions.txt") #Import demo data myGD=read.table(file="http://zzlab.net/GAPIT/data/mdp_numeric.txt",head=T) myGM=read.table(file="http://zzlab.net/GAPIT/data/mdp_SNP_information.txt",head=T) myCV=read.table(file="http://zzlab.net/GAPIT/data/mdp_env.txt",head=T) #myGD=read.table(file="~/Dropbox/Current/ZZLab/WSUCourse/CROPS545/Demo/mdp_numeric.txt",head=T) #myGM=read.table(file="~/Dropbox/Current/ZZLab/WSUCourse/CROPS545/Demo/mdp_SNP_information.txt",head=T) #Simultate 10 QTN on the first half chromosomes X=myGD[,-1] index1to5=myGM[,2]<6 X1to5 = X[,index1to5] taxa=myGD[,1] set.seed(99164) GD.candidate=cbind(taxa,X1to5) source("~/Dropbox/GAPIT/Functions/GAPIT.Phenotype.Simulation.R") mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=10, effectunit =.95,QTNDist="normal",CV=myCV,cveff=c(.51,.51)) setwd("~/Desktop/temp") source("~/Dropbox/GAPIT/Functions/gapit_functions.txt") source("~/Dropbox/GAPIT/Functions/GAPIT.Main.R") #Prediction by PCA and ENV #---------------------------------------------------------------------------------------- myGAPIT <- GAPIT( Y=mySim$Y, GD=myGD, GM=myGM, PCA.total=3, CV=myCV, #CV.Inheritance=6, group.from=7, group.to=7, group.by=10, QTN.position=mySim$QTN.position, #SNP.test=FALSE, memo="GLM5", ) str(myGAPIT) ry2=cor(myGAPIT$Pred[,8],mySim$Y[,2])^2 ru2=cor(myGAPIT$Pred[,8],mySim$u)^2 par(mfrow=c(2,1), mar = c(3,4,1,1)) plot(myGAPIT$Pred[,8],mySim$Y[,2]) mtext(paste("R square=",ry2,sep=""), side = 3) plot(myGAPIT$Pred[,8],mySim$u) mtext(paste("R square=",ru2,sep=""), side = 3) #Adding the top associated SNPs ntop=200 index=order(myGAPIT$P) top=index[1:ntop] myQTN=cbind(myGAPIT$PCA[,1:4], myCV[,2:3],myGD[,c(top+1)]) source("~/Dropbox/GAPIT/Functions/GAPIT.EMMAxP3D.R") source("~/Dropbox/GAPIT/Functions/GAPIT.Main.R") myGAPIT2 <- GAPIT( Y=mySim$Y, GD=myGD, GM=myGM, #PCA.total=3, CV=myQTN, group.from=1, group.to=1, group.by=10, QTN.position=mySim$QTN.position, SNP.test=FALSE, memo="GLM+QTN", ) ry2=cor(myGAPIT2$Pred[,8],mySim$Y[,2])^2 ru2=cor(myGAPIT2$Pred[,8],mySim$u)^2 par(mfrow=c(2,1), mar = c(3,4,1,1)) plot(myGAPIT2$Pred[,8],mySim$Y[,2]) mtext(paste("R square=",ry2,sep=""), side = 3) plot(myGAPIT2$Pred[,8],mySim$u) mtext(paste("R square=",ru2,sep=""), side = 3) #Cross validation set.seed(99164) n=nrow(mySim$Y) training=sample(n,round(n/2),replace=F) myGAPIT3 <- GAPIT( Y=mySim$Y[training,], GD=myGD, GM=myGM, CV=myCV, PCA.total=3, group.from=1, group.to=1, group.by=10, QTN.position=mySim$QTN.position, #SNP.test=FALSE, memo="Training", ) ry2=cor(myGAPIT3$Pred[,8],mySim$Y[,2])^2 ru2=cor(myGAPIT3$Pred[,8],mySim$u)^2 par(mfrow=c(2,1), mar = c(3,4,1,1)) plot(myGAPIT3$Pred[,8],mySim$Y[,2]) mtext(paste("R square=",ry2,sep=""), side = 3) plot(myGAPIT3$Pred[,8],mySim$u) mtext(paste("R square=",ru2,sep=""), side = 3) #Adding the top associated SNPs ntop=10 index=order(myGAPIT3$P) top=index[1:ntop] myQTN=cbind(myGAPIT3$PCA[,1:4], myCV[,2:3],myGD[,c(top+1)]) source("~/Dropbox/GAPIT/Functions/GAPIT.Main.R") source("~/Dropbox/GAPIT/Functions/GAPIT.EMMAxP3D.R") myGAPIT4 <- GAPIT( Y=mySim$Y[training,], GD=myGD, GM=myGM, #PCA.total=3, CV=myQTN, group.from=1, group.to=1, group.by=1, QTN.position=mySim$QTN.position, SNP.test=FALSE, memo="GLM+QTN", ) ry2=cor(myGAPIT4$Pred[training,8],mySim$Y[training,2])^2 ru2=cor(myGAPIT4$Pred[training,8],mySim$u[training])^2 par(mfrow=c(2,1), mar = c(3,4,1,1)) plot(myGAPIT4$Pred[training,8],mySim$Y[training,2]) mtext(paste("R square=",ry2,sep=""), side = 3) plot(myGAPIT4$Pred[training,8],mySim$u[training]) mtext(paste("R square=",ru2,sep=""), side = 3) quartz() #cor(myQTN[,-1]) heatmap(cor(myQTN[,-1])) #Adding the top associated SNPs ntop=5 index=order(myGAPIT3$P) top=index[1:ntop] myQTN=cbind(myGAPIT3$PCA[,1:4], myCV[,2:3],myGD[,c(top+1)]) #calculate prediction effect=rbind(myGAPIT3$effect.cv, as.matrix(myGAPIT3$effect.snp[top])) X=as.matrix(cbind(1, myCV[,2:3],myGAPIT3$PCA[,-1],myGD[,c(top+1)])) Pred=X%*%effect ry2=cor(Pred[training],mySim$Y[training,2])^2 ru2=cor(Pred[training],mySim$u[training])^2 par(mfrow=c(2,1), mar = c(3,4,1,1)) plot(Pred[training],mySim$Y[training,2]) mtext(paste("R square=",ry2,sep=""), side = 3) plot(Pred[training],mySim$u[training]) mtext(paste("R square=",ru2,sep=""), side = 3)