As the doctor gone rogue

May 17, 2010

Alternate color manhattan plot

Filed under: genetics, R — hypotheses @ 8:31 pm

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")
}
}
Advertisements

1 Comment »

  1. Hi sorry is it possible to post sample usage and sample data ?

    I got up to stuck up to this point

    > plot.select with(subset(data,plot.select==1),points(COOR,nlogp,pch=16,col=color[1]))
    Error in plot.xy(xy.coords(x, y), type = type, …) :
    object ‘color’ not found
    Calls: with … eval -> eval -> points -> points.default -> plot.xy
    Execution halted

    Comment by kevin — May 21, 2012 @ 12:49 am


RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Blog at WordPress.com.

%d bloggers like this: