rm(list=ls()) library(compiler) #required for cmpfun source("http://www.zzlab.net/GAPIT/gapit_functions.txt") 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) setwd("~/Dropbox/Current/ZZLab/WSUCourse/CROPS545/R/2016") source("G2P.R") source("GWASbyCor.R") X=myGD[,-1] index1to5=myGM[,2]<6 X1to5 = X[,index1to5] set.seed(99164) mySim=G2P(X= X1to5,h2=.75,alpha=1,NQTN=10,distribution="norm") p= GWASbyCor(X=X,y=mySim$y) index=order(p) top5=index[1:5] detected=intersect(top5,mySim$QTN.position) falsePositive=setdiff(top5, mySim$QTN.position) top5 mySim$QTN.position detected length(detected) falsePositive color.vector <- rep(c("deepskyblue","orange","forestgreen","indianred3"),10) m=nrow(myGM) plot(t(-log10(p))~seq(1:m),col=color.vector[myGM[,2]]) abline(v=mySim$QTN.position, lty = 2, lwd=2, col = "black") abline(v= falsePositive, lty = 2, lwd=2, col = "red") index.null=!index1to5 & !is.na(p) p.null=p[index.null] m.null=length(p.null) index.sort=order(p.null) p.null.sort=p.null[index.sort] head(p.null.sort) tail(p.null.sort) seq=seq(1:m.null) table=cbind(seq, p.null.sort, seq/m.null) head(table,15) tail(table,15) set.seed(99164) mySim=G2P(X= myGD[,-1],h2=.75,alpha=1,NQTN=10,distribution="norm") p= GWASbyCor(X=X,y=mySim$y) plot(t(-log10(p))~seq(1:m),col=color.vector[myGM[,2]]) abline(v=mySim$QTN.position, lty = 2, lwd=2, col = "black") bigNum=1e9 resolution=100000 bin=round((myGM[,2]*bigNum+myGM[,3])/resolution) result=cbind(myGM,t(p),bin) head(result) QTN.bin=result[mySim$QTN.position,] QTN.bin index.qtn.p=order(QTN.bin[,4]) QTN.bin[index.qtn.p,] myStat=GAPIT.FDR.TypeI( WS=c(1e0,1e3,1e4,1e5), GM=myGM, seqQTN=mySim$QTN.position, GWAS=result) str(myStat) 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=100 set.seed(99164) statRep=replicate(nrep, { mySim=G2P(X=myGD[,-1],h2=.75,alpha=1,NQTN=10,distribution="norm") p=p= GWASbyCor(X=myGD[,-1],y=mySim$y) seqQTN=mySim$QTN.position myGWAS=cbind(myGM,t(p),NA) myStat=GAPIT.FDR.TypeI(WS=c(1e0,1e3,1e4,1e5), GM=myGM,seqQTN=mySim$QTN.position,GWAS=myGWAS,maxOut=100,MaxBP=1e10) }) str(statRep) power=statRep[[2]] s.fdr=seq(3,length(statRep),7) fdr=statRep[s.fdr] fdr.mean=Reduce ("+", fdr) / length(fdr) s.tp1=seq(4,length(statRep),7) tp1=statRep[s.tp1] tp1.mean=Reduce ("+", tp1) / length(tp1) s.auc.fdr=seq(6,length(statRep),7) auc.fdr=statRep[s.auc.fdr] auc.fdr.mean=Reduce ("+", auc.fdr) / length(auc.fdr) #quartz() 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]) } #barplot on AUC barplot(auc.fdr.mean, names.arg=c("1bp", "1K", "10K","100K"), xlab="Resolution", ylab="AUC") nrep=100 set.seed(99164) #h2=25% statRep25=replicate(nrep, { mySim=G2P(X=myGD[,-1],h2=.25,alpha=1,NQTN=10,distribution="norm") p=p= GWASbyCor(X=myGD[,-1],y=mySim$y) seqQTN=mySim$QTN.position myGWAS=cbind(myGM,t(p),NA) myStat=GAPIT.FDR.TypeI(WS=c(1e0,1e3,1e4,1e5), GM=myGM,seqQTN=mySim$QTN.position,GWAS=myGWAS,maxOut=100,MaxBP=1e10)}) #h2=75% statRep75=replicate(nrep, { mySim=G2P(X=myGD[,-1],h2=.75,alpha=1,NQTN=10,distribution="norm") p=p= GWASbyCor(X=myGD[,-1],y=mySim$y) seqQTN=mySim$QTN.position myGWAS=cbind(myGM,t(p),NA) myStat=GAPIT.FDR.TypeI(WS=c(1e0,1e3,1e4,1e5), GM=myGM,seqQTN=mySim$QTN.position,GWAS=myGWAS,maxOut=100,MaxBP=1e10)}) power25=statRep25[[2]] s.t1=seq(4,length(statRep25),7) t1=statRep25[s.t1] t1.mean.25=Reduce ("+", t1) / length(t1) power75=statRep75[[2]] s.t1=seq(4,length(statRep75),7) t1=statRep75[s.t1] t1.mean.75=Reduce ("+", t1) / length(t1) plot(t1.mean.25[,4],power25, type="b", col="blue",xlim=c(0,1)) lines(t1.mean.75[,4], power75, type="b", col= "red")