## [1] "Document was last updated at 2018-12-08."
Table of contents
NSS
###1 Generate nss
Dir=list.dirs()
Dir=Dir[-1]
for(i in 1:length(Dir)) {
setwd(Dir[i])
bed=list.files(pattern = ".bed")
datIn=strsplit(bed, ".bed")[[1]]
datOut = paste0(datIn, "_ht")
phe=list.files(pattern = ".phe")
covar=list.files(pattern = ".covar")
nss=paste(gear, " nss --bfile ", datIn,
" --pheno ", phe, " --mpheno 1 ",
" --covar ", covar, " --covar-number 1 2 3 4 5",
" --out ", datOut)
print(nss)
system(nss)
setwd(CurDir) #return back
}
## [1] "java -jar /Users/gc5k/Documents/workspace/FromSVN/GEAR/gear.jar nss --bfile Edinburgh500 --pheno Edinburgh500.phe --mpheno 1 --covar Edinburgh500.covar --covar-number 1 2 3 4 5 --out Edinburgh500_ht"
## [1] "java -jar /Users/gc5k/Documents/workspace/FromSVN/GEAR/gear.jar nss --bfile Oxford500 --pheno Oxford500.phe --mpheno 1 --covar Oxford500.covar --covar-number 1 2 3 4 5 --out Oxford500_ht"
PPSR

Projected PC
pp1kg=paste(gear, " profile --bfile ../1KG_chr22/1kg_chr22 --score ../1KG_chr22/1kg_PPSR.score --out 1kg")
system(pp1kg)
ppOx=paste(gear, " profile --bfile Oxford/Oxford500 --score ../1KG_chr22/1kg_PPSR.score --out Oxford500")
system(ppOx)
ppEd=paste(gear, " profile --bfile Edinburgh/Edinburgh500 --score ../1KG_chr22/1kg_PPSR.score --out Edinburgh500")
system(ppEd)
refPC=read.table("1kg.profile", as.is = T, header = T)
oxPC=read.table("Oxford500.profile", as.is = T, header = T)
edPC=read.table("Edinburgh500.profile", as.is = T, header = T)
plot(refPC[,3], refPC[,4], col="grey", pch=16, cex=0.5, bty='n', xlab="E 1", ylab="E 2")
points(oxPC[,3], oxPC[,4], col="red", pch=16, cex=0.5)
points(edPC[,3], edPC[,4], col="green", pch=16, cex=1)
legend("topright", legend = c("Ox", "Ed"), col=c("red", "green"), pch=16, cex=0.5, bty ='n')

meta\(F_{st}\)
mPCA
mpc=paste(gear, "mpc --meta-batch mpc.txt --out mpc")
system(mpc)
layout(matrix(1:2, 1, 2))
mval=read.table("mpc.mval", as.is = T)
barplot(mval[,1], border = F)
mvec=read.table("mpc.mvec", as.is = T)
plot(mvec[,1], mvec[,2], bty='n', xlab="mpc 1", ylab="mpc 2", pch=16)

OATH
Dir=list.dirs()
Dir=Dir[-1]
#parameter [need further automation]
Dat=c("Edinburgh500", "Oxford500")
Size=c(500, 500)
model=c("1", "2","3")
model_code="1_2_3"
for(i in 1:length(Dir)) {
setwd(Dir[i])
oath=paste0(gear, " oath --cm ", Dat[i], "_ht.m.nss ",
" --nss-gz-batch ", Dat[i], "_ht.list.nss ",
" -n ", Size[i],
" --keep-nss ", paste(model, collapse = " "),
" --out ", Dat[i], "_", paste(model, collapse = "_"))
print(oath)
system(oath)
setwd(CurDir) #return back
}
## [1] "java -jar /Users/gc5k/Documents/workspace/FromSVN/GEAR/gear.jar oath --cm Edinburgh500_ht.m.nss --nss-gz-batch Edinburgh500_ht.list.nss -n 500 --keep-nss 1 2 3 --out Edinburgh500_1_2_3"
## [1] "java -jar /Users/gc5k/Documents/workspace/FromSVN/GEAR/gear.jar oath --cm Oxford500_ht.m.nss --nss-gz-batch Oxford500_ht.list.nss -n 500 --keep-nss 1 2 3 --out Oxford500_1_2_3"
## [1] "java -jar /Users/gc5k/Documents/workspace/FromSVN/GEAR/gear.jar oath --cm NA_ht.m.nss --nss-gz-batch NA_ht.list.nss -n NA --keep-nss 1 2 3 --out NA_1_2_3"
## [1] "java -jar /Users/gc5k/Documents/workspace/FromSVN/GEAR/gear.jar oath --cm NA_ht.m.nss --nss-gz-batch NA_ht.list.nss -n NA --keep-nss 1 2 3 --out NA_1_2_3"
GWAMA

