# HG changeset patch
# User dktanwar
# Date 1513649381 18000
# Node ID 3a53b8e91ede6630052726335f3720cb65c9f3e0

Uploaded

diff -r 000000000000 -r 3a53b8e91ede GSEA.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/GSEA.R	Mon Dec 18 21:09:41 2017 -0500
@@ -0,0 +1,61 @@
+## 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()
\ No newline at end of file