Table of contents

1 Arabidopsis

pheIdx=36
mod=read.table(paste0("./arab/P", pheIdx,".oath.mod"), as.is = T)
info=read.table(paste0("./arab/P", pheIdx,".oath.info"), as.is = T, header = T)
betaG=as.matrix(read.table(paste0("./arab/P", pheIdx,".oath.beta"), as.is = T, header = T))
seG=as.matrix(read.table(paste0("./arab/P", pheIdx,".oath.se"), as.is = T, header = T))
pMatG=-log10((1-pnorm(abs(betaG/seG)))*2)

pRange=matrix(0, nrow(pMatG), 3)
pRange[,1]=apply(pMatG, 1, min)
pRange[,2]=apply(pMatG, 1, max)
pRange[,3]=pRange[,2]-pRange[,1]
idx=which(pRange[,2]> -log10(0.05/nrow(pRange)) & pRange[,3]>3)

if (length(idx)>0) {
  layout(matrix(c(1,1,2),3,1))
  for(i in 1:length(idx)){
    od=order(pMatG[idx[i],])
    barplot(main=info$SNP[idx[i]], names.arg = "", xlab="p-value", pMatG[idx[i],od], border = NA, col=ifelse(pMatG[idx[i],od]> -log10(0.05/nrow(pRange)), "pink", "grey90"))
    abline(h=-log10(0.05/nrow(pMatG)), col="pink", lty=2)
    barplot(names.arg = rep("", nrow(mod)), xlab="", t(mod[od,]), border = NA, col=c("black", "red", "green", "blue", "cyan"))
  }
}

2 Grass

File="./barny/barny_cut_var.oath"
mod=read.table(paste0(File,".mod"), as.is = T)
info=read.table(paste0(File, ".info"), as.is = T, header = T)
betaG=as.matrix(read.table(paste0(File, ".beta"), as.is = T, header = T))
seG=as.matrix(read.table(paste0(File, ".se"), as.is = T, header = T))
pM=pchisq((betaG^2/seG^2), 1, lower.tail = F)
gc=qchisq(apply(pM, 2, median), 1, lower.tail = F)/qchisq(0.5, 1, F)
pMatG=-log10(pchisq((betaG^2/seG^2)/gc, 1, lower.tail = F))

pRange=matrix(0, nrow(pMatG), 3)
pRange[,1]=apply(pMatG, 1, min)
pRange[,2]=apply(pMatG, 1, max)
pRange[,3]=pRange[,2]-pRange[,1]
idx=which(pRange[,2] > -log10(0.05/nrow(pRange)) & pRange[,3]>24)

if (length(idx)>0) {
  layout(matrix(c(1,1,2),3,1))
  for(i in 1:length(idx)) {
    od=order(pMatG[idx[i],])
    barplot(main=info$SNP[idx[i]], names.arg = "", xlab="p-value", pMatG[idx[i],od], border = NA, col=ifelse(pMatG[idx[i],od]> -log10(0.05/nrow(pRange)), "pink", "grey90"))
    abline(h=-log10(0.05/nrow(pMatG)), col="pink", lty=2)
    barplot(names.arg = rep("", nrow(mod)), xlab="", t(mod[od,]), border = NA, col=c("red", "blue", "pink", "yellow", "black"))
  }
}