---
title: "Lab12"
author: "Zhou"
date: "4/13/2021"
output: html_document
---

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

## 1. MAS

## 1.1 Simulate phenotype

```{r}
set.seed(12)
setwd("/Users/ZhouTang/Downloads/WSU_courses/2021_spring/545/2021/Lab/Lab12")
library(compiler) #required for cmpfun
source("http://www.zzlab.net/GAPIT/gapit_functions.txt")
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)

#Simultate 10 QTN on the first half chromosomes
X=myGD[,-1]
index1to5=myGM[,2]<6
X1to5 = X[,index1to5]
taxa=myGD[,1]
GD.candidate=cbind(taxa,X1to5)
nqtn=5

# Phenotype
mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=nqtn, effectunit =.95,QTNDist="normal",CV=myCV,cveff=c(.01,.01))

# PCA
pcs = prcomp(myGD[,-1])
pc3 = pcs$x[,1:3]
```

## 1.2 Split sample into training and testing parts 

20% samples as testing, 80% samples as training

```{r}
n=nrow(mySim$Y)
# 20% samples as testing 
testing=sample(n,round(n/5),replace=F)
training=-testing

```

## 1.3 GWAS on training set

```{r}
#GWAS on training to find QTNs
myGAPIT_train <- GAPIT(
Y=mySim$Y[training,],
GD=myGD,GM=myGM,CV=myCV,
PCA.total=3,model="GLM",QTN.position=mySim$QTN.position, 
memo="GWAS_train_2")
```

## 1.4 Prediction with MAS

```{r}
ntop= 3
index=order(myGAPIT_train$GWAS[,4])
top=index[1:ntop]
myQTN=cbind(myCV[,1:3], pc3, myGD[,c(top+1)])

#Estimate effect and prediction
myGAPIT_pre <- GAPIT(
Y=mySim$Y[training,],
GD=myGD,GM=myGM,
CV=myQTN,model="GLM", 
SNP.test=FALSE,
memo="GLM_testing",)

```

## 1.5 Accuracy in training set

```{r}
order=match(mySim$Y[,1],myGAPIT_pre$Pred[,1])
myPred=myGAPIT_pre$Pred[order,]

ry2=cor(myPred[training,8],mySim$Y[training,2])^2
ru2=cor(myPred[training,8],mySim$u[training])^2
par(mfrow=c(2,1), mar = c(3,4,1,1))
plot(myPred[training,8],mySim$Y[training,2])
mtext(paste("R square=",ry2,sep=""), side = 3)
plot(myPred[training,8],mySim$u[training])
mtext(paste("R square=",ru2,sep=""), side = 3)

```

## 1.6 Accuracy in testing set

```{r}
ry2=cor(myPred[testing,8],mySim$Y[testing,2])^2
ru2=cor(myPred[testing,8],mySim$u[testing])^2
par(mfrow=c(2,1), mar = c(3,4,1,1))
plot(myPred[testing,8],mySim$Y[testing,2])
mtext(paste("R square=",ry2,sep=""), side = 3)
plot(myPred[testing,8],mySim$u[testing])
mtext(paste("R square=",ru2,sep=""), side = 3)

```

## 2. gBLUP by GAPIT

```{r}
myGAPIT_gBLUP <- GAPIT(
  Y=mySim$Y[training,],
  GD=myGD,
  GM=myGM,
  PCA.total=3,
  CV=myCV,
  model="gBLUP",
  SNP.test=FALSE,
  memo="gBLUP")

```

## 2.1 Accuracy in training set

```{r}
order=match(mySim$Y[,1],myGAPIT_gBLUP$Pred[,1])
myPred=myGAPIT_gBLUP$Pred[order,]

ry2=cor(myPred[training,8],mySim$Y[training,2])^2
ru2=cor(myPred[training,8],mySim$u[training])^2
ry2.blup=cor(myPred[training,5],mySim$Y[training,2])^2
ru2.blup=cor(myPred[training,5],mySim$u[training])^2
par(mfrow=c(2,2), mar = c(3,4,1,1))
plot(myPred[training,8],mySim$Y[training,2])
mtext(paste("R square=",ry2,sep=""), side = 3)
plot(myPred[training,8],mySim$u[training])
mtext(paste("R square=",ru2,sep=""), side = 3)
plot(myPred[training,5],mySim$Y[training,2])
mtext(paste("R square=",ry2.blup,sep=""), side = 3)
plot(myPred[training,5],mySim$u[training])
mtext(paste("R square=",ru2.blup,sep=""), side = 3)

```

## 2.2 Accuracy in testing set

```{r}
order=match(mySim$Y[,1],myGAPIT5$Pred[,1])
myPred=myGAPIT5$Pred[order,]

ry2=cor(myPred[testing,8],mySim$Y[testing,2])^2
ru2=cor(myPred[testing,8],mySim$u[testing])^2
ry2.blup=cor(myPred[testing,5],mySim$Y[testing,2])^2
ru2.blup=cor(myPred[testing,5],mySim$u[testing])^2
par(mfrow=c(2,2), mar = c(3,4,1,1))
plot(myPred[testing,8],mySim$Y[testing,2])
mtext(paste("R square=",ry2,sep=""), side = 3)
plot(myPred[testing,8],mySim$u[testing])
mtext(paste("R square=",ru2,sep=""), side = 3)
plot(myPred[testing,5],mySim$Y[testing,2])
mtext(paste("R square=",ry2.blup,sep=""), side = 3)
plot(myPred[testing,5],mySim$u[testing])
mtext(paste("R square=",ru2.blup,sep=""), side = 3)

```