From 801090f1888041e458f1be9b030c501906569df1 Mon Sep 17 00:00:00 2001 From: Ben Laufer Date: Tue, 25 Sep 2018 15:30:59 -0700 Subject: [PATCH] Fix for option parsing --- DM.R | 49 ++++++++++++++++++++++++++++++------------------- DM.R.sh | 7 ++++--- 2 files changed, 34 insertions(+), 22 deletions(-) diff --git a/DM.R b/DM.R index cab7edf..895c88c 100644 --- a/DM.R +++ b/DM.R @@ -7,6 +7,7 @@ # Last update: September 22nd 2018 # Install and Update ------------------------------------------------------ + rm(list=ls()) options(scipen=999) @@ -16,7 +17,8 @@ CRAN <- c("BiocManager", "remotes") new.packages <- CRAN[!(CRAN %in% installed.packages()[,"Package"])] if(length(new.packages)>0){ install.packages(new.packages, repos ="https://cloud.r-project.org") - stopifnot(suppressMessages(sapply(CRAN, require, character.only= TRUE)))} +} +stopifnot(suppressMessages(sapply(CRAN, require, character.only= TRUE))) packages <- c("gplots","RColorBrewer","openxlsx","CMplot","stringr","ggplot2","plyr","devtools","dmrseq", "BiocParallel","annotatr","liftOver","rGREAT","enrichR","ChIPseeker", @@ -29,9 +31,9 @@ new.packages <- packages[!(packages %in% installed.packages()[,"Package"])] if(length(new.packages)>0){ library("BiocManager") BiocManager::install(new.packages, ask = FALSE) - packages <- gsub("vqv/", "", packages) - stopifnot(suppressMessages(sapply(packages, require, character.only= TRUE)))} - +} +packages <- gsub("vqv/", "", packages) +stopifnot(suppressMessages(sapply(packages, require, character.only= TRUE))) suppressWarnings(BiocManager::valid(fix=TRUE, ask=FALSE)) # Global variables -------------------------------------------------------- @@ -42,7 +44,7 @@ option_list <- list( make_option(c("-g", "--genome"), type = "character", default = NULL, help = "Choose a genome (hg38, mm10, rn6, rheMac8) [required]"), make_option(c("-x", "--coverage"), type = "integer", default = 1, - help = "Choose a CpG coverage cutoff [default = 1]"), + help = "Choose a CpG coverage cutoff [default = %default]"), make_option(c("-t", "--testCovariate"), type = "character", default = NULL, help = "Choose a test covariate [required]"), make_option(c("-a", "--adjustCovariate"), type = "character", default = NULL, @@ -50,19 +52,28 @@ option_list <- list( make_option(c("-m", "--matchCovariate"), type = "character", default = NULL, help = "Choose covariate to balance permutations [default = NULL]"), make_option(c("-c", "--cores"), type = "integer", default = 2, - help = "Choose number of cores [default = 2]") + help = "Choose number of cores [default = %default]") ) - opt <- parse_args(OptionParser(option_list=option_list)) +cat("\n[DM.R] Assigning arguments to global variables \t\t", format(Sys.time(), "%d-%m-%Y %X"), "\n") + stopifnot(!is.null(opt$genome)) stopifnot(!is.null(opt$testCovariate)) genome <- as.character(opt$genome) coverage <- as.numeric(opt$coverage) testCovariate <- as.character(opt$testCovariate) -adjustCovariate <- as.character(opt$adjustCovariate) -matchCovariate <- as.character(opt$matchCovariate) +if(!is.null(opt$adjustCovariate)){ + adjustCovariate <- as.character(opt$adjustCovariate) +}else if(is.null(opt$adjustCovariate)){ + adjustCovariate <- opt$adjustCovariate +} +if(!is.null(opt$matchCovariate)){ + matchCovariate <- as.character(opt$matchCovariate) +}else if(is.null(opt$matchCovariate)){ + matchCovariate <- opt$matchCovariate +} cores <- as.numeric(opt$cores) # Setup Annotation Databases ---------------------------------------------- @@ -90,7 +101,7 @@ names <- gsub( "_.*$","", cov) bs<- read.bismark(files = cov, sampleNames = names, - rmZeroCov = FALSE, + rmZeroCov = TRUE, strandCollapse = TRUE, fileType = "cytosineReport", verbose = TRUE, @@ -121,7 +132,7 @@ head(getCoverage(bs.filtered, type="Cov")) cat("\n[DM.R] Saving Rdata \t\t\t\t\t", format(Sys.time(), "%d-%m-%Y %X"), "\n") bismark_env <- ls(all=TRUE) -save(list = bismark_env, file = "bismark.RData") +#save(list = bismark_env, file = "bismark.RData") #load("bismark.RData") # Distribtuion plots ------------------------------------------------------ @@ -129,16 +140,16 @@ save(list = bismark_env, file = "bismark.RData") cat("\n[DM.R] Plotting Empirical Distribution of CpGs \t\t", format(Sys.time(), "%d-%m-%Y %X"), "\n") # Raw -#plotEmpiricalDistribution(bs, testCovariate=testCovariate) -#plotEmpiricalDistribution(bs, testCovariate=testCovariate, type="Cov", bySample=TRUE) +#plotEmpiricalDistribution(bs, testCovariate = testCovariate) +#plotEmpiricalDistribution(bs, testCovariate= testCovariate, type = "Cov", bySample = TRUE) # Filtered pdf("Methylation_Distribution.pdf", height = 7.50, width = 11.50) -plotEmpiricalDistribution(bs.filtered, testCovariate=testCovariate) +plotEmpiricalDistribution(bs.filtered, testCovariate = testCovariate) dev.off() pdf("Methylation_Coverage.pdf", height = 7.50, width = 11.50) -plotEmpiricalDistribution(bs.filtered, testCovariate=testCovariate, type="Cov", bySample=TRUE) +plotEmpiricalDistribution(bs.filtered, testCovariate = testCovariate, type = "Cov", bySample = TRUE) dev.off() # DMRs -------------------------------------------------------------------- @@ -151,7 +162,7 @@ regions <- dmrseq(bs=bs.filtered, cutoff = 0.05, minNumRegion = 5, maxPerms = 10, - testCovariate=testCovariate, + testCovariate = testCovariate, adjustCovariate = adjustCovariate, matchCovariate = matchCovariate) @@ -170,8 +181,8 @@ dev.off() cat("\n[DM.R] Extracing raw differnces for DMRs \t\t", format(Sys.time(), "%d-%m-%Y %X"), "\n") # Get Raw differences -rawDiff <- meanDiff(bs.filtered, dmrs=regions, testCovariate=testCovariate) -sigRawDiff <- meanDiff(bs.filtered, dmrs=sigRegions, testCovariate=testCovariate) +rawDiff <- meanDiff(bs.filtered, dmrs = regions, testCovariate = testCovariate) +sigRawDiff <- meanDiff(bs.filtered, dmrs = sigRegions, testCovariate = testCovariate) # Add Raw differences regions$RawDiff <- rawDiff sigRegions$RawDiff <- sigRawDiff @@ -686,7 +697,7 @@ cat("\n[DM.R] Testing for large blocks (PMDs/HMDs) \t\t", format(Sys.time(), "%d register(MulticoreParam(1)) blocks <- dmrseq(bs=bs.filtered, cutoff = 0.05, - maxPerms = 100, + maxPerms = 10, testCovariate=testCovariate, adjustCovariate = adjustCovariate, matchCovariate = matchCovariate, diff --git a/DM.R.sh b/DM.R.sh index c4c452b..129bd92 100644 --- a/DM.R.sh +++ b/DM.R.sh @@ -13,7 +13,7 @@ # Author: Ben Laufer # Email: blaufer@ucdavis.edu # Last Update Date: 09-20-2018 -# Version: 0.98 +# Version: 0.99.2 # # DMR inference and data visualization for CpG_Me output and bismark cytosine reports # @@ -58,10 +58,11 @@ aklog call="Rscript \ --vanilla \ -/share/lasallelab/programs/CpG_Me/DM.R \ +/share/lasallelab/programs/DM.R \ --genome hg38 \ --coverage 1 \ ---adjustCovariate Diagnosis \ +--testCovariate Diagnosis \ +--adjustCovariate Age \ --cores 2" echo $call