---
title: "Lab_13"
author: "Zhou"
date: "4/20/2021"
output: html_document
---

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

## 1 rrBLUP

### 1.1 Simulate phenotype

```{r}
library(compiler) #required for cmpfun
source("http://www.zzlab.net/GAPIT/gapit_functions.txt")
setwd("/Users/ZhouTang/Downloads/WSU_courses/2021_spring/545/2021/Lab/Lab13")
set.seed(13)

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)
mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=2, effectunit =.95,QTNDist="normal",CV=myCV,cveff=c(.01,.01))

```

### 1.2 Ridge regression & gBLUP

```{r}
#Import rrBLUP
#install.packages("rrBLUP") 
library(rrBLUP)

#prepare data
y <- mySim$Y[,2]
M=as.matrix(X)

#Ridge Regression
ans1 <- mixed.solve(y=y,Z=M)

#gBLUP
# the kinship
K <- tcrossprod(M) #K = MM'
# By introducing the kinship matrix, the random effect is set to be individual effect
ans2 <- mixed.solve(y=y,K=K)

#Compare GEBV
plot(M%*%ans1$u, ans2$u)

```

### 1.3 rrBLUP & GAPIT

```{r}
myGAPIT <- GAPIT(
Y=mySim$Y,
GD=myGD,
GM=myGM,
model="gBLUP",
file.output=F)

order.raw=match(taxa,myGAPIT$Pred[,1])
plot(ans2$u, myGAPIT$Pred[order.raw,5]) 

```