Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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)
- }
Add Comment
Please, Sign In to add comment