rm(list=ls()) setwd('/Users/Zhiwu/Dropbox/Current/ZZLab/WSUCourse/CROPS545/mlmm-master') source('mlmm_cof.r') 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("/Users/Zhiwu/Dropbox//GAPIT/functions/gapit_functions.txt") setwd("/Users/Zhiwu/Dropbox/Current/ZZLab/WSUCourse/CROPS512/Demo") myGD <- read.table("mdp_numeric.txt", head = TRUE) myGM <- read.table("mdp_SNP_information.txt" , head = TRUE) #for PC and K setwd("~/Desktop/temp") myGAPIT0=GAPIT(GD=myGD,GM=myGM,PCA.total=3,) myPC=as.matrix(myGAPIT0$PCA[,-1]) myK=as.matrix(myGAPIT0$kinship[,-1]) myX=as.matrix(myGD[,-1]) #Siultate 10 QTN on the first chromosomes X=myGD[,-1] index1to5=myGM[,2]<6 X1to5 = X[,index1to5] taxa=myGD[,1] set.seed(99164) GD.candidate=cbind(taxa,X1to5) mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=10,QTNDist="norm") myy=as.numeric(mySim$Y[,-1]) myMLMM<-mlmm_cof(myy,myX,myPC[,1:2],myK,nbchunks=2,maxsteps=20) myP=myMLMM$pval_step[[1]]$out[,2] myGI.MP=cbind(myGM[,-1],myP) setwd("~/Desktop/temp") GAPIT.Manhattan(GI.MP=myGI.MP,seqQTN=mySim$QTN.position) GAPIT.QQ(myP) myGWAS=cbind(myGM,myP,NA) myStat=GAPIT.FDR.TypeI(WS=c(1e0,1e3,1e4,1e5),GM=myGM,seqQTN=mySim$QTN.position,GWAS=myGWAS) par(mfrow=c(1,2),mar = c(5,2,5,2)) plot(myStat$FDR[,1],myStat$Power,type="b") plot(myStat$TypeI[,1],myStat$Power,type="b") nrep=10 set.seed(99164) statRep=replicate(nrep, { mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=10,QTNDist="norm") myy=as.numeric(mySim$Y[,-1]) myMLMM<-mlmm_cof(myy,myX,myPC[,1:2],myK,nbchunks=2,maxsteps=20) myP=myMLMM$pval_step[[1]]$out[,2] myGWAS=cbind(myGM,myP,NA) myStat=GAPIT.FDR.TypeI(WS=c(1e0,1e3,1e4,1e5),GM=myGM,seqQTN=mySim$QTN.position,GWAS=myGWAS) }) power=statRep[[2]] #FDR s.fdr=seq(3,length(statRep),7) fdr=statRep[s.fdr] fdr.mean=Reduce ("+", fdr) / length(fdr) #AUC: power vs. FDR s.auc.fdr=seq(6,length(statRep),7) auc.fdr=statRep[s.auc.fdr] auc.fdr.mean=Reduce ("+", auc.fdr) / length(auc.fdr) theColor=rainbow(4) plot(fdr.mean[,1],power , type="b", col=theColor [1],xlim=c(0,1)) for(i in 2:ncol(fdr.mean)){ lines(fdr.mean[,i], power , type="b", col= theColor [i]) }