|
| 1 | +#!/usr/bin/env Rscript |
| 2 | +#' @(#)GWAS_manhattanplots 2017-10-05 A.Douglas and A.J.Travis |
| 3 | + |
| 4 | +suppressMessages(library(getopt)) |
| 5 | +suppressMessages(library(optparse)) |
| 6 | + |
| 7 | +#' debug in RStudio |
| 8 | + |
| 9 | +# command-line arguments |
| 10 | +opt1 <- make_option(c("-b", "--bh"), action = "store", default = 0) |
| 11 | +opt2 <- make_option(c("-c", "--chromosome"), action = "store", default = "1:23") |
| 12 | +opt3 <- make_option(c("-i", "--input"), action = "store", default = NULL) |
| 13 | +opt4 <- make_option(c("-o", "--output"), action = "store", default = "manhattan") |
| 14 | +opt5 <- make_option(c("-q", "--qqplot"), action = "store_true", default = FALSE) |
| 15 | +opt6 <- make_option(c("-v", "--verbose"), action = "store_true", default = FALSE) |
| 16 | +opt7 <- make_option(c("-y", "--ymax"), action = "store", default = "max") |
| 17 | + |
| 18 | +#' debug in RStudio |
| 19 | +#' args <- c('-b', '0.1', '-q', '-y', '10', '-i', 'files.txt') |
| 20 | +#' opt <- parse_args(OptionParser(option_list = list(opt1, opt2, opt3, opt4, opt5, opt6, |
| 21 | +#' opt7)), args = args, positional_arguments = TRUE) |
| 22 | + |
| 23 | +opt <- parse_args(OptionParser(option_list = list(opt1, opt2, opt3, opt4, opt5, opt6, |
| 24 | + opt7)), positional_arguments = TRUE) |
| 25 | + |
| 26 | +# limit number of chromosomes |
| 27 | +chr.name <- parse(text = opt$options$chromosome) |
| 28 | +chr.limit <- eval(parse(text = chr.name)) |
| 29 | +chr.class <- class(chr.limit) |
| 30 | +if (chr.class == "numeric") { |
| 31 | + chr.limit <- as.vector(chr.limit) |
| 32 | +} else { |
| 33 | + if (chr.class != "integer") { |
| 34 | + stop("list of chromosomes required e.g. 1:23") |
| 35 | + } |
| 36 | +} |
| 37 | + |
| 38 | +# PNG output file |
| 39 | +output.file <- paste(sep = "", opt$options$output, ".png") |
| 40 | + |
| 41 | +# check usage |
| 42 | +if (is.null(opt$options$input)) { |
| 43 | + nplots <- length(opt$args) |
| 44 | +} else { |
| 45 | + title <- read.table(opt$options$input, col.names = c("filename", "title"), colClasses = "character", |
| 46 | + sep = "\t") |
| 47 | + nplots <- length(title$filename) |
| 48 | +} |
| 49 | + |
| 50 | +if (nplots == 0) { |
| 51 | + stop("one or more '.manh' files required") |
| 52 | +} |
| 53 | + |
| 54 | +# script location |
| 55 | +script <- get_Rscript_filename() |
| 56 | +basename <- dirname(script) |
| 57 | + |
| 58 | +# source library |
| 59 | +lib.name <- paste(sep = "/", basename, "GWAS_manhattanplots.R") |
| 60 | +source(lib.name) |
| 61 | + |
| 62 | +# Required for Benjamini-Hochberg correction |
| 63 | +suppressMessages(library(multtest)) # include check if installed |
| 64 | + |
| 65 | +# sides |
| 66 | +BOTTOM <- 1 |
| 67 | +LEFT <- 2 |
| 68 | +TOP <- 3 |
| 69 | +RIGHT <- 4 |
| 70 | + |
| 71 | +# outer margins (inches) |
| 72 | +OM_BOTTOM <- 0 |
| 73 | +OM_LEFT <- 0.4 |
| 74 | +OM_TOP <- 0.2 |
| 75 | +OM_RIGHT <- 0.4 |
| 76 | + |
| 77 | +# inner margins (inches) |
| 78 | +IM_BOTTOM <- 0 |
| 79 | +IM_LEFT <- 0.4 |
| 80 | +IM_TOP <- 0.2 |
| 81 | +IM_RIGHT <- 0 |
| 82 | +IM_LAST <- 0.8 |
| 83 | + |
| 84 | +# overall plot height (inches) |
| 85 | +PLOTH <- 1.4 |
| 86 | + |
| 87 | +# Manhattan plot aspect ratio (width / height) |
| 88 | +MPLOTR <- 5 |
| 89 | + |
| 90 | +# Manhattan plot dimensions (inches) |
| 91 | +MPLOTH <- IM_TOP + PLOTH + IM_BOTTOM |
| 92 | +MPLOTW <- IM_LEFT + PLOTH * MPLOTR + IM_RIGHT |
| 93 | + |
| 94 | +# QQ plot dimensions (inches) |
| 95 | +QPLOTH <- MPLOTH |
| 96 | +QPLOTW <- IM_LEFT + PLOTH + IM_RIGHT |
| 97 | + |
| 98 | +#' markers <- read.table('marker_pos.txt', header = TRUE) |
| 99 | +#' blocks <- read.table('qtl_blocks.txt', header = TRUE) |
| 100 | + |
| 101 | +# page layout |
| 102 | +page.height <- OM_TOP + (MPLOTH * nplots) + IM_LAST + OM_BOTTOM |
| 103 | +plot.height <- c(rep.int(1, nplots - 1), (MPLOTH + IM_LAST)/MPLOTH) |
| 104 | +if (opt$options$qqplot) { |
| 105 | + page.width <- OM_LEFT + MPLOTW + QPLOTW + OM_RIGHT |
| 106 | + png(output.file, width = page.width, height = page.height, unit = "in", res = 300) |
| 107 | + layout(mat = matrix(c(1:(nplots * 2)), nplots, 2, byrow = TRUE), widths = c(MPLOTW, |
| 108 | + QPLOTW), heights = plot.height) |
| 109 | +} else { |
| 110 | + page.width <- OM_LEFT + MPLOTW + OM_RIGHT |
| 111 | + png(output.file, width = page.width, height = page.height, unit = "in", res = 300) |
| 112 | + layout(mat = matrix(c(1:(nplots)), nplots, 1, byrow = TRUE), heights = plot.height) |
| 113 | +} |
| 114 | + |
| 115 | +# outer margins: c(bottom, left, top, right) |
| 116 | +par(omi = c(OM_BOTTOM, OM_LEFT, OM_TOP, OM_RIGHT), lend = "square") |
| 117 | + |
| 118 | +# stack plots |
| 119 | +for (i in 1:nplots) { |
| 120 | + |
| 121 | + # only show axis on last stacked plot |
| 122 | + stack <- i < nplots |
| 123 | + if (stack) { |
| 124 | + par(mai = c(IM_BOTTOM, IM_LEFT, IM_TOP, IM_RIGHT)) |
| 125 | + } else { |
| 126 | + par(mai = c(IM_BOTTOM + IM_LAST, IM_LEFT, IM_TOP, IM_RIGHT)) |
| 127 | + } |
| 128 | + |
| 129 | + # read SNPs from file |
| 130 | + if (is.null(opt$options$input)) { |
| 131 | + snp <- read.table(opt$args[i], header = TRUE) |
| 132 | + } else { |
| 133 | + snp <- read.table(title$filename[i], header = TRUE) |
| 134 | + } |
| 135 | + |
| 136 | + # exclude P <= 0 (error flag?) |
| 137 | + snp <- snp[snp$P > 0, ] |
| 138 | + |
| 139 | + # calculate BH (Benjamini-Hochberg) adjusted P values |
| 140 | + mt <- mt.rawp2adjp(snp$P, proc = "BH") |
| 141 | + ind <- mt$index |
| 142 | + snp <- cbind(snp, BH = mt$adjp[order(ind), 2]) |
| 143 | + |
| 144 | + # Manhattan plot |
| 145 | + manhattan(snp, markers = NULL, blocks = NULL, limitchromosomes = chr.limit, ymax = opt$options$ymax, |
| 146 | + suggestiveline = 4, genomewideline = FALSE, highlight = opt$options$bh, stack = stack) |
| 147 | + |
| 148 | + # sub-plot title (left-justified) |
| 149 | + if (!is.null(opt$options$input)) { |
| 150 | + text <- title$title[i] |
| 151 | + } else { |
| 152 | + text <- sub(".manh", "", basename(opt$args[i])) |
| 153 | + } |
| 154 | + mtext(text, side = TOP, line = -0.8, adj = 0.02, cex = 0.8) |
| 155 | + |
| 156 | + # QQ plot |
| 157 | + if (opt$options$qqplot) { |
| 158 | + qq(snp, limitchromosomes = chr.limit, ymax = opt$options$ymax, highlight = opt$options$bh, |
| 159 | + stack = stack) |
| 160 | + } |
| 161 | +} |
| 162 | + |
| 163 | +# add axis labels in outer margins |
| 164 | +mtext(expression(" " - log[10](italic(p))), outer = TRUE, side = LEFT, line = 0) |
| 165 | + |
| 166 | +# check device was closed |
| 167 | +if (dev.off() != 1) { |
| 168 | + stop("error: Graphics device still open") |
| 169 | +} |
0 commit comments