if(!require("ggplot2")) {
print("Please install the ggplot plackage")
## install.packages('ggplot2', dep=T)
}
require(ggplot2)
d <- read.table("simulated_assoc.assoc",header=T)
d <- na.omit(d)
maxy <- 15
sig <- -log10(5e-8)
d$logp <- -log10(d$P)
d$x <- NA
ticks <- NULL
lastbase <- 0
order <- order(d$CHR,d$BP)
d <- d[order,]
d[d$logp > maxy,'logp'] <- maxy
numchrs <- length(unique(d$CHR))
for( i in unique(d$CHR) ) {
if( i==1 ) {
d[d$CHR==i,]$x <- d[d$CHR==i,]$BP
} else {
lastbase <- lastbase+tail(subset(d,CHR==i-1)$BP,1)
d[d$CHR==i,]$x <- d[d$CHR==i,]$BP+lastbase
}
ticks <- c(ticks,d[d$CHR==i,]$x[floor(length(d[d$CHR==i,]$x)/2)+1])
}
ticklim <- c(min(d$x),max(d$x))
if(numchrs > 1) {
cols <- rep(c("gray","black"),max(d$CHR))
} else {
cols <- "gray"
}
if (maxy=="max") maxy=ceiling(max(d$logp)) else maxy=maxy
if (maxy<8) maxy=8
pdf("manhattan_plot.pdf")
plot <- qplot(x,logp,data=d,
ylab=expression(-log[10](italic(p))) , colour=factor(CHR))
plot <- plot + scale_x_continuous(name="Chromosome",
breaks=ticks,
labels=(unique(d$CHR)))
plot <- plot + scale_y_continuous(limits=c(0,maxy),
breaks=1:maxy)
plot <- plot + scale_colour_manual(values=cols)
plot <- plot + theme(legend.position="none")
plot <- plot + labs(title="Manhattan Plot")
plot <- plot + theme(panel.background=element_blank(),
panel.grid.minor=element_blank())
plot <- plot + geom_abline(intercept=sig,slope=0, col="black")
