R/mhtplot2.R

In gap: Genetic Analysis Package

#' Manhattan plot with annotations
#'
#' @param data a data frame with three columns representing chromosome, position and p values.
#' @param control A list produced by [mht.control()] controlling plot behaviour.
#' @param hcontrol A list produced by [hmht.control()] defining highlighted markers or regions and labels.
#' @param ... Additional graphical arguments passed to [graphics::plot()]
#'   (e.g. `pch`, `bg`, `xlab`, `ylab`, `ylim`).
#'   If `mar` is supplied, the default margins are not modified.
#'
#' @details
#' 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. The example uses only chromosomes 14 and 20 of \insertCite{denhoed13}{gap}.
#'
#' @export
#' @return shown on or saved to the appropriate device.
#' @examples
#' \dontrun{
#' 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)
#' hops <- hmht.control(data=subset(mdata,!is.na(gene)),boxed=TRUE)
#' 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()
#'
#' # To exemplify the use of chr, pos and p without gene annotation
#' # in response to query from Vallejo, Roger <Roger.Vallejo@ARS.USDA.GOV>
#' opar <- par()
#' par(cex=0.4)
#' ops <- mht.control(colors=rep(c("lightgray","lightblue"),11),yline=2.5,xline=2)
#' mhtplot2(data.frame(mhtdata[,c("chr","pos","p")],gene=NA,color=NA),ops,xlab="",ylab="")
#' axis(2,at=1:16)
#' title("data in mhtplot used by mhtplot2")
#' par(opar)
#' }
#' @references
#' \insertRef{denhoed13}{gap}
#'
#' @author Jing Hua Zhao
#'
mhtplot2 <- function (data, control = mht.control(), hcontrol = hmht.control(), ...)
{
    for(p in c("grid")) {
       if (length(grep(paste("^package:", p, "$", sep=""), search())) == 0) {
          if (!requireNamespace(p, quietly = TRUE))
          warning(paste("mhtplot2 needs package `", p, "' to be fully functional; please install", sep=""))
       }
    }
    nonmiss <- !apply(is.na(data[,1:3]), 1, any)
    data2 <- data[nonmiss, ]
    chr <- data2[, 1]
    pos <- newpos <- data2[, 2]
    p <- data2[, 3]
    rc <- data2[, 4]
    tablechr <- table(chr)
    allchr <- as.vector(tablechr)
    n.chr <- length(allchr)
    type <- control$type
    usepos <- control$usepos
    logscale <- control$logscale
    base <- control$base
    cutoffs <- control$cutoffs
    colors <- control$colors
    labels <- control$labels
    srt <- control$srt
    gap <- control$gap
    pcex <- control$cex
    yline <- control$yline
    xline <- control$xline
    colorlist <- colors()
    if (is.null(colors)) 
        colors <- sample(colorlist, n.chr)
    tablechr <- unique(chr)
    if (is.null(labels)) 
        labels <- tablechr
    if (is.null(gap)) 
        gap <- 0
    if (!is.null(hcontrol$data)) {
        hdata <- hcontrol$data
        hdata2 <- hdata[!apply(is.na(hdata), 1, any), ]
        hchr <- hdata2[, 1]
        hpos <- hnewpos <- hdata2[, 2]
        hp <- hdata2[, 3]
        hname <- hdata2[, 4]
        hcol <- hdata2[, 5]
        namecol <- table(hname,hcol)
        namecol[namecol!=0] <- 1
        hcolors <- colnames(namecol)
        gmat <- namecol%*%(1:length(hcolors))
        hyoffs <- hcontrol$yoffs
        hboxed <- hcontrol$boxed
        hcex <- hcontrol$cex
        htablechr <- unique(hchr)
        hn.chr <- length(htablechr)
        hlabels <- unique(hname)
        htablechrname <- unique(data.frame(hchr, hname))
        if (is.null(hcolors)) 
            hcolors <- rep("red", length(hlabels))
        else hcolors <- hcontrol$colors
    }
    CMindex <- cumsum(allchr)
    for (i in 1:n.chr) {
        u <- CMindex[i]
        l <- CMindex[i] - allchr[i] + 1
        chr <- l:u
        if (usepos) 
            d <- diff(pos[chr])
        else d <- rep(1, allchr[i] - 1)
        newpos[chr] <- c(gap, d)
    }
    CM <- cumsum(as.numeric(newpos))
    args <- list(...)
    if ("ylim" %in% names(args)) 
        dp <- seq(args$ylim[1], args$ylim[2], length = sum(allchr))
    else dp <- seq(min(p), max(p), length = sum(allchr))
    if (logscale) 
        y <- -log(dp, base)
    else y <- dp
    y1 <- min(y)
    par(xaxt = "n", yaxt = "n")
    xy <- xy.coords(CM, y)
    plot(xy$x, xy$y, type = "n", ann = FALSE, axes = FALSE, ...)
    axis(1)
    axis(2)
    par(xaxt = "s", yaxt = "s")
    for (i in 1:n.chr) {
        u <- CMindex[i]
        l <- CMindex[i] - allchr[i] + 1
        chr <- l:u
        cat("Plotting points ", l, "-", u, "\n")
        if (logscale)
            y <- -log(p[chr], base)
        else y <- p[chr]
        col.chr <- colors()[rc[chr]]
        col.chr[is.na(rc[chr])]<-colors[i]
        if (type == "l")
            lines(CM[chr], y, col = col.chr, cex = pcex, ...)
        else points(CM[chr], y, col = col.chr, cex = pcex, ...)
        text(ifelse(i == 1, CM[1]+(CM[u]-CM[l])/2, CM[l]+(CM[u]-CM[l])/2), y1, pos = 1, offset = 1.5, 
             labels[i], srt = srt, ...)
    }
    for (i in 1:n.chr) {
        u <- CMindex[i]
        l <- CMindex[i] - allchr[i] + 1
        chr <- l:u
        if (logscale) 
            y <- -log(p[chr], base)
        else y <- p[chr]
        col.chr <- colors()[rc[chr]]
        if (type =="l"&all(is.na(rc[chr])))
          lines(CM[chr], y, col = col.chr, cex = pcex, ...)
        else points(CM[chr], y, col = col.chr, cex = pcex, ...)
        if (!is.null(hcontrol$data)) {
            chrs <- htablechrname[tablechr[i] == htablechrname[, 
                1], ]
            if (dim(chrs)[1] > 0) {
                hchrs <- as.character(chrs[, 2])
                for (k in 1:length(hchrs)) {
                  hregion <- hpos[hchr == chrs[k, 1] & hname == 
                    hchrs[k]]
                  hl <- chr[pos[chr] == hregion[1]]
                  hu <- chr[pos[chr] == hregion[length(hregion)]]
                  cat("  ... highlighting", hl, "-", hu, hchrs[k], 
                    "\n")
                  l1 <- hl - l + 1
                  l2 <- hu - l + 1
                  col.index <- as.integer(colnames(namecol)[gmat[rownames(gmat)==hchrs[k]]])
                  col.label <- colors()[col.index]
                  idx <- chr
                  hx <- CM[idx][l1]
                  hy <- max(y[l1:l2]) + hyoffs
                  lab <- hchrs[k]
                  if (hboxed) {
                     usr <- par("usr")
                     pin <- par("pin")
                     xinch <- diff(usr[1:2]) / pin[1]
                     yinch <- diff(usr[3:4]) / pin[2]
                     w <- strwidth(lab, units="inches") * xinch * 1.3
                     h <- strheight(lab, units="inches") * yinch * 1.6
                     rect(hx - w/2, hy - h/2, hx + w/2, hy + h/2,
                          col="white", border="black")
                     text(hx, hy, lab, cex=hcex, font=3)
                  } else {
                     text(hx, hy, lab, cex=hcex, font=3)
                  }
                }
            }
        }
    }
    if (!is.null(cutoffs)) 
        abline(h = cutoffs)
    if ("ylab" %in% names(args)) 
        mtext(args$ylab, 2, line = yline, las = 0, cex=0.5, ...)
    else mtext(ifelse(logscale, paste("-log", base, "(Observed value)", 
        sep = ""), "Observed value"), 2, line = yline, las = 0, cex=0.5, ...)
    if ("xlab" %in% names(args)) 
        xlabel <- args$xlab
    else xlabel <- ifelse(is.null(names(chr)), "Chromosome", 
        names(chr))
    mtext(xlabel, 1, line = xline, las = 0, cex=0.5, ...)
}