Mercurial > repos > dktanwar > fgsea_test
view GSEA.R @ 0:3a53b8e91ede draft
Uploaded
author | dktanwar |
---|---|
date | Mon, 18 Dec 2017 21:09:41 -0500 |
parents | |
children |
line wrap: on
line source
## 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()