Using msigdf with enrichment tools

Abstract

Recipes for reshaping the msigdf tibbles into the formats expected by common gene set enrichment tools: fgsea, clusterProfiler, and anything else that wants either a named list of gene sets or a two-column TERM2GENE frame.

library(dplyr)
library(msigdf)

Scope

msigdf ships the gene sets only. It deliberately does not depend on any enrichment package, so fgsea and clusterProfiler are not installed alongside it.

The chunks that reshape msigdf data are evaluated below, so the output you see is real. The chunks that call an enrichment package are marked eval=FALSE and are shown as recipes to copy — install the relevant Bioconductor package first:

# install.packages("BiocManager")
BiocManager::install(c("fgsea", "clusterProfiler"))

A named list of gene sets

fgsea::fgsea() and several other tools expect gene sets as a named list: one element per gene set, each holding a character vector of gene identifiers. split() produces exactly that.

Using the human hallmark collection (category_code == "h"):

hallmark <- msigdf.human %>%
  filter(category_code == "h")

pathways <- split(hallmark$symbol, hallmark$geneset)

length(pathways)

## [1] 50

lengths(pathways) %>% head()

##        HALLMARK_ADIPOGENESIS HALLMARK_ALLOGRAFT_REJECTION 
##                          200                          200 
##   HALLMARK_ANDROGEN_RESPONSE        HALLMARK_ANGIOGENESIS 
##                          101                           36 
##     HALLMARK_APICAL_JUNCTION      HALLMARK_APICAL_SURFACE 
##                          200                           44

pathways[["HALLMARK_NOTCH_SIGNALING"]] %>% head(10)

##  [1] "APH1A"  "ARRB1"  "CCND1"  "CUL1"   "DLL1"   "DTX1"   "DTX2"   "DTX4"  
##  [9] "FBXW11" "FZD1"

Feed that straight to fgsea, together with your own named vector of ranking statistics:

library(fgsea)

# `stats` is a named numeric vector: names are gene symbols, values are the
# ranking metric (log fold change, signed -log10 p-value, ...).
res <- fgsea(pathways = pathways, stats = stats)

res %>%
  as_tibble() %>%
  arrange(padj) %>%
  select(pathway, pval, padj, NES, size)

Any collection works the same way. To restrict to one sub-collection, filter on category_subcode as well:

reactome <- msigdf.human %>%
  filter(category_code == "c2", category_subcode == "cp.reactome")

reactome_pathways <- split(reactome$symbol, reactome$geneset)
length(reactome_pathways)

## [1] 1839

A TERM2GENE data frame

clusterProfiler::enricher() and GSEA() take a two-column data frame: gene set name first, gene identifier second. That is the msigdf layout minus the collection columns.

term2gene <- msigdf.human %>%
  filter(category_code == "h") %>%
  select(term = geneset, gene = symbol)

term2gene

## # A tibble: 7,322 × 2
##    term                  gene 
##    <chr>                 <chr>
##  1 HALLMARK_ADIPOGENESIS ABCA1
##  2 HALLMARK_ADIPOGENESIS ABCB8
##  3 HALLMARK_ADIPOGENESIS ACAA2
##  4 HALLMARK_ADIPOGENESIS ACADL
##  5 HALLMARK_ADIPOGENESIS ACADM
##  6 HALLMARK_ADIPOGENESIS ACADS
##  7 HALLMARK_ADIPOGENESIS ACLY 
##  8 HALLMARK_ADIPOGENESIS ACO2 
##  9 HALLMARK_ADIPOGENESIS ACOX1
## 10 HALLMARK_ADIPOGENESIS ADCY6
## # ℹ 7,312 more rows

Then:

library(clusterProfiler)

# Over-representation on a character vector of gene symbols.
ora <- enricher(gene = my_genes, TERM2GENE = term2gene)

# Or GSEA on a sorted, named numeric vector.
gsea <- GSEA(geneList = sorted_stats, TERM2GENE = term2gene)

Mouse

The mouse tibble has the same shape. Note the collection codes differ: hallmark is mh rather than h, and the numbered collections run m1-m8.

msigdf.mouse %>%
  distinct(category_code) %>%
  arrange(category_code)

## # A tibble: 7 × 1
##   category_code
##   <chr>        
## 1 m1           
## 2 m2           
## 3 m3           
## 4 m5           
## 5 m7           
## 6 m8           
## 7 mh

mouse_hallmark <- msigdf.mouse %>%
  filter(category_code == "mh")

mouse_pathways <- split(mouse_hallmark$symbol, mouse_hallmark$geneset)

length(mouse_pathways)

## [1] 50

mouse_pathways[["HALLMARK_NOTCH_SIGNALING"]] %>% head(10)

##  [1] "Aph1a"  "Arrb1"  "Ccnd1"  "Cul1"   "Dll1"   "Dtx1"   "Dtx2"   "Dtx4"  
##  [9] "Fbxw11" "Fzd1"

Mouse gene sets carry mouse symbols, so no ortholog mapping is needed if your data is already mouse.

Linking results back to MSigDB

msigdf.urls and msigdf.mouse.urls map every gene set name to its page on the MSigDB website. Join enrichment output against them to get clickable links in a report.

top_sets <- msigdf.human %>%
  filter(category_code == "h") %>%
  count(geneset, name = "n_genes") %>%
  arrange(n_genes) %>%
  head(5)

top_sets %>%
  inner_join(msigdf.urls, by = "geneset") %>%
  select(geneset, n_genes, url)

## # A tibble: 5 × 3
##   geneset                             n_genes url                               
##   <chr>                                 <int> <chr>                             
## 1 HALLMARK_NOTCH_SIGNALING                 32 http://software.broadinstitute.or…
## 2 HALLMARK_ANGIOGENESIS                    36 http://software.broadinstitute.or…
## 3 HALLMARK_HEDGEHOG_SIGNALING              36 http://software.broadinstitute.or…
## 4 HALLMARK_PANCREAS_BETA_CELLS             40 http://software.broadinstitute.or…
## 5 HALLMARK_WNT_BETA_CATENIN_SIGNALING      42 http://software.broadinstitute.or…

The same join applied to an fgsea result:

res %>%
  as_tibble() %>%
  arrange(padj) %>%
  head(20) %>%
  inner_join(msigdf.urls, by = c("pathway" = "geneset")) %>%
  select(pathway, padj, NES, url)

Session info