1 Overview

The {sparrow} package was built to facilitate the use of gene sets in the analysis of high throughput genomics data (primarily RNA-seq). It does so by providing these top-line functionalities:

The initial GSEA methods that sparrow wrapped were the ones provided by limma and edgeR. As such, many analyses using sparrow expect you to re-use the same data objects used for differential expression analysis, namely:

Other methods only require the user to provide a ranked vector of statistics that represent some differential expression statistic per gene, and the GSEA is performed by analyzing the ranks of genes within this vector.

The user can invoke one seas() call that can orchestrate multiple analyses of any type.

Currently supported gene set enrichment methods include:

##            method test_type package
## 1          camera preranked   limma
## 2        cameraPR preranked   limma
## 3           fgsea preranked   fgsea
## 4             ora       ora     ora
## 5             fry preranked   limma
## 6           roast preranked   limma
## 7           romer preranked   limma
## 8           goseq       ora   goseq
## 9     geneSetTest preranked   limma
## 10          logFC preranked   limma
## 11 svdGeneSetTest      meta sparrow

When using these methods in analyses that lead to publication, please cite the original papers that developed these methods and cite sparrow when its functionality assisted in your interpretation and analysis.

The sparrow package provides a small example expression dataset extracted from the TCGA BRCA dataset, which is available via the exampleExpressionSet function. In this vignette we will explore differential expression and gene set enrichment analysis by examining differences between basal and her2 PAM50 subtypes.

2 Standard Workflow

Let’s begin by setting up our work environment for exploratory analysis using the sparrow package.


Internally, sparrow leverages the data.table package for fast indexing and manipulation over data.frames. All functions that return data.frame looking objects back have converted it from an data.table prior to return. All such functions take an as.dt argument, which is set to FALSE by default that controls this behavior. If you want {sparrow} to return a data.table back to you from some function, try adding an as.dt = TRUE argument to the end of the function call.

2.1 Data Setup

sparrow is most straightforward to use when our data objects and analysis are performed with either the edgeR or voom/limma pipelines and when we use standard gene identifiers (like esnemble) as rownames() to these objects.

The exampleExpressionSet function gives us just such an object. We call it below in a manner that gives us an object that allows us to explore expression differences between different subtypes of breast cancer.

vm <- exampleExpressionSet(dataset = "tumor-subtype", do.voom = TRUE)

Below you’ll find the $targets data.frame of the voomed EList

vm$targets %>%
  select(Patient_ID, Cancer_Status, PAM50subtype)
##                                Patient_ID Cancer_Status PAM50subtype
## TCGA-A2-A0CM-01A-31R-A034-07 TCGA-A2-A0CM         tumor        Basal
## TCGA-BH-A0RX-01A-21R-A084-07 TCGA-BH-A0RX         tumor        Basal
## TCGA-BH-A18Q-01A-12R-A12D-07 TCGA-BH-A18Q         tumor        Basal
## TCGA-B6-A0RU-01A-11R-A084-07 TCGA-B6-A0RU         tumor        Basal
## TCGA-BH-A18P-01A-11R-A12D-07 TCGA-BH-A18P         tumor         Her2
## TCGA-C8-A275-01A-21R-A16F-07 TCGA-C8-A275         tumor         Her2
## TCGA-C8-A12Z-01A-11R-A115-07 TCGA-C8-A12Z         tumor         Her2
## TCGA-A2-A0T1-01A-21R-A084-07 TCGA-A2-A0T1         tumor         Her2
## TCGA-AC-A3OD-01A-11R-A21T-07 TCGA-AC-A3OD         tumor         LumA
## TCGA-AN-A0XS-01A-22R-A109-07 TCGA-AN-A0XS         tumor         LumA
## TCGA-A2-A0EM-01A-11R-A034-07 TCGA-A2-A0EM         tumor         LumA
## TCGA-AR-A24O-01A-11R-A169-07 TCGA-AR-A24O         tumor         LumA
## TCGA-D8-A4Z1-01A-21R-A266-07 TCGA-D8-A4Z1         tumor         LumA

Note that there are many tutorials online that outline how to generate expression matrices for use with differential expression and analysis, such as the one that is returned from the exampleExpressionSet function. Summarizing assay data into such a format is out of scope for this vignette, but you can reference the airway vignette for full details (among others).

2.2 Data Analysis

We will identify the genes and genesets that are differentially expressed between the basal and her2 subtypes. The vm object has already been voomd using this design:

##                              Basal Her2 LumA
## TCGA-A2-A0CM-01A-31R-A034-07     1    0    0
## TCGA-BH-A0RX-01A-21R-A084-07     1    0    0
## TCGA-BH-A18Q-01A-12R-A12D-07     1    0    0
## TCGA-B6-A0RU-01A-11R-A084-07     1    0    0
## TCGA-BH-A18P-01A-11R-A12D-07     0    1    0
## TCGA-C8-A275-01A-21R-A16F-07     0    1    0
## TCGA-C8-A12Z-01A-11R-A115-07     0    1    0
## TCGA-A2-A0T1-01A-21R-A084-07     0    1    0
## TCGA-AC-A3OD-01A-11R-A21T-07     0    0    1
## TCGA-AN-A0XS-01A-22R-A109-07     0    0    1
## TCGA-A2-A0EM-01A-11R-A034-07     0    0    1
## TCGA-AR-A24O-01A-11R-A169-07     0    0    1
## TCGA-D8-A4Z1-01A-21R-A266-07     0    0    1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$PAM50subtype
## [1] "contr.treatment"

We can test for differences between basla and her2 subtypes using the following contrast:

(cm <- makeContrasts(BvH=Basal - Her2, levels=vm$design))
##        Contrasts
## Levels  BvH
##   Basal   1
##   Her2   -1
##   LumA    0

2.2.1 Differential Gene Expression

In this section, we first show you the straightforward analysis you would do if you were only testing for differential gene expression.

With the data we have at hand, you would simply do the following:

fit <- lmFit(vm, vm$design) %>% %>%
tt <- topTable(fit, 'BvH', n=Inf,'none')

2.2.2 Gene Set Enrichment Analysis

Given that we now have all of the pieces of data required for a differential expression analysis, performing GSEA is trivial using the seas wrapper function. We simply need to now define (1) the battery of gene sets we want to test against, and (2) the GSEA methods we want to explore.

2.2.3 Gene Sets to Test

The sparrow package provides a GeneSetDb class to store collections of gene sets. The GeneSetDb object is used heavily for the internal functionality of {sparrow}, however you can provide sparrow with collections of gene sets using other containers from the bioconductor universe, like a BiocSet::BiocSet or a GSEABase::GeneSetCollection. This package provides convenience methods to convert between these different types of gene set containers. Please refer to The GeneSetDb Class section for more details.

The {sparrow} package also provides convenience methods to retrieve gene set collections from different sourckes, like MSigDB, PANTHER, KEGG, etc. These methods are named using the following pattern: get<CollectionName>Collection() to return a BiocSet with the gene sets from the collection, or get<CollectionName>GeneSetDb() to get a GeneSetDb of the same.

We’ll use the getMSigGeneSetDb convenience function provided by the sparrow package to load the hallmark ("h") and c2 (curated) ("c2") gene set collections from MSigDB.

gdb <- getMSigGeneSetDb(c("h", "c2"), "human", id.type = "entrez")

To retrieve a BiocSet of these same collections, you could do:

bsc <- getMSigCollection(c("h", "c2"), "human", id.type = "entrez")

You can view a table of the gene sets defined inside a GeneSetDb (gdb)object via its geneSets(gdb) accessor:

geneSets(gdb) %>%
  head %>%
##   collection                             name active  N
## 1         C2         ABBUD_LIF_SIGNALING_1_DN  FALSE 28
## 2         C2         ABBUD_LIF_SIGNALING_1_UP  FALSE 43
## 3         C2         ABBUD_LIF_SIGNALING_2_DN  FALSE  7
## 4         C2         ABBUD_LIF_SIGNALING_2_UP  FALSE 13
## 5         C2       ABDELMOHSEN_ELAVL4_TARGETS  FALSE 16

2.3 Running sparrow

Performing multiple gene set enrichment analyses over your contrast of interest simply requires you to provide a GeneSetDb (or BiocSet) object along with your data and an enumeration of the methods you want to use in your analysis.

The call to seas() will perform these analyses and return a SparrowResult object which you can then use for downstream analysis.

mg <- seas(
  vm, gdb, c('camera', 'fry', 'ora'),
  design = vm$design, contrast = cm[, 'BvH'],
  # these parameters define which genes are differentially expressed
  feature.max.padj = 0.05, feature.min.logFC = 1,
  # for camera:
  inter.gene.cor = 0.01,
  # specifies the numeric covariate to bias-correct for
  # "size" is found in the vm$genes data.frame, which makes its way to the
  # internal DGE statistics table ... more on that later
  feature.bias = "size")

We will unpack the details of the seas() call shortly …

2.4 Implicit Differential Expression

First, let’s note that in addition to running a plethora of GSEA’s over our data we’ve also run a standard differential expression analysis. If you’ve passed a matrix, ExpressionSet or EList into seas(), a limma-based lmFit %>% (eBayes|treat) %>% (topTable|topTreat) pipeline was run. If a DGEList was passed, then seas utilizes the edgeR-based glmQLFit %>% (glmQLFTest | glmTreat) %>% topTags pipeline.

The result of the internally run differential expression analysis is accessible via a call to logFC function on the SparrowResult object:

lfc <- logFC(mg)
lfc %>%
  select(symbol, entrez_id, logFC, t, pval, padj) %>%
##         symbol entrez_id       logFC           t      pval      padj
## 1         A1BG         1  0.67012895  1.07951394 0.2982819 0.6858344
## 2          ADA       100  0.53844094  0.92401125 0.3708544 0.7415607
## 3         CDH2      1000 -0.08180996 -0.09901074 0.9225083 0.9795974
## 4         AKT3     10000  0.58338138  1.29502525 0.2158892 0.6125318
## 5 LOC100009676 100009676 -0.09581391 -0.26985709 0.7911366 0.9398579
## 6         MED6     10001  0.04505155  0.15082239 0.8822288 0.9701384

We can confirm that the statistics generated internally in seas() mimic our explicit analysis above by verifying that the t-statistics generated by both approaches are identical.

comp <- tt %>%
  select(entrez_id, logFC, t, pval=P.Value, padj=adj.P.Val) %>%
  inner_join(lfc, by='entrez_id', suffix=c('.tt', '.mg'))
all.equal(comp$, comp$
## [1] TRUE

The internally performed differential expression analysis within the seas() call can be customized almost as extensively as an explicitly performed analysis that you would run using limma or edgeR by sending more parameters through seas()’s ... argument.

See the Custom Differential Expression section further in the vignette as well as the help available in ?calculateIndividualLogFC (which is called inside the seas() function) for more information.

2.5 Explicit GSEA

We also have the results of all the GSEA analyses that we specified to our seas call via the methods parameter.

## SparrowResult (max FDR by collection set to 0.20%)
## --------------------------------------------------- 
##    collection   method geneset_count sig_count sig_up sig_down
## 1          C2   camera          6150       349    206      143
## 2           H   camera            50         6      5        1
## 3          C2      fry          6150        95     33       62
## 4           H      fry            50         0      0        0
## 5          C2      ora          6150        96     33       63
## 6           H      ora            50         3      1        2
## 7          C2 ora.down          6150        73      6       67
## 8           H ora.down            50         2      0        2
## 9          C2   ora.up          6150        24     21        3
## 10          H   ora.up            50         0      0        0

The table above enumerates the different GSEA methods run over each geneset collection in the rows. The columns enumerate the number of genesets that the collection has in total (geneset_count), and how many were found significant at a given FDR, which is set to 20% by default. The show command for the SparrowResult object simply calls the tabulateResults() function, which you can call directly with the value of max.p that you might find more appropriate.

2.6 Exploring Results

GSEA results can be examined interactively via the command line, or via a shiny application. You can use the resultNames function to find out what GSEA methods were run, and therefore available to you, within the the SparrowResult object:

## [1] "camera"   "fry"      "ora"      "ora.down" "ora.up"

Note that when running an “over representation analysis” "ora" (or "goseq"), it will be run three different ways. The tests will be run first by testing all differentially expressed genes that meet a given set of min logFC and max FDR thresholds, then separately for only genes that go up in your contrast, and a third time for only the genes that go down.

The individual gene set statistics generated by each method are available via the result function (or several can be returned with results):

cam.res <- result(mg, 'camera')
cam.go.res <- results(mg, c('camera', 'ora.up'))

You can identify genesets with the strongest enrichment by filtering and sorting against the appropriate columns. We can, for instance, identify which hallmark gene sets show the strongest enrichment as follows:

cam.res %>%
  filter(padj < 0.1, collection == 'H') %>%
  arrange(desc(mean.logFC)) %>%
  select(name, n, mean.logFC, padj) %>%
##                                 name   n mean.logFC         padj
## 1            HALLMARK_MYC_TARGETS_V2  58  0.4461105 0.0002612790
## 2 HALLMARK_INTERFERON_ALPHA_RESPONSE  96  0.3916716 0.0874709010
## 3               HALLMARK_E2F_TARGETS 200  0.3465703 0.0001892151
## 4            HALLMARK_MYC_TARGETS_V1 200  0.2092836 0.0234431144

You can also list the members of a geneset and their individual differential expression statistics for the contrast under test using the geneSet function.

  select(symbol, entrez_id, logFC, pval, padj) %>% 
##   symbol entrez_id      logFC       pval      padj
## 1  HDAC5     10014  0.8984691 0.02253974 0.2522754
## 2 CSNK1E      1454 -0.1793725 0.52104817 0.8317753
## 3 CTNNB1      1499  0.2577554 0.54741640 0.8467181
## 4   JAG1       182  0.7293432 0.02496690 0.2625306
## 5   DVL2      1856  0.4921509 0.24186744 0.6362028
## 6   DKK1     22943  0.6567652 0.66735589 0.8982828

The results provided in the table generated from a call to geneSet are independant of GSEA method. The statistics appended to the gene set members are simply the ones generated from a differential expression analysis.

2.7 Plotting

{sparrow} provides a number of interactive plotting facilities to explore the enrichment of a single geneset under the given contrast. In the boxplots and density plots shown below, the log fold changes (logFCs) (or t-statistics) for all genes under the contrast are visualized in the “background” set, and these same values are shown for the desired geneset under the “geneset” group.

The logFC (or t-statistics) of the genes in the gene set are plotted as points, which allow you to hover to identify the identity of the genes that land in the regions of the distributions you care about.

Including interactive plots increases the size of the vignette’s by a lot and will be rejected by the bioconductor build servers, so all plots included in this vignette are static snapshots of the javascript enabled plots you would normally get from iplot().


      type = 'boxplot', value = 'logFC')

boxplot of geneset log2FC’s


      type = 'density', value = 'logFC')

density plot of geneset log2FC’s

GSEA plot

      type = 'gsea', value = 'logFC')

gsea plot of geneset log2FC’s

2.7.1 Interactive Exploration

A sister {sparrow.shiny} package is available that can be used to interactively explore SparrowResult objects to help you try to make sense of the enrichment hits you get (or not!). The application can be invoked as follows: