---
title: "Lab_6"
author: "Zhou"
date: "2/23/2021"
output: html_document
---

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

## 1. GWAS by correlation

```{r}
# Load 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)
source("http://zzlab.net/StaGen/2020/R/G2P.R")
source("http://zzlab.net/StaGen/2020/R/GWASbyCor.R")
X=myGD[,-1]

# phenotype simulation
set.seed(99164)
mySim=G2P(X= myGD[,-1],h2=.75,alpha=1,NQTN=10,distribution="normal")

# GWAS by correlation
p= GWASbyCor(X=X,y=mySim$y)

# Manhattan plot
color.vector <- rep(c("deepskyblue","orange","forestgreen","indianred3"),10)
m=nrow(myGM)
plot(t(-log10(p))~seq(1:m),col=color.vector[myGM[,2]])

# True QTN position
abline(v=mySim$QTN.position, lty = 2, lwd=2, col = "gray")
mySim$QTN.position

# Top 10 SNPs
index=order(p)
top10=index[1:10]
top10
abline(h=-log10(max(p[top10])), lty = 1, lwd=2, col = "black")

# Detected QTNs within top 10
detected=intersect(top10,mySim$QTN.position)
detected
abline(v= detected, lty = 2, lwd=2, col = "blue")

# False positive within top 10
falsePositive=setdiff(top10, mySim$QTN.position)
falsePositive
abline(v= falsePositive, lty = 2, lwd=2, col = "red")
```

## 2. Bin approach

```{r}
bigNum=1e9
resolution=100000
bin=round((myGM[,2]*bigNum+myGM[,3])/resolution)
result=cbind(myGM,t(p),bin)
head(result)
tail(result)

# Bins of QTNs
QTN.bin=result[mySim$QTN.position,]
QTN.bin

# Sort bins of QTNs
index.qtn.p=order(QTN.bin[,4])
QTN.bin[index.qtn.p,]

```

## 3. Area under Curve (AUC)

### 3.1 Power vs FDR and TypeI error

```{r}
library(compiler) #required for cmpfun
source("http://www.zzlab.net/GAPIT/gapit_functions.txt")
myStat=GAPIT.FDR.TypeI(
# WS- window size
WS=c(1e0,1e3,1e4,1e5), 
#Input: GM - m by 3  matrix for SNP name, chromosome and BP
GM=myGM,
#Input: seqQTN - s by 1 vecter for index of QTN on GM (+1 for GDP column wise)
seqQTN=mySim$QTN.position,
#Input: GWAS - SNP,CHR,BP,P,MAF
GWAS=result)

str(myStat)

# Area Under Curve (AUC)
par(mfrow=c(1,2),mar = c(5,5,5,2))
plot(myStat$FDR[,1],myStat$Power,type="b",  xlab = "FDR", ylab = "Power" )
plot(myStat$TypeI[,1],myStat$Power,type="b", xlab = "Type I error", ylab="Power")

```

### 3.2 Average power vs Average FDR and TypeI error

```{r}
# Replicate 100 times and plot the means
nrep=100
set.seed(99164)
statRep=replicate(nrep, {
mySim=G2P(X=myGD[,-1],h2=.5,alpha=1,NQTN=10,distribution="norm")
p=p= GWASbyCor(X=myGD[,-1],y=mySim$y)
seqQTN=mySim$QTN.position
myGWAS=cbind(myGM,t(p),NA)
myStat=GAPIT.FDR.TypeI(WS=c(1e0,1e3,1e4,1e5), GM=myGM,seqQTN=mySim$QTN.position,GWAS=myGWAS,maxOut=100,MaxBP=1e10)
})
str(statRep)

# Power
power=statRep[[2]]
# FDR
s.fdr=seq(3,length(statRep),7)
fdr=statRep[s.fdr]
fdr.mean=Reduce ("+", fdr) / length(fdr)

# Plots of power vs FDR
theColor=rainbow(4)
# theColor=c("red","blue","green")
plot(fdr.mean[,1],power , type="b", col=theColor [1],xlim=c(0,1))
for(i in 2:ncol(fdr.mean)){
  lines(fdr.mean[,i], power , type="b", col= theColor [i])
}

#AUC
s.auc.fdr=seq(6,length(statRep),7)
auc.fdr=statRep[s.auc.fdr]
auc.fdr.mean=Reduce ("+", auc.fdr) / length(auc.fdr)
# Plots of AUC
barplot(auc.fdr.mean, names.arg=c("1bp", "1K", "10K","100K"), xlab="Resolution", ylab="AUC")

```


