comparison GSEA.R @ 0:3a53b8e91ede draft

Uploaded

0:3a53b8e91ede

## How to execute this tool
# $Rscript GSEA.R --input ranked_genes_list.rnk --input Mus_musculus_GSEA_GO_sets_all_symbols_highquality_April_2015.gmt 
# --output GSEA_results.txt --output 

# Send R errors to stderr
options(show.error.messages = F, error = function(){cat(geterrmessage(), file = stderr()); q("no", 1, F)})

# Avoid crashing Galaxy with an UTF8 error on German LC settings
loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")

# Import library
library("getopt")
library("fgsea")
library("Rcpp")
library("data.table")

options(stringAsfactors = FALSE, useFancyQuotes = FALSE)

# Take in trailing command line arguments
args <- commandArgs(trailingOnly = TRUE)

# Get options using the spec as defined by the enclosed list
# Options are read from the default: commandArgs(TRUE)
option_specification = matrix(c(
  'input1', 'i1', 2, 'character',
  'input2', 'i2', 2, 'character',
  'output', 'o', 2, 'character'
), byrow = TRUE, ncol = 4);

# Parse options
options = getopt(option_specification);

# Print options to stderr for debugging
# cat("\n input: ", options$input1)
# cat("\n input: ", options$input2)
# cat("\n output: ", options$output)

# Rank file
ranks <- fread(options$input1, header=T, stringsAsFactors = F)
ranks <- data.frame(ranks)
r <- ranks[abs(ranks$logFC) >= 0.5 & ranks$PValue <= 0.05,]
r <- r[,c(1, 2)]
ranks <- setNames(r[,2], r[,1])

# Pathways database
pathways <- gmtPathways(options$input2)

# running analysis
fgseaRes <- fgsea(pathways, ranks, minSize=10, maxSize=500, nperm=1000)
res <- as.data.frame(fgseaRes[order(pval), ], stringsAsFactors = F)

# save results
write.table(x = res[,1:7], file = options$output, quote = F, row.names = F, sep = "\t")
# 
# topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
# topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n=10), pathway]
# topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
# 
# pdf(paste0(options$output, ".pdf"), width = 8.5, height = 11)
# plotGseaTable(pathways[topPathways], ranks, fgseaRes,gseaParam = 0.5)
# z <- dev.off()