#Study of FDR and Power
#Written By Qishan Wang, Meng Li and Zhiwu Zhang
#Last update: October 26, 2014

#Install packages (Do this section only for new installation of R)
#-------------------------------------------------------------------------------
#source("http://www.bioconductor.org/biocLite.R") 
#biocLite("multtest")
#install.packages("gplots")

#Step 0: Import library and GAPIT functions run this section each time to start R)
#######################################################################################
library('MASS') # required for ginv
library(multtest)
library(gplots)
library(compiler) #required for cmpfun
source("http://www.zzlab.net/GAPIT/emma.txt")
source("http://www.zzlab.net/GAPIT/gapit_functions.txt")
source("/Users/Zhiwu/Dropbox/Current/revolutionr/gapit/gapit_functions.txt")
#############################################################################################

#Step 1: Prepare data 
#download tutorial data and save them in myGAPIT directory under C drive and  run tutorials
setwd("C:\\myGAPIT")
myG <- read.delim("mdp_genotype_test.hmp.txt" , head = FALSE)

#Set working directory
setwd("/Users/zhiwu/temp")

#Convert to numerical fomrat
myGAPIT0=GAPIT(G=myG,PCA.total=20,file.output=FALSE)
myGD=myGAPIT0$GD
myGM=myGAPIT0$GI
myKI=myGAPIT0$kinship
myCV=myGAPIT0$PCA[,1:4]

#Step 1: Define simulation parameters
#Parameters
#-------------------------------------------------------------------------------
model="QOnly" #This is for documentation
theRep=10 ###########change the number of rep here
h2=0.75 ######change heritability here
NQTN=10 ###########change the number of QTN here
effectunit=1
WS=c(1e0,1e3,1e4,1e5,1e6,1e7)    #window size series
alpha=c(.01,.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)
theSeed=12345     #seed for sampling
MaxBP=1e10
QTNDist="geometry" #options are "normal" and "geometry"
maxOut=100 #Number of rows in power FDR report
acceleration=0
maxLoop=10
method.sub="reward"#options are penalty, reward, mean and median
method.sub.final="reward"#options are penalty, reward, mean and median
alpha=c(.01,.05,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)
WS=c(1e0,1e3,1e4,1e5,1e6,1e7)
maxOut=100
memo=paste(model,method.sub,method.sub.final,acceleration,theRep,h2,NQTN,sep="_")


#Step 2: Warm up (Skip this section in final simulation)
myQTN.position=sample(nrow(myGM),10,replace=F)
myGWAS=cbind(myGM,runif(nrow(myGM)),1)
myGWAS[myQTN.position,4]=myGWAS[myQTN.position,4]/100000
myPower=GAPIT.Power(WS=WS, alpha=alpha, maxOut=maxOut,seqQTN=myQTN.position,GM=myGM,GWAS=myGWAS)
head(myPower$Power.Alpha)
head(myPower$FDR)
head(myPower$Power)
#Testing Manhattan and QQ plots
GAPIT.Manhattan(GI.MP = myGWAS[,2:4], name.of.trait = "Trait",plot.type = "Genomewise",
                DPP=50000,cutOff=0.01,band=5,seqQTN=myQTN.position)
GAPIT.QQ(P.values=myGWAS[,4], plot.type = "log_P_values", name.of.trait = "Trait",DPP=50000)


#Step 3: Excuse simulation
set.seed(theSeed)
print("Iterating on rep...")
i=1
for(i in 1:theRep){
infor=paste("$$$$$$$$$$$$$$$$$$ ",i," $$$$$$$$$$$$$$$$$$$$$$$",sep="")
print(infor)

#Simulating phenotype
print("Simulating phenotype...")
#source("/Users/zhiwu/Dropbox/Current/paper/BigData/BUS/FARM-CPU/Develope.Phenotype.Simulation.R")
myPheno=GAPIT.Phenotype.Simulation(GD=myGD,h2=h2,NQTN=NQTN,QTNDist=QTNDist,effectunit=effectunit,theSeed=theSeed)
myY=myPheno$Y
myU=myPheno$effect
myE=myPheno$residual
myQTN.position=myPheno$QTN.position
myEffect=myPheno$addeffect
colnames(myY)=c("taxa",paste("Sim",i,sep=""))

#Hiding QTNs (Do not perform this on big.matrix as it change file on disk)
myGD1=myGD
#myGD1[,(myQTN.position+1)]=0 #Remove variation

#marker density
#index=c(1,seq(2,ncol(myGD),by=2))

#Shuffle phenotype
#n=nrow(myY)
#s=sample(1:n,n)    
#myY=myY[s,]


print("GWAS by GAPIT...")
#GWAS with GAPIT (modeling here)
#-------------------------------------------------------------------------------
#Step 1: Run BUSS
myGAPIT=GAPIT(
Y=myY[,c(1,2)],
GD=myGD1,
GM=myGM,
KI=myKI,
CV=myCV,
#PCA.total=3,
#file.output=FALSE,
group.from=1,
group.to=1,
group.by=1,

memo=memo,
QTN.position=myQTN.position,
threshold.output=0.001,
iteration.output=TRUE,
) 

#Save and calculate average power and FDR
#-------------------------------------------------------------------------------
if(i==1){
  myPower.Alpha=myGAPIT$Power.Alpha
  myPower=myGAPIT$Power
  myFDR=myGAPIT$FDR
}else{
  myPower.Alpha=((i-1)*myPower.Alpha+myGAPIT$Power.Alpha)/i
  myPower=((i-1)*myPower+myGAPIT$Power)/i
  myFDR=((i-1)*myFDR +   myGAPIT$FDR)/i
}

}#end of rep loop 

#Step 4: Ouput simulation results
#Final report on means
#-------------------------------------------------------------------------------

colnames(myFDR)=  paste("FDR(",WS,")",sep="")
colnames(myPower)=paste("Power(",WS,")",sep="")
colnames(myPower.Alpha)=paste("Power(",WS,")",sep="")
write.csv(cbind(myFDR,myPower),paste("GAPIT.Power.by.FDR.",memo,".csv",sep=""))
write.csv(cbind(alpha,myPower.Alpha),paste("GAPIT.Power.by.TypeI.",memo,".csv",sep=""))

print("Simulation accoplished")