-
Notifications
You must be signed in to change notification settings - Fork 3
/
Debuge.R
executable file
·70 lines (63 loc) · 1.94 KB
/
Debuge.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#Import library
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")
#Import data
setwd("~/Dropbox/gapit/GAPIT_Tutorial_Data")
myGD <- read.table("mdp_numeric.txt", head = TRUE)
myGM <- read.table("mdp_SNP_information.txt" , head = TRUE)
myG <- read.delim("mdp_genotype_test.hmp.txt", head = FALSE)
#Simulate phentype
set.seed(99164)
myPheno=GAPIT.Phenotype.Simulation(GD=myGD,h2=.5,NQTN=10,QTNDist="normal",effectunit=.299)
myY=myPheno$Y
colnames(myY)=c("taxa", "SimPheno")
#Override with modified functions
source('~/Dropbox/gapit/Functions/GAPIT.R', chdir = TRUE)
source('~/Dropbox/gapit/Functions/GAPIT.Main.R', chdir = TRUE)
source('~/Dropbox/gapit/Functions/GAPIT.Manhattan.R', chdir = TRUE)
source('~/Dropbox/gapit/Functions/GAPIT.MAF.R', chdir = TRUE)
setwd("~/Desktop/temp")
#GAPIT on G
myGAPIT <- GAPIT(
Y=myY,
G=myG,
group.from=1,
group.to=1,
group.by=10,
PCA.total=3,
memo="Test"
)
#GAPIT on GD
setwd("~/Desktop/temp")
myGAPIT=GAPIT(
plot.style="Beach",#options: Rainbow, Rushville, FarmCPU, Congress, Ocean, PLINK, Beach
Y=myY,
GD=myGD,
GM=myGM,
QTN.position= myPheno$QTN.position,
PCA.total=3,
group.from=1,
group.to=1,
group.by=10,
file.out=T,
#memo="Beach"
)
#BLINK
setwd("~/Desktop/temp")
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("~/Dropbox/zhanglab/BLINK/mac/blink --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=myPheno$QTN.position)
GAPIT.QQ(myP)