fgsea
is an R-package for fast preranked gene set enrichment analysis (GSEA). The performance is achieved by using an algorithm for cumulative GSEA-statistic calculation. This allows to reuse samples between different gene set sizes. See the preprint for algorithmic details.
Loading example pathways and gene-level statistics:
data(examplePathways)
data(exampleRanks)
Running fgsea:
fgseaRes <- fgsea(pathways = examplePathways,
stats = exampleRanks,
minSize=15,
maxSize=500,
nperm=10000)
The resulting table contains enrichment scores and p-values:
head(fgseaRes[order(pval), ])
## pathway pval padj ES
## 1: 5990980_Cell_Cycle 0.0001236552 0.002275657 0.5388497
## 2: 5990979_Cell_Cycle,_Mitotic 0.0001260239 0.002275657 0.5594755
## 3: 5991210_Signaling_by_Rho_GTPases 0.0001320132 0.002275657 0.4238512
## 4: 5991454_M_Phase 0.0001377790 0.002275657 0.5576247
## 5: 5991023_Metabolism_of_carbohydrates 0.0001396648 0.002275657 0.4944766
## 6: 5991209_RHO_GTPase_Effectors 0.0001408649 0.002275657 0.5248796
## NES nMoreExtreme size leadingEdge
## 1: 2.686525 0 369 66336,66977,12442,107995,66442,19361,
## 2: 2.753317 0 317 66336,66977,12442,107995,66442,12571,
## 3: 2.016900 0 231 66336,66977,20430,104215,233406,107995,
## 4: 2.564862 0 173 66336,66977,12442,107995,66442,52276,
## 5: 2.250319 0 160 11676,21991,15366,58250,12505,20527,
## 6: 2.380726 0 157 66336,66977,20430,104215,233406,107995,
It takes about ten seconds to get results with significant hits after FDR correction:
sum(fgseaRes[, padj < 0.01])
## [1] 73
One can make an enrichment plot for a pathway:
plotEnrichment(examplePathways[["5991130_Programmed_Cell_Death"]],
exampleRanks) + labs(title="Programmed Cell Death")
Or make a table plot for a bunch of selected pathways:
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))
plotGseaTable(examplePathways[topPathways], exampleRanks, fgseaRes,
gseaParam = 0.5)
Please, be aware that fgsea
function takes about O(nk^{3/2}) time, where n is number of permutations and k is a maximal size of the pathways. That means that setting maxSize
parameter with a value of ~500 is strongly recommended.
Also, fgsea
is parallelized using BiocParallel
package. By default the first registered backend returned by bpparam()
is used. To tweak the parallelization one can either specify BPPARAM
parameter used for bclapply
of set nproc
parameter, which is a shorthand for setting BPPARAM=MulticoreParam(workers = nproc)
.
For convenience there is reactomePathways
function that obtains pathways from Reactome for given set of genes. Package reactome.db
is required to be installed.
pathways <- reactomePathways(names(exampleRanks))
fgseaRes <- fgsea(pathways, exampleRanks, nperm=1000, maxSize=500)
head(fgseaRes)
## pathway
## 1: Meiotic Synapsis
## 2: Rora activates gene expression
## 3: Bmal1:Clock,Npas2 activates circadian gene expression
## 4: Translocation of Glut4 to the Plasma Membrane
## 5: Endocrine-committed (Ngn3+) progenitor cells
## 6: Late stage (branching morphogenesis) pancreatic bud precursor cells
## pval padj ES NES nMoreExtreme size
## 1: 0.5527638 0.7975585 0.2885754 0.9397103 329 27
## 2: 0.8375000 0.9271420 -0.3087414 -0.6713139 401 5
## 3: 0.4351852 0.7344723 0.4209054 1.0459978 234 9
## 4: 0.6928105 0.8771748 0.2387284 0.8444781 423 39
## 5: 0.4949698 0.7687686 0.6477746 1.0155043 245 2
## 6: 0.9623762 0.9839168 -0.3460577 -0.5580225 485 2
## leadingEdge
## 1: 15270,12189,71846,19357
## 2: 20787,328572,12753,11865
## 3: 20893,59027,19883
## 4: 17918,19341,20336,22628,22627,20619,
## 5: 18088,18506
## 6: 15205,11925
One can also start from .rnk
and .gmt
files as in original GSEA:
rnk.file <- system.file("extdata", "naive.vs.th1.rnk", package="fgsea")
gmt.file <- system.file("extdata", "mouse.reactome.gmt", package="fgsea")
Loading ranks:
ranks <- read.table(rnk.file,
header=TRUE, colClasses = c("character", "numeric"))
ranks <- setNames(ranks$t, ranks$ID)
str(ranks)
## Named num [1:12000] -63.3 -49.7 -43.6 -41.5 -33.3 ...
## - attr(*, "names")= chr [1:12000] "170942" "109711" "18124" "12775" ...
Loading pathways:
pathways <- gmtPathways(gmt.file)
str(head(pathways))
## List of 6
## $ 1221633_Meiotic_Synapsis : chr [1:64] "12189" "13006" "15077" "15078" ...
## $ 1368092_Rora_activates_gene_expression : chr [1:9] "11865" "12753" "12894" "18143" ...
## $ 1368110_Bmal1:Clock,Npas2_activates_circadian_gene_expression : chr [1:16] "11865" "11998" "12753" "12952" ...
## $ 1445146_Translocation_of_Glut4_to_the_Plasma_Membrane : chr [1:55] "11461" "11465" "11651" "11652" ...
## $ 186574_Endocrine-committed_Ngn3+_progenitor_cells : chr [1:4] "18012" "18088" "18506" "53626"
## $ 186589_Late_stage_branching_morphogenesis_pancreatic_bud_precursor_cells: chr [1:4] "11925" "15205" "21410" "246086"
And runnig fgsea:
fgseaRes <- fgsea(pathways, ranks, minSize=15, maxSize=500, nperm=1000)
head(fgseaRes)
## pathway
## 1: 1221633_Meiotic_Synapsis
## 2: 1445146_Translocation_of_Glut4_to_the_Plasma_Membrane
## 3: 442533_Transcriptional_Regulation_of_Adipocyte_Differentiation_in_3T3-L1_Pre-adipocytes
## 4: 508751_Circadian_Clock
## 5: 5334727_Mus_musculus_biological_processes
## 6: 573389_NoRC_negatively_regulates_rRNA_expression
## pval padj ES NES nMoreExtreme size
## 1: 0.52810903 0.7033452 0.2885754 0.9468611 309 27
## 2: 0.68561873 0.8294954 0.2387284 0.8436265 409 39
## 3: 0.09133489 0.2267892 -0.3640706 -1.3687852 38 31
## 4: 0.78685613 0.8811827 0.2516324 0.7267117 442 17
## 5: 0.38814815 0.5802419 0.2469065 1.0552212 261 106
## 6: 0.40852575 0.6037390 0.3607407 1.0418153 229 17
## leadingEdge
## 1: 15270,12189,71846,19357
## 2: 17918,19341,20336,22628,22627,20619,
## 3: 20602,327987,59024,67381,70208,12537,
## 4: 20893,59027,19883
## 5: 60406,19361,15270,20893,12189,68240,
## 6: 60406,20018,245688,20017