---
title: "Lab 9"
output: html_notebook
editor_options: 
  chunk_output_type: console
---

# GWAS by GLM

```{r}
library('MASS') # required for ginv
library(multtest)
library(gplots)
library(compiler) #required for cmpfun
library("scatterplot3d")

source("http://www.zzlab.net/GAPIT/emma.txt")
source("http://www.zzlab.net/GAPIT/gapit_functions.txt")
source("http://zzlab.net/StaGen/2020/R/G2P.R")

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)

X=myGD[,-1]
index1to5=myGM[,2]<6
X1to5 = X[,index1to5]

set.seed(20210324)
mySim=G2P(X=X1to5,h2=.75,alpha=1,NQTN=10,distribution="normal")
myY=cbind(as.data.frame(myGD[,1]), mySim$y)

myGAPIT=GAPIT(
Y=myY,
GD=myGD,
GM=myGM,
QTN.position=mySim$QTN.position,
PCA.total=3,
model="GLM")
```

# GWAS by MLM

```{r}
myGAPIT=GAPIT(
Y=myY,
GD=myGD,
GM=myGM,
QTN.position=mySim$QTN.position,
PCA.total=3,
model="MLM")

```


# GWAS by CMLM

```{r}
myGAPIT=GAPIT(
Y=myY,
GD=myGD,
GM=myGM,
QTN.position=mySim$QTN.position,
PCA.total=3,
group.from=1,
group.to=1000,
group.by=10,
model="CMLM")

```

# GWAS by MLMM

```{r}
# pheno = as.data.frame(cbind(myGD[,1], mySim$y))
# colnames(pheno) = c("taxa","Y1")

# pheno = data.frame(taxa = myGD[,1], Y1 = mySim$y)

myGAPIT_MLM <- GAPIT(
# Y=pheno,
# Y = myY,
Y = mySim$y,
GD=myGD,#Genotype
GM=myGM,#Map information
PCA.total=3,
QTN.position=mySim$QTN.position,
model="MLMM",# Blink
file.output=T
)

```

# GWAS by multiple model

Try by yourself (ref: https://www.zzlab.net/GAPIT/gapit_help_document.pdf)

