#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") library("EMMREML") 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) #myCV=read.table(file="~/Dropbox/Current/ZZLab/WSUCourse/CROPS545/Demo/mdp_env.txt",head=T) #Simultate 10 QTN on the first half chromosomes X=myGD[,-1] index1to5=myGM[,2]<6 X1to5 = X[,index1to5] taxa=myGD[,1] GD.candidate=cbind(taxa,X1to5) set.seed(99164) mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=20, effectunit =.95,QTNDist="normal",CV=myCV,cveff=c(.51,.51)) #order phenotype by ID index=order(mySim$Y[,1]) myY=mySim$Y[index,] colnames(myY)=c("Taxa","trait") n=nrow(myY) #Get PC and merge to CV myGAPIT0=GAPIT(GD=myGD,GM=myGM,PCA.tota =3,file.out=FALSE) myPCCV=cbind(myGAPIT0$PCA, myCV[,-1]) myKI=myGAPIT0$kinship #Configration setwd("~/Desktop/temp") #compare cblup and gblup nrep=30 nfold=3 set.seed(99164) t=0 for (rep in 1:nrep){ #setup fold myCut=cut(1:n,nfold,labels=FALSE) fold=sample(myCut,n) for (f in 1:nfold){ t=t+1 testing=(fold==f) training=!testing myGAPIT <- GAPIT( Y=myY[training,], GM=myGM, GD=myGD, KI=myKI, CV=myPCCV, group.from=1000, group.to=1000, group.by=10, SNP.test=FALSE, SUPER_GS=TRUE, file.output=FALSE ) #Calculate correlation r.ref=cor(myGAPIT$Pred[training,5],mySim$u[training])^2 r.inf=cor(myGAPIT$Pred[!training,5],mySim$u[!training])^2 acc=cbind(r.ref, r.inf) #calculate average across fold and rep print(c(rep, f)) if(rep==1 & f==1){acc.avg=acc}else{acc.avg=acc/t+acc.avg*((t-1)/t)} }#fold }#rep acc.sblup=acc.avg acc.gblup=acc.avg acc.cblup=acc.avg