|
0
|
1 ## How to execute this tool
|
|
|
2 # $Rscript GSEA.R --input ranked_genes_list.rnk --input Mus_musculus_GSEA_GO_sets_all_symbols_highquality_April_2015.gmt
|
|
|
3 # --output GSEA_results.txt --output
|
|
|
4
|
|
|
5 # Send R errors to stderr
|
|
|
6 options(show.error.messages = F, error = function(){cat(geterrmessage(), file = stderr()); q("no", 1, F)})
|
|
|
7
|
|
|
8 # Avoid crashing Galaxy with an UTF8 error on German LC settings
|
|
|
9 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
|
|
|
10
|
|
|
11 # Import library
|
|
|
12 library("getopt")
|
|
|
13 library("fgsea")
|
|
|
14 library("Rcpp")
|
|
|
15 library("data.table")
|
|
|
16
|
|
|
17 options(stringAsfactors = FALSE, useFancyQuotes = FALSE)
|
|
|
18
|
|
|
19 # Take in trailing command line arguments
|
|
|
20 args <- commandArgs(trailingOnly = TRUE)
|
|
|
21
|
|
|
22 # Get options using the spec as defined by the enclosed list
|
|
|
23 # Options are read from the default: commandArgs(TRUE)
|
|
|
24 option_specification = matrix(c(
|
|
|
25 'input1', 'i1', 2, 'character',
|
|
|
26 'input2', 'i2', 2, 'character',
|
|
|
27 'output', 'o', 2, 'character'
|
|
|
28 ), byrow = TRUE, ncol = 4);
|
|
|
29
|
|
|
30 # Parse options
|
|
|
31 options = getopt(option_specification);
|
|
|
32
|
|
|
33 # Print options to stderr for debugging
|
|
|
34 # cat("\n input: ", options$input1)
|
|
|
35 # cat("\n input: ", options$input2)
|
|
|
36 # cat("\n output: ", options$output)
|
|
|
37
|
|
|
38 # Rank file
|
|
|
39 ranks <- fread(options$input1, header=T, stringsAsFactors = F)
|
|
|
40 ranks <- data.frame(ranks)
|
|
|
41 r <- ranks[abs(ranks$logFC) >= 0.5 & ranks$PValue <= 0.05,]
|
|
|
42 r <- r[,c(1, 2)]
|
|
|
43 ranks <- setNames(r[,2], r[,1])
|
|
|
44
|
|
|
45 # Pathways database
|
|
|
46 pathways <- gmtPathways(options$input2)
|
|
|
47
|
|
|
48 # running analysis
|
|
|
49 fgseaRes <- fgsea(pathways, ranks, minSize=10, maxSize=500, nperm=1000)
|
|
|
50 res <- as.data.frame(fgseaRes[order(pval), ], stringsAsFactors = F)
|
|
|
51
|
|
|
52 # save results
|
|
|
53 write.table(x = res[,1:7], file = options$output, quote = F, row.names = F, sep = "\t")
|
|
|
54 #
|
|
|
55 # topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
|
|
|
56 # topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n=10), pathway]
|
|
|
57 # topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
|
|
|
58 #
|
|
|
59 # pdf(paste0(options$output, ".pdf"), width = 8.5, height = 11)
|
|
|
60 # plotGseaTable(pathways[topPathways], ranks, fgseaRes,gseaParam = 0.5)
|
|
|
61 # z <- dev.off() |
