args = commandArgs(trailingOnly=TRUE)
library(qqman)
RUNNAME <- args[1]
THRESHOLD <- args[2]
reportListN <- paste("RAiSD_ReportList.", RUNNAME,".txt", sep="")
reportList <- read.table(reportListN)[,1]
reportList
data <- data.frame(pos=c(), value=c(), chr=c())
for (i in 1:length(reportList[])) {
d <- read.table(paste(reportList[i]), header=F, skip=0)[,1:7]
tmp.dat <- data.frame(pos=d[,1], value=d[,7]*1, chr=rep(i, length(d[,1])))
data <- rbind(data, tmp.dat)}
output_file <- paste("RAiSD_ManhattanPlot.", RUNNAME,".pdf", sep="")
pdf(output_file)
snp <- 1:dim(data)[1]
mydf <- data.frame(snp, data)
thres<-as.numeric(THRESHOLD)
topQ<-thres*100
threshold <- quantile(x=data$value, probs = thres)
title_msg <- paste("Manhattan plot for ", RUNNAME,"
threshold=",threshold," (", topQ,"%)")
manhattan(mydf, chr="chr", bp="pos", cex = 0.5, p="value", snp="snp", logp=F, ylab=bquote(mu ~ "statistic"), genomewideline = FALSE, suggestiveline = FALSE,pos=0, col=c("blue2", "darkorange1"), ylim=c(0.0, max(data$value, na.rm = TRUE)*1.4))
title(title_msg)
abline(h=threshold, col="red", lw=2)
dev.off()