QQ plot for a vector of pvalues

qqplot_pvalue.r

qq = function(pvector, title="Quantile-quantile plot of p-values",red=NULL) {
  # adapted from http://www.gettinggeneticsdone.com/2009/11/qq-plots-of-p-values-in-r-using-ggplot2.html
  pvector[pvector==0]=min(pvector[pvector>0])
  obs <- -log10(pvector)
  N <- length(obs) ## number of p-values
  null <- -log(1:N/N,10)
  MAX <- max(c(obs,null))
  c95 <- rep(0,N)
  c05 <- rep(0,N)
  for(i in 1:N){
    c95[i] <- qbeta(0.95,i,N-i+1)
    c05[i] <- qbeta(0.05,i,N-i+1)
  }
  plot(null, -log(c95,10), ylim=c(0,MAX), xlim=c(0,MAX), type="l", 
       axes=FALSE, xlab="", ylab="",lty="dashed")
  par(new=T)
  plot(null, -log(c05,10), ylim=c(0,MAX), xlim=c(0,MAX), type="l", 
       axes=FALSE, xlab="", ylab="",lty="dashed")
  abline(0,1,col="red")
  par(new=T)
  qqplot(null,obs, xlab=expression(Expected~~-log[10](italic(p))), ylab=expression(Observed~~-log[10](italic(p))),xlim=c(0,MAX),ylim=c(0,MAX), main=title,pch=19)
}