R/plot.clumps.R

##' Plot clumps of markers on Manhattan plot
##' 
##' @description Plot clumps resulting from running the \code{\link[cgmisc]{clump.markers}} function.
##' @author Marcin Kierczak <\email{Marcin.Kierczak@@slu.se}>
##' @param an an object of the \code{\link[GenABEL]{gwaa.data-class}} 
##' @param clumps a result of \code{\link[cgmisc]{clump.markers}} 
##' @param chr chromosome to display
##' @return NULL
##' @examples
##'  \dontrun{example}
##' @keywords plot clumps clumping
##' @seealso \code{\link[cgmisc]{clump.markers}}
##' @export
plot.clumps <- function(an, clumps, chr, region) {
  par(mfrow=c(1,1))
  minCoord <- region[1]
  maxCoord <- region[2]
  region <- which(an@annotation$Chromosome == chr & an@annotation$Position > minCoord & an@annotation$Position < maxCoord)
  pvals <- -log10(an@results$P1df[region])
  coords <- an@annotation$Position[region]
  plot(coords, pvals, pch=19, cex=.5, ann = FALSE, xaxt='n', yaxt='n', bty='n', type='n', ylim = c(-length(clumps)-1,max(pvals)))
  step <- (max(coords) - min(coords)) / 5
  axis(1, at = seq(min(coords), max(coords), by=step), labels=format(seq(min(coords), max(coords), by=step)/1e6, scientific=F, digits=3))
  axis(2, at=seq(0, max(pvals), 1), labels=T)
  #axis(2, at=seq(-1, -1-length(clumps), -1),labels=abs(seq(-1, -1-length(clumps), -1)))
  mtext('Position (Mb)', 1, 3)
  #mtext(expression(clumps~'            '~-log[10](p-value)), 2, 3)
  mtext(expression(-log[10](p-value)), 2, 3)
  points(coords, pvals, pch=19, cex=.5)
  mycols <- colorRampPalette(colors=c("slateblue","grey","red"))
  cols <- mycols(length(clumps))
  for (i in 1:length(clumps)) {
    x <- rep(min(clumps[[i]]$coord), times=2)
    y <- c(-1-i, max(pvals))
    #lines(x,y, col="tomato", lty=2)
    lines(c(min(clumps[[i]]$coord), max(clumps[[i]]$coord)), c(-1-i, -1-i), col=cols[i])
    #points(clumps[[i]]$coord, -log10(clumps[[i]]$pval), col=cols[i], type='l', lwd=2)
    points(clumps[[i]]$coord, -log10(clumps[[i]]$pval), col=cols[i], cex=.8, pch=19)
    #abline(v=(min(clumps[[i]]$coord)+max(clumps[[i]]$coord))/2, col="tomato", lty=2)
    points(clumps[[i]]$coord, rep(-1-i, times=length(clumps[[i]]$coord)), col=cols[i], pch=19, cex=0.8)
  }
}

Try the cgmisc package in your browser

Any scripts or data that you put into this service are public.

cgmisc documentation built on May 2, 2019, 5:50 p.m.