#Demonstration code to use FarmCPU (Fixed and random model Circulating Probability Unification)
#Designed by Xiaolei Liu and Zhiwu Zhang
#Writen by Xiaolei Liu and Zhiwu Zhang
#Last update: January 1, 2016

#Installation of required R package (once for a computer, details see FarmCPU user manual.)
#  install.packages("bigmemory")
#  install.packages("biganalytics")

rm(list=ls())

#Import library (each time to start R)
library(bigmemory)
library(biganalytics)
require(compiler) #for cmpfun

#Import GAPIT
source("http://www.zzlab.net/GAPIT/gapit_functions.txt")#web source code
#source("gapit_functions.txt")#users can download the source code and save in a .txt file
#Import FarmCPU
source("http://www.zzlab.net/FarmCPU/FarmCPU_functions.txt")#web source code
#source("FarmCPU_functions.txt")#users can download the source code and save in a .txt file
#import data
#FarmCPU uses four data sets: genotype, phenotype, map and covariates
#Phenotype: details see FarmCPU user manual 3.1
#Genotype: details see FarmCPU user manual 3.2
#Map: details see FarmCPU user manual 3.3
#Covariates: details see FarmCPU user manual 3.4
myY <- read.table("mdp_traits_validation.txt", head = TRUE)
myGD <- read.big.matrix("mdp_numeric.txt",type = "char",sep = "\t",header = T)#Attention the type, "char" or "double"
myGM <-  read.table("mdp_SNP_information.txt" , head = TRUE)
#myCV <- read.table("PCA.txt", head = TRUE)#You can use myGD and myGM data and perform GAPIT to generate a principal components file PCA.txt

#GWAS with FarmCPU...
myFarmCPU=FarmCPU(
Y=myY[,c(1,8)],#Phenotype, FarmCPU only accept single phenotype
GD=myGD,#Genotype
GM=myGM,#Map information
#CV=myCV[,2:4], #FarmCPU does not need the taxa information
#method.bin="optimum", #options are "static" and "optimum", default is static and this gives the fastest speed. If you want to use random model to optimize possible QTNs selection, use method.bin="optimum"
#bin.size=c(5e5,5e6,5e7),#bin.size is used to set for the bin size
#bin.selection=seq(10,100,10),#bin.selection is used to set possible QTN number
#Prior=NULL,#Prior is used to add "QTNs" as prior knowledge, details see FarmCPU user manual
#maxLoop=5,#maxLoop is used to set the maximum iterations you want
#iteration.output=TRUE,#iteration.output=TRUE means to output results of every iteration
#threshold.output=1,#FarmCPU can only output SNPs with a more significant p value than threshold.output, this is used to save time
#MAF.calculate=TRUE,#FarmCPU will calcualte the MAF value of each SNP and filter the SNPs that have a lower MAF value than "maf.threshold"
#maf.threshold=0,
#QTN.position=c(1257,1973,2676,1754,1533,1644,2232,2663,1045,2286),#this is only for mark the QTN in manhattan plot for the 7th phenotype: myY[,c(1,8)]
)