---
title: "Lab14"
author: "Zhou"
date: "4/28/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Data and simulation

```{r}
setwd("/Users/ZhouTang/Downloads/WSU_courses/2021_spring/545/2021/Lab/Lab14")
source("http://zzlab.net/GAPIT/GAPIT.library.R")
source("http://zzlab.net/GAPIT/gapit_functions.txt")

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

X=myGD[,-1]
taxa=myGD[,1]
index1to5=myGM[,2]<6
X1to5 = X[,index1to5]
GD.candidate=cbind(as.data.frame(taxa),X1to5)

# simulation of phenotype
set.seed(99164)
mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=20, effectunit =.95,QTNDist="normal",CV=myCV,cveff=c(.2,.2),a2=.5,adim=3,category=1,r=.4)
myY=mySim$Y
colnames(myY)=c("Taxa","Trait1")

# Training set and testing set
nfold=5
# Another way to sample data
sets=sample(cut(1:nrow(myY),nfold,labels=FALSE),nrow(myY))
training=myY[,c(1,2)]
test_index=sets==1
training[test_index,2]=NA
testing=cbind(myY[test_index,c(1,2)], mySim$u[test_index])

```

## gBLUP through GAPIT

```{r}
myGAPIT5 <- GAPIT(
  Y=training,
  GD=myGD,
  GM=myGM,
  PCA.total=3,
  CV=myCV,
  model="gBLUP",
  file.output=F
  )

# Merge data here.
gapit5=merge(myY,myGAPIT5$Pred[,c(1,3,5,8)],by.x="Taxa",by.y="Taxa")
# Accuracy on training set
fit.GAPIT=cor(gapit5[gapit5[,3]==1,2],gapit5[gapit5[,3]==1,5]) 
fit.GAPIT
# Accuracy on testing set
acc.GAPIT=cor(gapit5[gapit5[,3]!=1,2],gapit5[gapit5[,3]!=1,5]) 
acc.GAPIT

```

## BayesA, B and Cpi through BAGS

```{r}
#Prepare BAGS
source('http://zzlab.net/sandbox/BAGS.R')
niter=200
m=ncol(X)
# training set
GR=myGD[!test_index,-1];YR=as.matrix(mySim$Y[!test_index,2])
# testing set
GI=myGD[test_index,-1];YI=as.matrix(mySim$Y[test_index,2])

#Bayes A:
myBayes_a=BAGS(X=GR,y=YR,pi=0,burn.in=100,burn.out=100,recording=T)

#Bayes B:
myBayes_b=BAGS(X=GR,y=YR,pi=.95,burn.in=100,burn.out=100,recording=T)

#Bayes Cpi:
myBayes_cpi=BAGS(X=GR,y=YR,pi=1,burn.in=100,burn.out=100,recording=T)

```

## BLR

```{r}
#install.packages("BLR")
library(BLR)
nIter=2000              #### number of iteration
burnIn=1500             #### burn in a part of iteration

PCA=myGAPIT5$PCA
pcEnv=PCA[,-1]

# Bayeian LASSO
myBLR =BLR(y=as.matrix(myY[!test_index,2]),
	XF=pcEnv[! test_index,-1],
	XL=as.matrix(myGD[! test_index,-1]),
	nIter=nIter,
	burnIn=burnIn)

# Accuracy on training 
pred.ref=as.matrix(myGD[!test_index,-1])%*%myBLR$bL
modelfit <- cor(myY[! test_index,2],pred.ref)
modelfit

# Accuracy on testing
pred.inf=as.matrix(myGD[test_index,-1])%*%myBLR$bL
accuracy <- cor(myY[test_index,2],pred.inf)
accuracy

```

## BGLR

```{r}
#install.packages("BGLR")
library(BGLR)
nIter=2000              #### number of iteration
burnIn=1500             #### burnin a part of iteration
myBGLR =BGLR(y=as.matrix(myY[! test_index,2]),
ETA=list(list(X=pcEnv[! test_index,-1],model='FIXED', saveEffects=TRUE),
	list(X=as.matrix(myGD[! test_index,-1]),model='BRR', saveEffects=TRUE)
	#list(X=as.matrix(myGD[! test_index,-1]),model='BL', saveEffects=TRUE)
	#list(X=as.matrix(myGD[! test_index,-1]),model='BayesA', saveEffects=TRUE)
	#list(X=as.matrix(myGD[! test_index,-1]),model='BayesB', saveEffects=TRUE)
	),
	nIter=nIter,
	burnIn=burnIn)

# Accuracy on training
pred.ref=as.matrix(myGD[! test_index,-1])%*%myBGLR$ETA[[2]]$b
modelfit <- cor(myY[! test_index,2],pred.ref)
modelfit

# Accuracy on testing
pred.inf=as.matrix(myGD[test_index,-1])%*%myBGLR$ETA[[2]]$b
accuracy <- cor(myY[test_index,2],pred.inf)
accuracy

```

