comparison GO_MWU.R @ 2:5acf9dfdfa27 draft default tip

planemo upload commit 66a856bcce69986d9a6f1a39820dd9b3f4f6b0db

1:f7287f82602f

loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")

library("getopt")
library("tools")
options(stringAsFactors = FALSE, useFancyQuotes = FALSE)
args <- commandArgs(trailingOnly = TRUE)

# get options, using the spec as defined by the enclosed list.
# we read the options from the default: commandArgs(TRUE).
spec <- matrix(c(
    "help", "h", 0, "logical",
    "scriptdir", "s", 1, "character",
    "input", "i", 1, "character",
    "goAnnotations", "a", 1, "character",
    "goDatabase", "g", 1, "character",
    "goDivision", "d", 1, "character",
    "largest", "o", 1, "numeric",
    "smallest", "m", 1, "numeric",
    "clusterheight", "c", 1, "numeric",
    "pcut", "p", 1, "numeric",
    "hcut", "t", 1, "numeric",
    "version", "v", 0, "character"
), byrow = TRUE, ncol = 4)
opt <- getopt(spec)

# if help was asked for print a friendly message

    cat(getopt(spec, usage = TRUE))
    q(status = 1)
}

if (!is.null(opt$version)) {
    cat("0.1.0\n")
    q(status = 1)
}

# enforce the following required arguments
if (is.null(opt$scriptdir)) {
    cat("'scriptdir' is required\n")
    q(status = 1)
}
if (is.null(opt$input)) {
    cat("'input' is required\n")
    q(status = 1)
}
if (is.null(opt$goDatabase)) {

if (is.null(opt$goDivision)) {
    cat("'goDivision' is required\n")
    q(status = 1)
}
if (is.null(opt$clusterheight)) {
    cat("'clusterheight' is required\n")
    q(status = 1)
}
if (is.null(opt$largest)) {
    cat("'largest' is required\n")
    q(status = 1)
}
if (is.null(opt$smallest)) {
    cat("'smallest' is required\n")
    q(status = 1)
}
if (is.null(opt$pcut)) {
    cat("'pcut' is required\n")
    q(status = 1)
}
if (is.null(opt$hcut)) {
    cat("'hcut' is required\n")
    q(status = 1)
}

source_local <- function(fname) {
    # argv <- commandArgs(trailingOnly = FALSE)
    # base_dir <- dirname(substring(argv[grep("--file = ", argv)], 8))
    base_dir <- opt$scriptdir
    source(paste(base_dir, fname, sep = "/"))
}

source_local("gomwu.functions.R")

# It might take a few minutes for MF and BP. Do not rerun it if you just want
# to replot the data with different cutoffs, go straight to gomwuPlot. If you
# change any of the numeric values below, delete the files that were generated
# in previos runs first.

gomwuStats(opt$input, opt$goDatabase, opt$goAnnotations, opt$goDivision, opt$scriptdir,
    # replace with full path to perl executable if it is not in your system's PATH already
    perlPath = "perl",
    # a GO category will not be considered if it contains more than this fraction of the total number of genes
    largest = opt$largest,
    smallest = opt$smallest, # a GO category should contain at least this many genes to be considered.

# ----------- Plotting results

# change this to a pdf output
pdf(paste0(dir,"/","Rplots.pdf"))

results <- gomwuPlot(opt$input, opt$goAnnotations, opt$goDivision,
    # genes with the measure value exceeding this will be counted as "good genes".
    # This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing
    # Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes").
    absValue = -log(0.05, 10),
    # 	absValue = 1, # un-remark this if you are using log2-fold changes
    # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue.
    level1 = 0.1,
    level2 = 0.05, # FDR cutoff to print in regular (not italic) font.
    level3 = 0.01, # FDR cutoff to print in large bold font.
    # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text)
    # for better "word cloud" effect
    txtsize = 1.2,
    treeHeight = 0.5, # height of the hierarchical clustering tree
    # 	colors = c("dodgerblue2","firebrick1","skyblue2","lightcoral")
    # these are default colors, uncomment and change if needed
)
# manually rescale the plot so the tree matches the text
# if there are too many categories displayed, try make it more stringent with level1 = 0.05,level2=0.01,level3=0.001.

# text representation of results, with actual adjusted p-values
write.table(results[[1]], paste0(dir, "/", "results.tsv"), sep = "\t")

    if (bestrr$pval[best] <= pcut) {
        annots <- c(annots, sub("\\d+\\/\\d+ ", "", row.names(bestrr)[best]))
    }
}

mwus <- read.table(paste0(dir,"/",paste("MWU", opt$goDivision, name, sep = "_"), ".", ext), header = T)
best_GOs <- mwus[mwus$name %in% annots, ]
write.table(best_GOs, paste0(dir, "/","best_go.tsv"), sep = "\t")
dev.off()

2:5acf9dfdfa27

loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")

library("getopt")
library("tools")
options(stringAsFactors = FALSE, useFancyQuotes = FALSE)

getScriptPath <- function() {
    cmd.args <- commandArgs()
    m <- regexpr("(?<=^--file=).+", cmd.args, perl = TRUE)
    script.dir <- dirname(regmatches(cmd.args, m))
    if (length(script.dir) == 0) stop("can't determine script dir: please call the script with Rscript")
    if (length(script.dir) > 1) stop("can't determine script dir: more than one '--file' argument detected")
    return(script.dir)
}

args <- commandArgs(trailingOnly = TRUE)

# get options, using the spec as defined by the enclosed list.
# we read the options from the default: commandArgs(TRUE).
spec <- matrix(c(
    "help", "h", 0, "logical",
    "input", "i", 1, "character",
    "goAnnotations", "a", 1, "character",
    "goDatabase", "g", 1, "character",
    "goDivision", "d", 1, "character",
    "largest", "o", 1, "numeric",
    "smallest", "m", 1, "numeric",
    "absValue", "k", 1, "numeric",
    "clusterheight", "c", 1, "numeric",
    "textsize", "e", 1, "numeric",
    "pcut", "p", 1, "numeric",
    "hcut", "t", 1, "numeric",
    "l1", "1", 1, "numeric",
    "l2", "2", 1, "numeric",
    "l3", "3", 1, "numeric",
    "version", "v", 0, "character"
), byrow = TRUE, ncol = 4)
opt <- getopt(spec)

# if help was asked for print a friendly message

    cat(getopt(spec, usage = TRUE))
    q(status = 1)
}

if (!is.null(opt$version)) {
    cat("0.3.0\n")
    q(status = 1)
}

scriptdir <- getScriptPath()

### if running in Rstudio ###
# Uncomment this block, change the values and run from here.
# opt = list()
# opt$input = "heats.csv"
# opt$goAnnotations = "amil_defog_iso2go.tab"
# opt$goDatabase = "go.obo"
# opt$goDivision = "BP"

### Optionally uncomment the options below that you want to change and change the default value.
# opt$clusterheight <- 0.25
# opt$textsize <- 1.2
# opt$largest <- 0.1
# opt$smallest <- 5
# opt$absValue <- -log(0.05,10)
# opt$l1 <- 0.1
# opt$l2 <- 0.05
# opt$l3 <- 0.01
# opt$pcut <- 1e-2
# opt$hcut <- 0.9


# enforce the following required arguments
if (is.null(opt$input)) {
    cat("'input' is required\n")
    q(status = 1)
}
if (is.null(opt$goDatabase)) {

if (is.null(opt$goDivision)) {
    cat("'goDivision' is required\n")
    q(status = 1)
}
if (is.null(opt$clusterheight)) {
    opt$clusterheight <- 0.25
}
if (is.null(opt$textsize)) {
    opt$textsize <- 1.2
}
if (is.null(opt$largest)) {
    opt$largest <- 0.1
}
if (is.null(opt$smallest)) {
    opt$smallest <- 5
}
if (is.null(opt$absValue)) {
    opt$absValue <- -log(0.05, 10)
}
if (is.null(opt$l1)) {
    opt$l1 <- 0.1
}
if (is.null(opt$l2)) {
    opt$l2 <- 0.05
}
if (is.null(opt$l3)) {
    opt$l3 <- 0.01
}
if (is.null(opt$pcut)) {
    opt$pcut <- 1e-2
}
if (is.null(opt$hcut)) {
    opt$hcut <- 0.9
}

source_local <- function(fname) {
    # argv <- commandArgs(trailingOnly = FALSE)
    # base_dir <- dirname(substring(argv[grep("--file = ", argv)], 8))
    base_dir <- scriptdir
    source(paste(base_dir, fname, sep = "/"))
}

source_local("gomwu.functions.R")

# It might take a few minutes for MF and BP. Do not rerun it if you just want
# to replot the data with different cutoffs, go straight to gomwuPlot. If you
# change any of the numeric values below, delete the files that were generated
# in previos runs first.

gomwuStats(opt$input, opt$goDatabase, opt$goAnnotations, opt$goDivision, scriptdir,
    # replace with full path to perl executable if it is not in your system's PATH already
    perlPath = "perl",
    # a GO category will not be considered if it contains more than this fraction of the total number of genes
    largest = opt$largest,
    smallest = opt$smallest, # a GO category should contain at least this many genes to be considered.

# ----------- Plotting results

# change this to a pdf output
pdf(paste0(dir, "/", "Rplots.pdf"), width = 7, height = 7)
# png(paste0(dir,"/","Rplots.png"),res=100)

results <- gomwuPlot(opt$input, opt$goAnnotations, opt$goDivision,
    # genes with the measure value exceeding this will be counted as "good genes".
    # This setting is for signed log-pvalues. Specify absValue=0.001 if you are doing
    # Fisher's exact test for standard GO enrichment or analyzing a WGCNA module (all non-zero genes = "good genes").
    absValue = opt$absValue,
    # 	absValue = 1, # un-remark this if you are using log2-fold changes
    # FDR threshold for plotting. Specify level1=1 to plot all GO categories containing genes exceeding the absValue.
    level1 = opt$l1,
    level2 = opt$l2, # FDR cutoff to print in regular (not italic) font.
    level3 = opt$l3, # FDR cutoff to print in large bold font.
    # decrease to fit more on one page, or increase (after rescaling the plot so the tree fits the text)
    # for better "word cloud" effect
    txtsize = opt$textsize,
    treeHeight = 0.5, # height of the hierarchical clustering tree
    # 	colors = c("dodgerblue2","firebrick1","skyblue2","lightcoral")
    # these are default colors, uncomment and change if needed
)
dev.off()
# manually rescale the plot so the tree matches the text
# if there are too many categories displayed, try make it more stringent with level1 = 0.05,level2=0.01,level3=0.001.

# text representation of results, with actual adjusted p-values
write.table(results[[1]], paste0(dir, "/", "results.tsv"), sep = "\t")

    if (bestrr$pval[best] <= pcut) {
        annots <- c(annots, sub("\\d+\\/\\d+ ", "", row.names(bestrr)[best]))
    }
}

mwus <- read.table(paste0(dir, "/", paste("MWU", opt$goDivision, name, sep = "_"), ".", ext), header = T)
best_GOs <- mwus[mwus$name %in% annots, ]
write.table(best_GOs, paste0(dir, "/", "best_go.tsv"), sep = "\t", row.names = FALSE)