Gene Set Enrichment Analysis with LIGER

In liger: Lightweight Iterative Geneset Enrichment

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.

Examples

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.

Testing individual gene sets

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.

Testing multiple gene sets

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 Session Info

```r sessionInfo() ````

liger documentation built on Jan. 25, 2021, 9:07 a.m.