Mercurial > repos > iuc > gwastools_manhattan_plot
diff manhattan.R @ 0:c4bf5e913c2e draft default tip
"planemo upload commit b1883bac95e73fc6ffe2a36db3115ad5e5a1eba4"
| author | iuc |
|---|---|
| date | Fri, 11 Oct 2019 17:31:09 -0400 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/manhattan.R Fri Oct 11 17:31:09 2019 -0400 @@ -0,0 +1,60 @@ +options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) + +# we need that to not crash galaxy with an UTF8 error on German LC settings. +loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") + +suppressPackageStartupMessages({ + library(GWASTools) + library(optparse) +}) +option_list <- list( + make_option(c("-f", "--file"), type="character", help="RInput GWAS file"), + make_option("--pval", type="integer", help="Pvalue column"), + make_option("--chromosome", type="integer", help="Chromosome column"), + make_option("--ymin", type="double", help="Min y value"), + make_option("--ymax", type="double", help="Max y value"), + make_option("--trunc", help="Show truncation lines", action="store_true"), + make_option("--sig", type="double", help="Significance level for lines"), + make_option("--thin", type="double", help="Thinning value", action="store_true", dest="thin"), + make_option("--ppb", type="integer", help="Points per bin, if thinning value is specified", action="store_true")) +args <- parse_args(OptionParser(option_list=option_list)) +file <- args$file +data <- read.table(args$file, header=TRUE) +pval <- data[,args$pval] +chromosome <- data[,args$chromosome] +if(!is.null(args$ymin) & !is.null(args$ymax)){ + ylimit <- c(args$ymin,args$ymax) +}else if(xor(!is.null(args$ymin), !is.null(args$ymax))){ + print("If specifying range, both ymin and ymax must be set") + ylimit <- NULL +}else{ + ylimit <- NULL +} +if(is.null(args$trunc)){ + trunc <- FALSE +}else{ + trunc <- TRUE +} +if(!is.null(args$sig)){ + sig <- args$sig +}else{ + sig <- NULL +} +if(!is.null(args$thin)){ + thin = args$thin + if(thin == 0){ + thin = NULL + } +}else{ + thin <- FALSE +} +if(!is.null(thin) & !is.null(args$ppb)){ + ppb = args$ppb +} +pdf("manhattan.pdf") +if(isFALSE(thin)){ + manhattanPlot(pval, chromosome, ylim=ylimit, trunc.lines=trunc, signif=sig) +}else{ + manhattanPlot(pval, chromosome, ylim=ylimit, trunc.lines=trunc, signif=sig, thinThreshold = thin, pointsPerBin = ppb) +} +invisible(dev.off())
