In stp4freedom/GWhEAT: Performed GWAS in a fast and efficient manner

if (!file.exists("GAPIT_Tutorial_Data.zip"))
{
  download.file("http://zzlab.net/GAPIT/GAPIT_Tutorial_Data.zip", destfile = "GAPIT_Tutorial_Data.zip")
  unzip("GAPIT_Tutorial_Data.zip")
}
download.file("http://zzlab.net/GAPIT/data/CROP545_Covariates.txt", destfile = "CROPS545_Covariates.txt")
download.file("http://zzlab.net/GAPIT/data/CROP545_Phenotype.txt", destfile = "CROPS545_Phenotype.txt")
# Import the GAPIT demo data genotypes
gt_scan <- data.frame(read.table("GAPIT_Tutorial_Data/mdp_numeric.txt", header = T, stringsAsFactors = F, sep = "\t", nrows = 1))
classes <- sapply(gt_scan, class)
genotypes <- data.frame(read.table("GAPIT_Tutorial_Data/mdp_numeric.txt", header = T, row.names = 1, colClasses = classes, stringsAsFactors = F, sep = "\t"))

GM <- read.table("GAPIT_Tutorial_Data/mdp_SNP_information.txt", header = T, stringsAsFactors = F, sep = "\t")
CV <- read.table("CROPS545_Covariates.txt", header = T, stringsAsFactors = F, sep = "\t")
phenotypes <- read.table("CROPS545_Phenotype.txt", header = T, stringsAsFactors = F, sep = "\t")


#Load dependencies
if (!require("pacman")) install.packages("pacman")
pacman::p_load(ggplot2,knitr,gridExtra,kableExtra,data.table)


#source("R/Function_qqplot.R")
#source("R/Function_manhattan_QTN.R")
#source("R/Function_manhattan.R")
#source("R/FUNCTION_powerFDRtype1error.R")

Question1#2#3

Please see GWAS_apply.R for the source code to these questions

Question 4

Demo with QTN

#Create QTNs
source("http://www.zzlab.net/StaGen/2020/R/G2P.R")
source("http://www.zzlab.net/StaGen/2020/R/GWASbyCor.R")
NQTN = 10 #specify number of QTN
h2 = 0.75 #heritability
alpha = 0.6 #alpha value
distribution = "normal" #specify distribution
#set.seed(1337)
#plotting
mySim=G2P(genotypes, h2, alpha, NQTN, distribution)

phenotypes_n_1=phenotypes[,2]
test=GWASapply(mySim$y,genotypes,CV,GM,PCA.M = 3,mySim$QTN.position,plots=TRUE,messages=TRUE)
hist(test$power.fdr.type1error.res$power)

Question 5

Simulated Data

Manhattan Plots

 #phenotype simulation
  p= GWASbyCor(X=genotypes,y=mySim$y)#GWAS by using the correlation methods
  p=data.frame(p)
  pres=power.fdr(p, QTN.position=mySim$QTN.position,cutoff=NULL)
  gcm=manhattan_plot(GM,p,cutoff=NULL,mySim$QTN.position,FP=pres$FP.fdr.power,TP=pres$TP.fdr.power, trait = "unknown",messages=TRUE)
  gcq=qq_plot(GM,p,mySim$QTN.position, trait = "unknown")
  test=GWASapply(mySim$y,genotypes,CV,GM,PCA.M = 3,QTN.position=mySim$QTN.position)
  gtm=manhattan_plot(GM,test$P.value.res,cutoff=NULL,mySim$QTN.position,FP=test$False.Positive,TP=test$True.Positive, trait = "unknown",messages=TRUE)
  gtq=qq_plot(GM,test$P.value.res,mySim$QTN.position, trait = "unknown")
  grid.arrange(gcm, gtm, nrow = 2, ncol = 1, top = "Variance of Principal Components")
  grid.arrange(gcq,gtq, nrow = 2, ncol = 1, top = "Variance of Principal Components")

Mean True and False Positives

gcrep=replicate (30,{
  mySim=G2P(genotypes, h2, alpha, NQTN, distribution) #phenotype simulation
  p= GWASbyCor(X=genotypes,y=mySim$y)#GWAS by using the correlation methods
  p=data.frame(p)
  pres=power.fdr(p, QTN.position=mySim$QTN.position,cutoff=NULL)
  FalP=length(pres$FP.fdr.power)
  TruP=length(pres$TP.fdr.power)
  num_output=c()
  num_output=c(FalP,TruP)
  })
View(gcrep)
print(rowMeans(gcrep))
# 5.16 (# QTN in top 10 p-value SNPs)  211.59 (# QTN with p-value lower than 7th highest)
print(apply(gcrep,1,sd))

garep=replicate(30,{
  mySim=G2P(genotypes, h2, alpha, NQTN, distribution)
  test_rep=GWASapply(mySim$y,genotypes,CV,GM,PCA.M = 3,QTN.position=mySim$QTN.position)
  FalP=length(test_rep$False.Positive)
  TruP=length(test_rep$True.Positive)
  num_output=c()
  num_output=c(FalP,TruP)
  })
print(rowMeans(garep))
# 5.16 (# QTN in top 10 p-value SNPs)  211.59 (# QTN with p-value lower than 7th highest)
print(apply(garep,1,sd))

Question 6

phenotypes_n_1=phenotypes[,2]
loop_start <- proc.time()
test=GWASapply_rapid(phenotypes_n_1,genotypes,CV,GM,PCA.M = 3)
loop_end <- proc.time()
loop_elapsed <- loop_end[3] - loop_start[3]

timet=c()
garep=replicate(30,{
loop_start <- proc.time()
test=GWASapply_rapid(phenotypes_n_1,genotypes,CV,GM,PCA.M = 3)
loop_end <- proc.time()
loop_elapsed <- loop_end[3] - loop_start[3]
timet=c(loop_elapsed)
  })
View(garep)
mean(garep)
manhattan_plot(GM,test$P.value.res,cutoff=NULL, trait = "unknown",messages=TRUE)