---
title: "Lab10"
author: "Zhou"
date: "3/30/2021"
output: html_document
---

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

## FarmCPU

```{r}
library(compiler) #required for cmpfun
source("http://www.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)

#Simultate 10 QTN on the first half chromosomes
X=myGD[,-1]
index1to5=myGM[,2]<6
X1to5 = X[,index1to5]
taxa=myGD[,1]
set.seed(20210331)
GD.candidate=cbind(taxa,X1to5)
mySim=GAPIT.Phenotype.Simulation(GD=GD.candidate,GM=myGM[index1to5,],h2=.5,NQTN=10, effectunit =.95,QTNDist="normal")

```

```{r}
#install.packages("bigmemory")
#install.packages("biganalytics")

library(bigmemory)
library(biganalytics)
source("http://www.zzlab.net/FarmCPU/FarmCPU_functions.txt")
setwd("/Users/ZhouTang/Downloads/WSU_courses/2021_spring/545/2021/Lab/Lab10")
myFarmCPU= FarmCPU(
    Y=mySim$Y,
    GD=myGD,
    GM=myGM ) 

myP=as.numeric(myFarmCPU$GWAS$P.value)
myGI.MP=cbind(myGM[,-1],myP)
GAPIT.Manhattan(GI.MP=myGI.MP,seqQTN=mySim$QTN.position)
GAPIT.QQ(myP)

```

### FarmCPU via GAPIT

```{r}
# FarmCPU through GAPIT
myGAPIT_MLM <- GAPIT(
Y=mySim$Y,
GD=myGD,#Genotype
GM=myGM,#Map information
PCA.total=3,
QTN.position=mySim$QTN.position,
model="FarmCPU",# Can choose MLM CMLM GLM SUPER MLMM FarmCPU Blink
file.output=T)

```


## Blink

### Blink C version

```{r}
myY=mySim$Y
colnames(myY)=c("taxa", "SimPheno")
myGD1=t(myGD[,-1])
write.table(myY,file="myData.txt",quote=F,row.name=F,col.name=T,sep="\t")
write.table(myGD1,file="myData.dat",quote=F,row.name=F,col.name=F,sep="\t")
write.table(myGM,file="myData.map",quote=F,row.name=F,col.name=T,sep="\t")
#run blink
system("./blink_mac --file myData --max_loop 10 --numeric --gwas")
#Extract p values
result<- read.table("SimPheno_GWAS_result.txt", head = TRUE)
myP=as.numeric(result[,5])
myGI.MP=cbind(myGM[,-1],myP)
GAPIT.Manhattan(GI.MP=myGI.MP,seqQTN=mySim$QTN.position)
GAPIT.QQ(myP)

```

### BLINK via GAPIT

```{r}
myGAPIT_MLM <- GAPIT(
Y=mySim$Y,
GD=myGD,#Genotype
GM=myGM,#Map information
PCA.total=3,
QTN.position=mySim$QTN.position,
model="Blink",# Can choose MLM CMLM GLM SUPER MLMM FarmCPU Blink
file.output=T
)

```