#Comparing models
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/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)

#Siultate 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)
mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=10,QTNDist="norm")

#Kinship from QTNs
setwd("~/Desktop/temp")
myGAPIT0=GAPIT(
  Y=mySim$Y,
  GD=GD.candidate,
  GM=myGM[index1to5,],
  SNP.test = F,
  group.from=5,
  group.to=5,
  group.by=10,
  memo="KbyQTN")

myGAPIT1=GAPIT(
  Y=mySim$Y,
  GD=myGD,
  GM=myGM,
  KI=myGAPIT0$kinship,
  QTN.position=mySim$QTN.position,
  PCA.total=3,
  #SNP.test = F,
  group.from=1000000,
  group.to=1000000,
  group.by=10,
  memo="GWASbyQTN")

cor(mySim$u,myGAPIT0$Pred$BLUP)
cor(mySim$u,myGAPIT1$Pred$BLUP)
par(mfrow=c(1,2))
plot(mySim$u,myGAPIT0$Pred$BLUP)
plot(mySim$u,myGAPIT1$Pred$BLUP)

#RUN t test
setwd("~/Desktop/temp")
myGAPIT=GAPIT(
  Y=mySim$Y,
  GD=myGD,
  GM=myGM,
  QTN.position=mySim$QTN.position,
  PCA.total=0,
  group.from=1,
  group.to=1,
  group.by=10,
  memo="ttest")

#RUN SUPER
myGAPIT=GAPIT(
  Y=mySim$Y,
  GD=myGD,
  GM=myGM,
  QTN.position=mySim$QTN.position,
  PCA.total=3,
  sangwich.top="MLM", #options are GLM,MLM,CMLM, FaST and SUPER 
  sangwich.bottom="SUPER", #options are GLM,MLM,CMLM, FaST and SUPER 
  LD=0.1,
  memo="SUPER")

#Power
myStat=GAPIT.FDR.TypeI(
WS=c(1e0,1e3,1e4,1e5),
GM=myGM,
seqQTN=mySim$QTN.position,
GWAS=myGAPIT$GWAS)

#Replicate
nrep=3
set.seed(99164)
statRep=replicate(nrep,{
mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=10,QTNDist="norm")
myGAPIT=GAPIT(Y=mySim$Y,  GD=myGD,  GM=myGM,  QTN.position=mySim$QTN.position,  PCA.total=3,  
sangwich.top="MLM", #options are GLM,MLM,CMLM, FaST and SUPER   
sangwich.bottom="SUPER", #options are GLM,MLM,CMLM, FaST and SUPER   
LD=0.1,  memo="SUPER")
myStat=GAPIT.FDR.TypeI(WS=c(1e0,1e3,1e4,1e5),GM=myGM,seqQTN=mySim$QTN.position,GWAS=myGAPIT$GWAS)})
  
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])
}