Alternate color manhattan plot

I am trying to make my own manhattan plot the way I want for a long time, and this one seems the be the shortest version that I come up with. Here’s the feature of my function.

1. Should be transportable. Only requires “CHR”, “BP”, “P” to run at a minimum.

2. You can choose to draw your own reference lines, using — cutoffs

3. You can define how tall your graph will be for comparison, –using maxy

4. It is based on basic plot function. So, you can add title to your plot later on. Add more color points, using points [need to convert your coordinate to the same way first].

The down sides.

1. I use the order on chromosome, instead of chromosomal position. So, it is not good to use this code if you only have one gene, and want to the signal location on the chromosome.

2. The signals are overlapping between chromosome in the low p-value region. However, this does not obscure your signal in the top. Improvement can be added.

3. Not all chromosomes are labelled on the X-axis. This is probably the most annoying things that need to be improved.

4. It doesn’t have much advantage in terms of speed. It still took 83 seconds to plot 860,000 signals. That’s roughly 10,000 points per second, though.

Besides all these it still seems pretty good.

manhattan <- function(data,color=c("black","grey"),maxy,cutoffs= c(5e-8)) {
if (! sum(c("CHR","BP","P") %in% names(data)) == 3) stop("Data must have CHR BP P column")
data <- data[order(data$CHR,data$BP),]
data <- data[!is.na(data$P),]
nsnps <- nrow(data)
data$COOR <- 1:nsnps
ticks <- tapply(data$COOR,data$CHR,median)
data$nlogp <- -log10(data$P)
chr <- sort(unique(data$CHR))
par(mar=c(5,4,4,2))
if (missing(maxy)) {
	plot(1:nsnps,data$nlogp,
	xaxt="n",type="n",xlab="chromosome",ylab="-log10(p)")
	} else {
	plot(1:nsnps,data$nlogp,ylim=c(0,maxy),
	xaxt="n",type="n",xlab="chromosome",ylab="-log10(p)")	}
axis(1,at=ticks,labels=chr)
plot.select <- ifelse(data$CHR%%2,1,0)
with(subset(data,plot.select==1),points(COOR,nlogp,pch=16,col=color[1]))
with(subset(data,plot.select==0),points(COOR,nlogp,pch=16,col=color[2]))
for (coi in 1:length(cutoffs)) {
	abline(h=-log10(cutoffs[coi]),lty="dashed",col="red")
}
}