VignetteGWASBAYES.R
In GWAS.BAYES: GWAS for Selfing Species

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  tidy = FALSE
)

## ---- warning=FALSE,message=FALSE---------------------------------------------
library(GWAS.BAYES)

## -----------------------------------------------------------------------------
data("vignette_lm_dat")
head(vignette_lm_dat[,1:10])

## -----------------------------------------------------------------------------
Y <- vignette_lm_dat$Phenotype
SNPs <- vignette_lm_dat[,-1]

## -----------------------------------------------------------------------------
SNPs <- standardize(SNPs = SNPs,method = "major-minor",number_cores = 1)
SNPs[1:6,1:10]

## -----------------------------------------------------------------------------
aggregate_list <- aggregate_SNPs(SNPs = SNPs, Y = Y)
SNPs <- aggregate_list$SNPs
Y <- aggregate_list$Y

## -----------------------------------------------------------------------------
level_list <- level_function(SNPs = SNPs, MAF = 0.01)
SNPs <- level_list$SNPs
level_list$SNPs_Dropped
dim(SNPs)

## -----------------------------------------------------------------------------
fullPreprocess <- preprocess_SNPs(SNPs = vignette_lm_dat[,-1],
  Y = vignette_lm_dat$Phenotype,MAF = 0.01,number_cores = 1,
  na.rm = FALSE)
all.equal(fullPreprocess$SNPs,SNPs)
all.equal(fullPreprocess$Y,Y)
fullPreprocess$SNPs_Dropped;level_list$SNPs_Dropped

## ----fig.cap="Percent Variation in the Response explained by the principal components",fig.wide = TRUE, fig.align="center", fig.width=6.75, fig.height=4.75----
principal_comp <- pca_function(SNPs = SNPs,number_components = 3,
  plot_it = TRUE)

## -----------------------------------------------------------------------------
principal_comp <- as.matrix(principal_comp[,1])

## -----------------------------------------------------------------------------
Significant_SNPs_Bonf <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, 
  principal_components = principal_comp,frequentist = TRUE,
  controlrate = "bonferroni",threshold = .05,kinship = FALSE, 
  info = FALSE) #Bonferroni Correction
  sum(Significant_SNPs_Bonf$Significant) #Five Significant SNPs
  which(Significant_SNPs_Bonf$Significant == 1)

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75-----
resids_diag(Y = Y,SNPs = SNPs,
  significant = Significant_SNPs_Bonf$Significant,
  kinship = FALSE,principal_components = principal_comp,
  plot_it = TRUE)

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75-----
cor_plot(SNPs = SNPs,
  significant = Significant_SNPs_Bonf$Significant,
  info = FALSE)

## -----------------------------------------------------------------------------
Significant_SNPs_FDR <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, 
  principal_components = principal_comp,frequentist = TRUE,
  controlrate = "BH",threshold = .05,kinship = FALSE,info = FALSE) 
sum(Significant_SNPs_FDR$Significant) #Five Significant SNPs
which(Significant_SNPs_FDR$Significant == 1)

## -----------------------------------------------------------------------------
Significant_SNPs_BFDR <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, 
  principal_components = principal_comp,frequentist = FALSE,nullprob = .5,
  alterprob = .5,threshold = .05,kinship = FALSE,info = FALSE)
sum(Significant_SNPs_BFDR$Significant) #Five Significant SNPs
which(Significant_SNPs_BFDR$Significant == 1)

## -----------------------------------------------------------------------------
GA_results <- postGWAS(Y = Y,SNPs = SNPs,number_cores = 1, 
  significant = Significant_SNPs_Bonf$Significant,
  principal_components = principal_comp,maxiterations = 100,
  runs_til_stop = 10,kinship = FALSE,info = FALSE)
GA_results

## -----------------------------------------------------------------------------
data("vignette_kinship_dat")
head(vignette_kinship_dat[,1:10])

## -----------------------------------------------------------------------------
Y <- vignette_kinship_dat$Phenotype
SNPs <- vignette_kinship_dat[,-1]
fullPreprocess <- preprocess_SNPs(SNPs = SNPs, Y = Y, MAF = 0.01,
  number_cores = 1,na.rm = FALSE)
SNPs <- fullPreprocess$SNPs
Y <- fullPreprocess$Y
fullPreprocess$SNPs_Dropped

## -----------------------------------------------------------------------------
library(rrBLUP,quietly = TRUE)
k <- A.mat(SNPs,n.core = 1)
dim(k)

## -----------------------------------------------------------------------------
Significant_SNPs_Bonf <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, 
  principal_components = FALSE,frequentist = TRUE,controlrate = "bonferroni",
  threshold = .05,kinship = k, info = FALSE)
sum(Significant_SNPs_Bonf$Significant)#Four Significant SNPs
which(Significant_SNPs_Bonf$Significant == 1)

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75-----
resids_diag(Y = Y,SNPs = SNPs,significant = Significant_SNPs_Bonf$Significant,
  kinship = k,principal_components = FALSE,plot_it = TRUE)

## -----------------------------------------------------------------------------
Significant_SNPs_FDR <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, 
  principal_components = FALSE,frequentist = TRUE,controlrate = "BH",
  threshold = .05,kinship = k, info = FALSE)
sum(Significant_SNPs_FDR$Significant) #Six Significant SNPs
which(Significant_SNPs_FDR$Significant == 1)

## -----------------------------------------------------------------------------
Significant_SNPs_BFDR <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, 
  principal_components = FALSE,frequentist = FALSE,nullprob = .5,
  alterprob = .5,threshold = .05,kinship = k, info = FALSE)
sum(Significant_SNPs_BFDR$Significant) #4 Significant SNPs
which(Significant_SNPs_BFDR$Significant == 1)

## -----------------------------------------------------------------------------
GA_results <- postGWAS(Y = Y,SNPs = SNPs,number_cores = 1, 
  significant = Significant_SNPs_FDR$Significant,
  principal_components = FALSE,kinship = k, info = FALSE)
GA_results

## -----------------------------------------------------------------------------
data("RealDataSNPs_Y")

Y <- RealDataSNPs_Y$Phenotype
SNPs <- subset(RealDataSNPs_Y,select = -c(Phenotype))

fullPreprocess <- preprocess_SNPs(SNPs = SNPs,Y = Y,
  MAF = 0.01,number_cores = 1,na.rm = FALSE)
SNPs <- fullPreprocess$SNPs
Y <- fullPreprocess$Y
fullPreprocess$SNPs_Dropped

## -----------------------------------------------------------------------------
data("RealDataInfo")
head(RealDataInfo[,1:6])

## -----------------------------------------------------------------------------
RealDataInfo <- RealDataInfo[,-fullPreprocess$SNPs_Dropped]

## -----------------------------------------------------------------------------
data("RealDataKinship")
kinship <- as.matrix(RealDataKinship)

## -----------------------------------------------------------------------------
Significant_SNPs <- preselection(Y = Y, SNPs = SNPs,number_cores = 1, 
  principal_components = FALSE,frequentist = TRUE,controlrate = "bonferroni",
  threshold = .05,kinship = kinship, info = RealDataInfo)
sum(Significant_SNPs$Significant)#11 Significant SNPs
Significant_SNPs[Significant_SNPs$Significant == 1,c(1,2)]

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75-----
resids_diag(Y = Y, SNPs = SNPs, significant = Significant_SNPs$Significant,
  kinship = kinship)

## -----------------------------------------------------------------------------
Significant_SNPs <- preselection(Y = log(Y), SNPs = SNPs,number_cores = 1, 
  principal_components = FALSE,frequentist = TRUE,controlrate = "bonferroni",
  threshold = .05,kinship = kinship,info = RealDataInfo)
sum(Significant_SNPs$Significant)#4 Significant SNPs
Significant_SNPs[Significant_SNPs$Significant == 1,c(1,2)]

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75-----
resids_diag(Y = log(Y), SNPs = SNPs, 
  significant = Significant_SNPs$Significant,kinship = kinship)

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75,warning=FALSE,message=FALSE----
library(qqman,quietly = TRUE) 
manhattan(Significant_SNPs,chr = "Chromosomes",bp = "Positions",
  p = "P_values",suggestiveline = FALSE,
  genomewideline = -log10(.05/nrow(Significant_SNPs)))

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75-----
qq(Significant_SNPs$P_values)

## -----------------------------------------------------------------------------
GA_results <- postGWAS(Y = log(Y),SNPs = SNPs,number_cores = 1, 
  significant = Significant_SNPs$Significant,principal_components = FALSE,
  kinship = kinship,info = RealDataInfo)
GA_results

## ---- fig.wide = TRUE,fig.align="center", fig.width=6.75, fig.height=4.75-----
cor_plot(SNPs = SNPs[,Significant_SNPs$Significant == 1],
  significant = c(1,1,1,1),
  info = RealDataInfo[,Significant_SNPs$Significant == 1])

## -----------------------------------------------------------------------------
1 - colMeans(SNPs[,Significant_SNPs$Significant == 1])

## -----------------------------------------------------------------------------
which(Significant_SNPs$Significant == 1)

## -----------------------------------------------------------------------------
summary(1 - colMeans(SNPs[,801:852]))