Gene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically significant, concordant differences between two biological states. The original algorithm is detailed in Subramanian, Tamayo, et al. with Java implementations available through the Broad Institute.
The liger
package provides a lightweight R implementation of this enrichment test on a list of values. Given a list of values, such as p-values or log-fold changes derived from differential expression analysis or other analyses comparing biological states, this package enables you to test a priori defined set of genes for enrichment to enable interpretability of highly significant or high fold-change genes.
Consider an example, simulated dataset.
library(liger) # load gene set data("org.Hs.GO2Symbol.list") # get universe universe <- unique(unlist(org.Hs.GO2Symbol.list)) # get a gene set gs <- org.Hs.GO2Symbol.list[[1]] # fake dummy example where everything in gene set is perfectly enriched vals <- rnorm(length(universe), 0, 10) names(vals) <- universe vals[gs] <- rnorm(length(gs), 100, 10) head(vals) # look at vals
Here, vals
can be seen as representing a list of log-fold changes derived from differential expression analysis on samples in two biological states. We want to interpret the set of differentially expressed genes with high positive fold changes using gene set enrichment analysis.
To test for enrichment of a particular gene set:
names(org.Hs.GO2Symbol.list)[[1]] gs # look at gs gsea(values=vals, geneset=gs, mc.cores=1, plot=TRUE, n.rand=500)
In this simulation, we created vals
such that gs
was obviously enriched. And indeed, we see that this gene set exhibits significant enrichment.
Now to test for enrichment of another gene set:
gs.new <- org.Hs.GO2Symbol.list[[2]] names(org.Hs.GO2Symbol.list)[[2]] head(gs.new) # look at gs.new gsea(values=vals, geneset=gs.new, mc.cores=1, n.rand=500)
In this simulation, we created vals
such that gs.new
was obviously not enriched. And indeed, we see that this gene set does not exhibit significant enrichment.
If we simulate a more ambiguous case:
# add some noise vals[sample(1:length(universe), 1000)] <- rnorm(1000, 100, 10) # test previously perfectly enriched gene set again gs <- org.Hs.GO2Symbol.list[[1]] gsea(values=vals, geneset=gs, mc.cores=1, n.rand=500)
The enrichment plots and p-values are affected as expected.
We can also test a number of gene sets:
bulk.gsea(values=vals, set.list=org.Hs.GO2Symbol.list[1:5], mc.cores=1, n.rand=500)
To save on computation time, we can also iterative assess significance:
iterative.bulk.gsea(values=vals, set.list=org.Hs.GO2Symbol.list[1:5], mc.cores=1, n.rand=500)
```r sessionInfo() ````
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.