#CROPS545_Lecture25_RR #Zhiwu Zhang #April 3, 2016 rm(list=ls()) #Import GAPIT #source("http://www.bioconductor.org/biocLite.R") #biocLite("multtest") #install.packages("EMMREML") #install.packages("gplots") #install.packages("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") source("~/Dropbox/GAPIT/Functions/GAPIT.Phenotype.PCA.View.R") #Import BLR #install.packages("BLR") library(BLR) #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) #Preparing data X=myGD[,-1] taxa=myGD[,1] index1to5=myGM[,2]<6 X1to5 = X[,index1to5] GD.candidate=cbind(as.data.frame(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(.2,.2),a2=.5,adim=3,category=1,r=.4) #Validation set.seed(99164) n=nrow(mySim$Y) pred=sample(n,round(n/5),replace=F) train=-pred #MLM myY=mySim$Y y <- mySim$Y[,2] setwd("~/Desktop/temp") #---------------------------------------------------------------------------------------- myGAPIT <- GAPIT( Y=myY[train,], GD=myGD, GM=myGM, PCA.total=3, CV=myCV, group.from=1000, group.to=1000, group.by=10, QTN.position=mySim$QTN.position, #SNP.test=FALSE, memo="MLM") order.raw=match(taxa,myGAPIT$Pred[,1]) pcEnv=cbind(myGAPIT$PCA[,-1],myCV[,-1]) acc.GAPIT=cor(myGAPIT$Pred[order.raw,5][pred],mySim$u[pred]) fit.GAPIT=cor(myGAPIT$Pred[order.raw,5][train],mySim$u[train]) acc.GAPIT fit.GAPIT nIter=2000 #### number of iteration burnIn=1500 #### burnin a part of iteration myBLR=BLR(y=as.matrix(myY[train,2]),XF=pcEnv[train,],XL=as.matrix(myGD[train,-1]),nIter=nIter,burnIn=burnIn) #myBLR=BLR(y=as.matrix(myY[train,2]),XF=pcEnv[train,],XR=as.matrix(myGD[train,-1]),nIter=nIter,burnIn=burnIn) #myBLR=BLR(y=as.matrix(myY[train,2]),XF=pcEnv[train,],XR=as.matrix(myGD[train,-1]),XL=as.matrix(myGD[train,-1]),nIter=nIter,burnIn=burnIn) pred.inf=as.matrix(myGD[pred,-1])%*%myBLR$bL pred.ref=as.matrix(myGD[train,-1])%*%myBLR$bL accuracy <- cor(myY[pred,2],pred.inf) modelfit <- cor(myY[train,2],pred.ref) hat_h2=var(pred.ref)/(var(pred.ref)+myBLR$varE) accuracy modelfit hat_h2 plot(myBLR$tau2) plot(-log10(1/(exp(10000*myBLR$tau2)))) plot(myBLR$SD.bL) plot(myBLR$SD.bL,myBLR$tau2) #Extract p values myP=1/(exp(10000*myBLR$tau2)) myGI.MP=cbind(myGM[,-1],myP) GAPIT.Manhattan(GI.MP=myGI.MP,seqQTN=mySim$QTN.position) GAPIT.QQ(myP)