manhattanPlus.R

manhattanPlus <- function (x, chr = "CHR", bp = "BP", p = "P", snp = "SNP",
                           col = c("gray10", "gray60"), chrlabs = NULL,
                           suggestiveline = -log10(1e-05), 
                           genomewideline = -log10(5e-08), highlight = NULL,
                           logp = TRUE, annotatePval = NULL, annotateTop = TRUE,
                           otherColors = rev(viridis::viridis(3)[1:2]),
                           otherDataColName = FALSE, otherCatName = FALSE,
                           ...) 
{
  CHR = BP = P = index = NULL
  if (!(chr %in% names(x))) 
    stop(paste("Column", chr, "not found!"))
  if (!(bp %in% names(x))) 
    stop(paste("Column", bp, "not found!"))
  if (!(p %in% names(x))) 
    stop(paste("Column", p, "not found!"))
  if (! identical(otherDataColName, FALSE) & !(otherDataColName %in% names(x)))
    stop(paste("Column", otherDataColName, "not found! If you don't want to use it, set it to FALSE, as default."))
  if (! identical(otherCatName, FALSE) & !(otherCatName %in% names(x)))
    stop(paste("Column", otherCatName, "not found! If you don't want to use it, set it to FALSE, as default."))
  if (!(snp %in% names(x))) 
    warning(paste("No SNP column found. OK unless you're trying to highlight."))
  if (!is.numeric(x[[chr]])) 
    stop(paste(chr, "column should be numeric. Do you have 'X', 'Y', 'MT', etc? If so change to numbers and try again."))
  if (!is.numeric(x[[bp]])) 
    stop(paste(bp, "column should be numeric."))
  if (!is.numeric(x[[p]])) 
    stop(paste(p, "column should be numeric."))
  d = data.frame(CHR = x[[chr]], BP = x[[bp]], P = x[[p]])
  if (! identical(otherDataColName, FALSE) & (otherDataColName %in% names(x)))
    d$score = x[[otherDataColName]]
  if (! identical(otherCatName, FALSE) & (otherCatName %in% names(x)))
    d$cat = x[[otherCatName]]
  if (!is.null(x[[snp]])) 
    d = transform(d, SNP = x[[snp]])
  d <- subset(d, (is.numeric(CHR) & is.numeric(BP) & is.numeric(P)))
  d <- d[order(d$CHR, d$BP), ]
  if (logp) {
    d$logp <- -log10(d$P)
  } else {
    d$logp <- d$P
  }
  d$pos = NA
  d$index = NA
  ind = 0
  for (i in unique(d$CHR)) {
    ind = ind + 1
    d[d$CHR == i, ]$index = ind
  }
  nchr = length(unique(d$CHR))
  if (nchr == 1) {
    d$pos = d$BP
    ticks = floor(length(d$pos))/2 + 1
    xlabel = paste("Chromosome", unique(d$CHR), "position")
    labs = ticks
  } else {
    lastbase = 0
    ticks = NULL
    for (i in unique(d$index)) {
      if (i == 1) {
        d[d$index == i, ]$pos = d[d$index == i, ]$BP
      }
      else {
        lastbase = lastbase + tail(subset(d, index == 
                                            i - 1)$BP, 1)
        d[d$index == i, ]$pos = d[d$index == i, ]$BP + 
          lastbase
      }
      ticks = c(ticks, (min(d[d$index == i, ]$pos) + max(d[d$index == 
                                                             i, ]$pos))/2 + 1)
    }
    xlabel = "Chromosome"
    labs <- unique(d$CHR)
  }
  xmax = ceiling(max(d$pos) * 1.03)
  xmin = floor(max(d$pos) * -0.03)
  dotargs <- list(...)
  if (! identical(otherDataColName, FALSE)){
    e <- d[! is.na(d$score),]
    def_args <- list(xaxt = "n", yaxt = "n", bty = "n", xaxs = "i", yaxs = "i", 
                     las = 1, xlim = c(xmin, xmax),
                     ylim = c(0, ceiling(max(e$score))), xlab = "",
                     ylab = "")
    do.call("plot", c(NA, dotargs, def_args[!names(def_args) %in% 
                                              names(dotargs)]))
    rect(xleft = e$pos-1, ybottom = 0, xright = e$pos+1,ytop = e$score,
         col = otherColors[e$cat], border = otherColors[e$cat], ...)
    axis(4, col = otherColors[1], col.axis = otherColors[1], las = 2)
    mtext("Importance score", side = 4, line = 2, cex = 0.7)
    legend("topright", legend = unique(e$cat), fill = rev(otherColors),
           bg = "transparent")
  }
  def_args <- list(xaxt = "n", bty = "n", xaxs = "i", yaxs = "i", 
                   las = 1, pch = 20, xlim = c(xmin, xmax),
                   ylim = c(0, ceiling(max(d$logp))), xlab = xlabel,
                   ylab = expression(-log[10](italic(p))))
  par(new=TRUE)
  do.call("plot", c(NA, dotargs, def_args[!names(def_args) %in% 
                                            names(dotargs)]))
  if (!is.null(chrlabs)) {
    if (is.character(chrlabs)) {
      if (length(chrlabs) == length(labs)) {
        labs <- chrlabs
      }
      else {
        warning("You're trying to specify chromosome labels but the number of labels != number of chromosomes.")
      }
    }
    else {
      warning("If you're trying to specify chromosome labels, chrlabs must be a character vector")
    }
  }
  if (nchr == 1) {
    axis(1, ...)
  } else {
    axis(1, at = ticks, labels = labs, ...)
  }
  col = rep(col, max(d$CHR))
  if (nchr == 1) {
    with(d, points(pos, logp, pch = 20, col = col[1], ...))
  } else {
    icol = 1
    for (i in unique(d$index)) {
      with(d[d$index == unique(d$index)[i], ], points(pos, 
                                                      logp, col = col[icol], pch = 20, ...))
      icol = icol + 1
    }
  }
  if (suggestiveline) 
    abline(h = suggestiveline, col = "blue")
  if (genomewideline) 
    abline(h = genomewideline, col = "red")
  if (!is.null(highlight)) {
    if (any(!(highlight %in% d$SNP))) 
      warning("You're trying to highlight SNPs that don't exist in your results.")
    d.highlight = d[which(d$SNP %in% highlight), ]
    with(d.highlight, points(pos, logp, col = "green3", pch = 20, 
                             ...))
  }
  if (!is.null(annotatePval)) {
    topHits = subset(d, P <= annotatePval)
    par(xpd = TRUE)
    if (annotateTop == FALSE) {
      with(subset(d, P <= annotatePval), textxy(pos, -log10(P), 
                                                offset = 0.625, labs = topHits$SNP, cex = 0.45), 
           ...)
    } else {
      topHits <- topHits[order(topHits$P), ]
      topSNPs <- NULL
      for (i in unique(topHits$CHR)) {
        chrSNPs <- topHits[topHits$CHR == i, ]
        topSNPs <- rbind(topSNPs, chrSNPs[1, ])
      }
      textxy(topSNPs$pos, -log10(topSNPs$P), offset = 0.625, 
             labs = topSNPs$SNP, cex = 0.5, ...)
    }
  }
  par(xpd = FALSE)
}