# Load necessary packages
library(parallel)
library(doParallel)
source("http://zzlab.net/GAPIT/gapit_functions.txt")

# Import demo data
myGD <- read.table(file = "http://zzlab.net/GAPIT/data/mdp_numeric.txt", head = TRUE)
myGM <- read.table(file = "http://zzlab.net/GAPIT/data/mdp_SNP_information.txt", head = TRUE)

# Simulate 10 QTN on the first half chromosomes
index1to5 <- myGM[, 2] < 6
set.seed(99164)
mySim <- GAPIT.Phenotype.Simulation(GD = myGD[, c(TRUE, index1to5)], GM = myGM[index1to5, ],
                                   h2 = 0.7, NQTN = 40, effectunit = 0.95, QTNDist = "normal")

# GWAS with GAPIT using parallel processing
num_cores <- detectCores()  # Detect the number of available CPU cores
cl <- makeCluster(num_cores)  # Create a cluster with all CPU cores
models = c("GLM", "MLM", "CMLM", "SUPER", "MLMM", "FarmCPU", "Blink")
registerDoParallel(cl)  # Register the cluster for parallel processing

myGAPIT <- foreach(model = c("GLM", "MLM", "CMLM", "SUPER", "MLMM", "FarmCPU", "Blink"), .combine = c) %dopar% {
  GAPIT(Y = mySim$Y, GD = myGD, GM = myGM, PCA.total = 3,
        QTN.position = mySim$QTN.position, model = model)
}

stopCluster(cl)  # Stop the cluster


### Output multiple Manhattans with GWAS resutls, which were read from the located direction

   GMM=GAPIT.Multiple.Manhattan(model_store=models,Y.names=colnames(mySim$Y)[-1],GM=myGM,seqQTN = mySim$QTN.position,plot.type=c("w","h"))
   GAPIT.Circle.Manhattan.Plot(band=1,r=3,GMM$multip_mapP,plot.type=c("c","q"),signal.line=1,xz=GMM$xz,threshold=0.01)
   GMM=GAPIT.Multiple.Manhattan(model_store=models,Y.names=colnames(mySim$Y)[-1],GM=myGM,seqQTN = mySim$QTN.position,plot.type=c("s"))