# Please do the following before you get started:
# 1. create a folder for your PLINK analysis and save it as your working directory
# 2. import data
# 3. load GAPIT functions or G2P function, we only need GAPIT for simulation


# 1. test PLINK with its test data----
setwd(MAIN)
setwd("./plink_mac")
# use the plink test file to run a GWAS
system("./plink --file toy --assoc")
# if good, you get a .log, and a .assoc

setwd(PLINK)
# PLINK has its own format requirement for files.
# We will get files ready before GWAS.
# Phenotype file is prepared by the following R code
# Genotype and map file are prepared by iPat

# Phenotype file:
pheno=GAPIT.Phenotype.Simulation(GD=myGD,GM=myGM,h2=.5,NQTN=10,QTNDist="norm")
data.Y = data.frame(FID = pheno$Y[,1], SID = pheno$Y[,1], phenotype = pheno$Y[,-1])
write.table(x = data.Y, file = "pheno.txt", row.names = F, quote = F, sep = "\t")

# GD and GM
write.table(myGD, file = "genotype.txt", row.names=F, quote = F, sep = "\t")
write.table(myGM, file = "mapinfo.txt", row.names = F, quote = F, sep = "\t")

# file.copy(file.path(MAIN,"plink_mac/plink"), to="plink")
# geno: missing rate, if less than 0.2, filter
# maf: minor allel frequency, if less than 0.05, filter
# system("./plink --file genotype_recode --allow-no-sex --assoc --geno 0 --maf 0.05 --pheno pheno.txt")
system("./plink --ped myGD_recode.ped --assoc --map myGD_recode.map --allow-no-sex --adjust --pheno pheno.txt --all-pheno -out outP --geno 0.2 --maf 0.5")

output = read.table("outP.P1.qassoc", header = T)
p=output$P
length(p)
plinkp = data.frame(myGM, P=doutput$P)
# system("./plink --ped myGD_recode.ped --assoc --map myGD_recode.map --allow-no-sex
#       --adjust --pheno pheno.txt --all-pheno -out outP --geno 0.2 --maf 0.05")