The sparrow package facilitates the use of gene sets in the analysis of high throughput genomics data. It provides simple execution and comparison of several GSEA approaches through a unified interface within the user’s workspace or interactively via a shiny application provided by the sparrow.shiny package. This package also provides an easy wrapper to single sample gene set scoring and geneset-centric heatmaps for visualization. sparrow package version: 1.12.0
sparrow 1.12.0
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:
seas
function is a wrapper that orchestrates the execution of any
number of user-specified gene set enrichment analyses (GSEA) over a particular
experimental contrast of interest. This will create a SparrowResult
object which stores the results of each GSEA method internally, allowing
for easy query and retrieval.{sparrow.shiny}
package provides an explore
function, which is invoked on SparrowResult
objects returned from a call to
seas
. The shiny application facilitates interactive exploration of these
GSEA results. This application can also be deployed to a shiny server and can
be initialized by uploading a serialized SparrowResult
*.rds
file.ora()
which wraps the biased
enrichment functionality found within limma::kegga
and generalizes it to
work against data.frame inputs with arbitrary genesets.scoreSingleSamples
function is a wrapper that enables the user to
generate single sample gene set scores using a variety of different
single sample gene set scoring methods.BiocSet
s
from widely used databases, like getMSigCollection()
(MSigDB),
getKeggCollection()
(KEGG), getPantherCollection()
(PANTHER database), and getReactomeCollection()
(Reactome) with support for different organisms and identifier
types (partially).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:
EList
, DGEList
, or expression matrix)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.
Let’s begin by setting up our work environment for exploratory analysis using the sparrow package.
library(sparrow)
library(magrittr)
library(dplyr)
library(ggplot2)
library(ComplexHeatmap)
library(circlize)
library(edgeR)
library(data.table)
theme_set(theme_bw())
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.
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.
Below you’ll find the $targets
data.frame of the voomed EList
## 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).
We will identify the genes and genesets that are differentially expressed
between the basal and her2 subtypes. The vm
object has already been voom
d
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:
## Contrasts
## Levels BvH
## Basal 1
## Her2 -1
## LumA 0
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:
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.
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.
To retrieve a BiocSet
of these same collections, you could do:
You can view a table of the gene sets defined inside a GeneSetDb
(gdb
)object
via its geneSets(gdb)
accessor:
## 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
## 6 C2 ABDULRAHMAN_KIDNEY_CANCER_VHL_DN FALSE 13
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 …
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:
## 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$t.tt, comp$t.mg)
## [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.
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.
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
):
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) %>%
head
## 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.
geneSet(mg, name = 'HALLMARK_WNT_BETA_CATENIN_SIGNALING') %>%
select(symbol, entrez_id, logFC, pval, padj) %>%
head()
## 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.
{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()
.
Boxplot
Density
GSEA plot
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: