gap (version 1.1-1)

mhtplot2: Manhattan plot with annotations

Description

To generate Manhattan plot with annotations. The function is generic and for instance could be used for genomewide p values or any random variable that is uniformly distributed. By default, a log10-transformation is applied. Note that with real chromosomal positions, it is also appropriate to plot and some but not all chromosomes.

It is possible to specify options such as xlab, ylim and font family when the plot is requested for data in other context.

To maintain back compatibility options as in mhtplot are used. The positions of the horizontal labels are now in the middle rather than at the beginning of their bands in the plot.

Usage

mhtplot2(data, control=mht.control(), hcontrol=hmht.control(), ...)

Arguments

data

a data frame with three columns representing chromosome, position and p values

control

A control function named mht.control() with the following arguments,

1. type. a flag with value "p" or "l" indicating if points or lines are to be drawn.
2. usepos. a flag to use real chromosomal positions as composed to ordinal positions wi

hcontrol

A control function named hmht.control() with the following arguments,

1. data. chunk of data to be highlighted with default value NULL
2. colors. colors for annotated genes
3. yoffset. offset above the data point showing most significan

...

other options in compatible with the R plot function

Value

• The plot is shown on or saved to the appropriate device.

Examples

# The following example uses only chromosomes 14 and 20 of the Nat Genet paper.
#
mdata <- within(hr1420,{
  c1<-colour==1
  c2<-colour==2
  c3<-colour==3
  colour[c1] <- 62
  colour[c2] <- 73
  colour[c3] <- 552
})
mdata <- mdata[,c("CHR","POS","P","gene","colour")]
ops <- mht.control(colors=rep(c("lightgray","gray"),11),yline=1.5,xline=2,srt=0)
hops <- hmht.control(data=subset(mdata,!is.na(gene)))
v <- "Verdana"
ifelse(Sys.info()['sysname']=="Windows", windowsFonts(ffamily=windowsFont(v)), ffamily <- v)
tiff("mh.tiff",width=.03937*189,height=.03937*189/2,units="in",res=1200,compress="lzw")
par(las=2, xpd=TRUE, cex.axis=1.8, cex=0.4)
mhtplot2(with(mdata,cbind(CHR,POS,P,colour)),ops,hops,pch=19,
         ylab=expression(paste(plain("-"),log[10],plain("p-value"),sep=" ")),
         family="ffamily")
axis(2,pos=2,at=seq(0,25,5),family="ffamily",cex=0.5,cex.axis=1.1)
dev.off()