---
title: "Reproducible GSEA Benchmarking"
author:
- name: Ludwig Geistlinger
affiliation: Center for Computational Biomedicine, Harvard Medical School
email: ludwig_geistlinger@hms.harvard.edu
package: GSEABenchmarkeR
abstract: >
The _GSEABenchmarkeR_ package implements an extendable framework for
reproducible evaluation of set- and network-based methods for enrichment
analysis of gene expression data. This includes support for the efficient
execution of these methods on comprehensive real data compendia (microarray
and RNA-seq) using parallel computation on standard workstations and
institutional computer grids. Methods can then be assessed with respect to
runtime, statistical significance, and relevance of the results for the
phenotypes investigated.
output:
BiocStyle::html_document:
toc: true
toc_depth: 2
vignette: >
% \VignetteIndexEntry{Reproducible GSEA Benchmarking}
% \VignetteEngine{knitr::rmarkdown}
---
```{r setup, echo=FALSE}
suppressPackageStartupMessages({
library(GSEABenchmarkeR)
library(EnrichmentBrowser)
})
```
# Purpose of the package
The purpose of the `r Biocpkg("GSEABenchmarkeR")` package is to compare the performance
of different methods for gene set enrichment analysis across many gene expression
datasets.
Users interested in conducting gene set enrichment analysis for a specific dataset
of choice are recommended to use the `r Biocpkg("EnrichmentBrowser")` package instead.
In other words,
- if you are interested in analysing a particular microarray or RNA-seq dataset,
e.g. a case-control study where you want to find out which GO terms / KEGG pathways
are enriched for differentially expressed genes, i.e. your primary goal is _biological
interpretation of a specific dataset under investigation_, then use the
`r Biocpkg("EnrichmentBrowser")` package.
- if you want to assess the performance (runtime, type I error rate, etc) of
different enrichment methods across many datasets and in certain simulated setups -
i.e. your primary goal is to _understand methodological aspects and compare methods against each other_,
then use the `r Biocpkg("GSEABenchmarkeR")` package.
# Setup
Although gene set enrichment analysis (GSEA) has become an integral part of
high-throughput gene expression data analysis, the assessment of enrichment
methods remains rudimentary and *ad hoc*.
In the absence of suitable gold standards, the evaluation is commonly restricted
to selected data sets and biological reasoning on the relevance of resulting
enriched gene sets.
However, this is typically incomplete and biased towards a novel method being
presented.
As the evaluation of GSEA methods is thus typically based on self-defined
standards, [Mitrea et al. (2013)](https://doi.org/10.3389/fphys.2013.00278)
identified the lack of gold standards for consistent assessment and comparison
of enrichment methods as a major bottleneck.
Furthermore, it is often cumbersome to reproduce existing assessments for
additional methods, as this typically involves considerable effort of data
processing and method collection.
Leveraging the representative and extendable collection of
enrichment methods available in the `r Biocpkg("EnrichmentBrowser")` package,
the `r Biocpkg("GSEABenchmarkeR")` package facilitates efficient execution
of these methods on comprehensive real data compendia.
The compendia are curated collections of microarray and RNA-seq datasets
investigating human diseases (mostly specific cancer types), for which
disease-relevant gene sets have been defined _a priori_.
Consistently applied to these datasets, enrichment methods can then be subjected
to a systematic and reproducible assessment of _(i)_ computational runtime,
_(ii)_ statistical significance, especially how the fraction of
significant gene sets relates to the fraction of differentially expressed genes,
and _(iii)_ phenotype relevance, i.e. whether enrichment methods produce gene set
rankings in which phenotype-relevant gene sets accumulate at the top.
In the following, we demonstrate how the package can be used to
- load specific pre-defined and user-defined data compendia,
- carry out differential expression analysis across datasets,
- apply enrichment methods to multiple datasets, and
- benchmark results with respect to the chosen criteria.
We start by loading the package.
```{r lib}
library(GSEABenchmarkeR)
```
# Expression data sources
The `r Biocpkg("GSEABenchmarkeR")` package implements a general interface for
loading compendia of expression datasets.
This includes
- the pre-defined GEO2KEGG microarray compendium that consists of 42 datasets
investigating a total of 19 different human diseases as collected by
Tarca et al. ([2012](https://doi.org/10.1186/1471-2105-13-136) and
[2013](https://doi.org/10.1371/journal.pone.0079217)),
- the pre-defined TCGA RNA-seq compendium, consisting of datasets from
[The Cancer Genome Atlas](https://cancergenome.nih.gov)
investigating a total of 34 different cancer types, and
- user-defined data from file.
In the following, we describe both pre-defined compendia in more detail and also
demonstrate how user-defined data can be incorporated.
## Microarray compendium
Although RNA-seq (read count data) has become the *de facto* standard for
transcriptomic profiling, it is important to know that many methods for
differential expression and gene set enrichment analysis have been originally
developed for microarray data (intensity measurements).
However, differences in data distribution assumptions (microarray: quasi-normal,
RNA-seq: negative binomial) have made adaptations in differential expression
analysis and, to some extent also in gene set enrichment analysis, necessary.
Nevertheless, the comprehensive collection and curation of microarray data in
online repositories such as [GEO](https://www.ncbi.nlm.nih.gov/geo) still
represent a valuable resource.
In particular, Tarca et al. ([2012](https://doi.org/10.1186/1471-2105-13-136)
and [2013](https://doi.org/10.1371/journal.pone.0079217)) compiled 42 datasets
from GEO, each investigating a human disease for which a specific
[KEGG](http://www.genome.jp/kegg) pathway exists.
These pathways are accordingly defined as the target pathways for the various
enrichment methods when applied to the respective datasets. For instance,
methods are expected to rank the
[Alzheimer's disease pathway](http://www.genome.jp/kegg-bin/show_pathway?hsa05010)
close to the top when applied to
[GSE1297](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1297),
a case-control study of Alzheimer's disease.
Furthermore, Tarca et al. made these datasets available in the _Bioconductor_
packages `r Biocpkg("KEGGdzPathwaysGEO")` and
`r Biocpkg("KEGGandMetacoreDzPathwaysGEO")`.
The `r Biocpkg("GSEABenchmarkeR")` package simplifies access to the compendium
and allows to load it into the workspace via
```{r maComp}
geo2kegg <- loadEData("geo2kegg")
names(geo2kegg)
```
A specific dataset of the compendium can be obtained via
```{r getDatasetProbe}
geo2kegg[[1]]
```
which returns, in this example, an `ExpressionSet` (documented in the
`r Biocpkg("Biobase")` package) that contains expression levels of
22,283 probe sets measured for 16 patients.
To prepare the datasets for subsequent analysis, the `r Biocpkg("GSEABenchmarkeR")`
package provides the function `maPreproc`.
The function invokes `EnrichmentBrowser::probe2gene` on each dataset
to summarize expression levels for probes annotated to the same gene.
Here, we apply the function to the first 5 datasets of the compendium.
```{r maPreproc}
geo2kegg <- maPreproc(geo2kegg[1:5])
```
Now,
```{r getDatasetGene}
geo2kegg[[1]]
```
returns a `r Biocpkg("SummarizedExperiment")` that contains the summarized
expression values for 12,994 genes. Furthermore, sample groups
are defined in the `colData` column **GROUP**, yielding here 7 cases (1) and
9 controls (0).
```{r maGroups}
se <- geo2kegg[[1]]
table(se$GROUP)
```
Note: The `maPreproc` returns datasets consistently mapped to NCBI Entrez
Gene IDs, which is compatible with most downstream applications. However,
mapping to a different ID type such as ENSEMBL IDs or HGNC symbols can also be
done using the function `EnrichmentBrowser::idMap`.
## RNA-seq compendium
The Cancer Genome Atlas ([TCGA](https://cancergenome.nih.gov)) project performed
a molecular investigation of various cancer types on an unprecedented scale
including various genomic high-throughput technologies.
In particular, transcriptomic profiling of the investigated cancer types has
comprehensively been carried out with RNA-seq in tumor and adjacent normal tissue.
Among the various resources that redistribute TCGA data,
[Rahman et al. (2015)](https://doi.org/10.1093/bioinformatics/btv377)
consistently preprocessed the RNA-seq data for 24 cancer types and made the data
available in the GEO dataset
[GSE62944](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944).
The GSE62944 compendium can be loaded using the `loadEData` function,
which provides the datasets ready for subsequent differential expression and
gene set enrichment analysis.
Here, we load the compendium into the workspace using only two of the datasets.
```{r rseqComp}
tcga <- loadEData("tcga", nr.datasets=2)
names(tcga)
```
For example, the breast cancer dataset contains RNA-seq read counts for roughly
20,000 genes measured in 1,119 tumor (1) and 113 adjacent normal (0) samples.
```{r brca}
brca <- tcga[[2]]
brca
table(brca$GROUP)
```
## User-defined data compendium
With easy and fast access to the GEO2KEGG and TCGA compendia, enrichment
methods can be directly applied and assessed on well-studied, standardized
expression datasets.
Nevertheless, benchmarking with the `r Biocpkg("GSEABenchmarkeR")` package is
designed to be extendable to additional datasets as well.
Therefore, the `loadEData` function also accepts a directory where datasets,
preferably of class `r Biocpkg("SummarizedExperiment")`, have been saved as
[RDS](https://stat.ethz.ch/R-manual/R-devel/library/base/html/readRDS.html)
files.
```{r userComp}
data.dir <- system.file("extdata", package="GSEABenchmarkeR")
edat.dir <- file.path(data.dir, "myEData")
edat <- loadEData(edat.dir)
names(edat)
edat[[1]]
```
# Differential expression
To perform differential expression (DE) analysis between sample groups for
selected datasets of a compendium, the `r Biocpkg("GSEABenchmarkeR")` package
provides the function `runDE`.
The function invokes `EnrichmentBrowser::deAna` on each dataset, which contrasts
the sample groups as defined in the **GROUP** variable.
Here, we apply the function to 5 datasets of the GEO2KEGG compendium.
```{r deAna}
geo2kegg <- runDE(geo2kegg, de.method="limma", padj.method="flexible")
rowData(geo2kegg[[1]], use.names=TRUE)
```
Note: DE studies typically report a gene as differentially expressed
if the corresponding DE *p*-value, corrected for multiple testing, satisfies
the chosen significance level.
Enrichment methods that work directly on the list of DE genes are then
substantially influenced by the multiple testing correction.
An example is the frequently used over-representation analysis (ORA), which
assesses the overlap between the DE genes and a gene set under study based
on the hypergeometric distribution (see the vignette of the
`r Biocpkg("EnrichmentBrowser")` package, Appendix A, for an introduction).
ORA is inapplicable if there are few genes satisfying the significance threshold,
or if almost all genes are DE.
Using `padj.method="flexible"` accounts for these cases by applying multiple
testing correction in dependence on the observed degree of differential expression:
- the correction method from
[Benjamini and Hochberg](http://sci2s.ugr.es/keel/pdf/specific/articulo/Shaffer95MHT.pdf)
(BH) is applied, if it renders $\ge 1\%$ and $\le 25\%$ of all measured
genes as DE,
- the *p*-values are left unadjusted, if the BH correction results in $<1\%$
DE genes, and
- the more stringent
[Bonferroni](http://sci2s.ugr.es/keel/pdf/specific/articulo/Shaffer95MHT.pdf)
correction is applied, if the BH correction results in $>25\%$ DE genes.
Note that resulting $p$-values are not further used for assessing the
statistical significance of DE genes within or between datasets.
They are solely used to determine which genes are included in the analysis with
ORA - where the flexible correction ensures that the fraction of included genes
is roughly in the same order of magnitude across datasets.
Alternative strategies could also be applied (such as taking a constant number of
genes for each dataset or generally excluding ORA methods from the assessment).
# Enrichment analysis
In the following, we demonstrate how to carry out enrichment analysis in a
benchmark setup.
Therefore, we use the collection of human KEGG gene sets as obtained
with `getGenesets` from the `r Biocpkg("EnrichmentBrowser")`
package.
```{r getGS}
library(EnrichmentBrowser)
kegg.gs <- getGenesets(org="hsa", db="kegg")
```
At the core of applying a specific enrichment method to a single dataset is the
`runEA` function, which delegates execution of the chosen method to either
`EnrichmentBrowser::sbea` (set-based enrichment analysis) or
`EnrichmentBrowser::nbea` (network-based enrichment analysis).
In addition, it returns CPU time used and allows saving results for subsequent
assessment.
Here, we carry out ORA on the first dataset of the GEO2KEGG compendium.
```{r runEA}
kegg.ora.res <- runEA(geo2kegg[[1]], method="ora", gs=kegg.gs, perm=0)
kegg.ora.res
```
The function `runEA` can also be used to carry out several methods on multiple
datasets.
As an example, we carry out ORA and
[CAMERA](https://doi.org/10.1093/nar/gks461)
on 5 datasets of the GEO2KEGG compendium saving the results in a temporary
directory.
```{r eaAll}
res.dir <- tempdir()
res <- runEA(geo2kegg, methods=c("ora", "camera"),
gs=kegg.gs, perm=0, save2file=TRUE, out.dir=res.dir)
res$ora[1:2]
```
Note: saving the results to file is typically recommended when carrying out
several methods on multiple datasets for subsequent assessment.
This makes results, potentially obtained from time-consuming computations,
persistent across _R_ sessions.
In case of unexpected errors, this also allows resumption from the point of failure.
## User-defined enrichment methods
User-defined enrichment methods can easily be plugged into the benchmarking framework.
For demonstration, we define a dummy enrichment method that randomly draws _p_-values
from a uniform distribution.
```{r}
method <- function(se, gs)
{
ps <- runif(length(gs))
names(ps) <- names(gs)
return(ps)
}
```
We then execute this method on two datasets of the GEO2KEGG compendium using `runEA`
as before.
```{r}
res <- runEA(geo2kegg[1:2], method, kegg.gs)
res
```
# Benchmarking
Once methods have been applied to a chosen benchmark compendium, they can be
subjected to a comparative assessment of runtime, statistical significance,
and phenotype relevance.
To demonstrate how each criterion can be evaluated, we consider the example of
the previous section where we applied ORA and CAMERA on 5 datasets of the
GEO2KEGG compendium.
However, note that this minimal example is used to illustrate the basic
functionality in a time-saving manner - as generally intended in a vignette.
To draw conclusions on the individual performance of both methods, a more
comprehensive assessment, involving application to the full compendium, should be
carried out.
## Runtime
Runtime, i.e. CPU time used, is an important measure of the
applicability of a method.
For enrichment methods, runtime mainly depends on whether methods rely on
permutation testing, and how computationally intensive recomputation of
the respective statistic in each permutation is (see
[Figure 4](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4721010/figure/Fig4)
in
[Geistlinger et al., 2016](https://doi.org/10.1186/s12859-016-0884-1)).
To obtain the runtime from the application of ORA and CAMERA to 5 datasets of the
GEO2KEGG compendium, we can use the `readResults` function as we have saved
results to the indicated result directory in the above call of `runEA`.
```{r readRT}
ea.rtimes <- readResults(res.dir, names(geo2kegg),
methods=c("ora", "camera"), type="runtime")
ea.rtimes
```
For visualization of assessment results, the `bpPlot` function can be used to
create customized boxplots for specific benchmark criteria.
```{r plotRuntime, fig.width=6, fig.height=6}
bpPlot(ea.rtimes, what="runtime")
```
As both methods are simple gene set tests without permutation, they are among the
fastest in the field - with CAMERA being roughly twice as fast as ORA.
```{r runtimeORAvsCAMERA, fig.width=6, fig.height=6}
mean(ea.rtimes$ora) / mean(ea.rtimes$camera)
```
## Fraction of significant gene sets
The statistical accuracy of the significance estimation in gene set tests has
been repeatedly debated.
For example, systematic inflation of statistical significance in ORA due to an
unrealistic independence assumption between genes is well-documented
([Goeman and Buehlmann, 2007](https://doi.org/10.1093/bioinformatics/btm051)).
On the other hand, the permutation procedure incorporated in many gene set tests
has been shown to be biased
([Efron and Tibshirani, 2007](https://doi.org/10.1214/07-AOAS101)),
and also inaccurate if permutation $p$-values are reported as zero
([Phipson and Smyth, 2010](https://doi.org/10.2202/1544-6115.1585)).
These shortcomings can lead to inappropriately large fractions of significant
gene sets, and can considerably impair prioritization of gene sets in practice.
It is therefore important to evaluate resulting fractions of significant gene
sets in comparison to other methods and with respect to the fraction of
differentially expressed genes as a baseline.
We use the `readResults` function to obtain the saved gene set rankings of ORA
and CAMERA when applied to 5 datasets of the GEO2KEGG compendium (see above call
of `runEA`).
```{r readRankings}
ea.ranks <- readResults(res.dir, names(geo2kegg),
methods=c("ora", "camera"), type="ranking")
lengths(ea.ranks)
ea.ranks$ora[1:2]
```
The `evalNrSigSets` calculates the percentage of significant gene sets given a
significance level `alpha` and a multiple testing correction method `padj`.
We can visualize assessment results as before using `bpPlot`, which demonstrates
here that CAMERA produces substantially larger fractions of significant
gene sets than ORA.
```{r plotAdjSigSets, fig.width=6, fig.height=6}
sig.sets <- evalNrSigSets(ea.ranks, alpha=0.05, padj="BH")
sig.sets
bpPlot(sig.sets, what="sig.sets")
```
## Phenotype relevance
As introduced above, Tarca et al. ([2012](https://doi.org/10.1186/1471-2105-13-136)
and [2013](https://doi.org/10.1371/journal.pone.0079217)) also assigned a target
pathway to each dataset of the GEO2KEGG compendium, which is considered
highly-relevant for the respective phenotype investigated.
However, the relation between dataset, investigated phenotype, and assigned target
pathway is not always clear-cut.
In addition, there is typically more than one pathway that is considered relevant
for the investigated phenotype.
On the other hand, evaluations of published enrichment methods often conclude on
phenotype relevance, if there is *any* association between top-ranked gene sets
and the investigated phenotype.
A more systematic approach is used in the
[MalaCards](https://www.malacards.org) database of human diseases.
Here, relevance of [GO](http://www.geneontology.org)
and [KEGG](http://www.genome.jp/kegg) gene sets is summarized from
(*i*) experimental evidence and (*ii*) co-citation with the respective disease
in the literature.
### MalaCards disease relevance rankings
The `r Biocpkg("GSEABenchmarkeR")` package provides MalaCards relevance rankings
for the diseases investigated in the datasets of the GEO2KEGG and TCGA compendia.
Here, we load the relevance rankings for KEGG gene sets and demonstrate how they
can be incorporated in the assessment of phenotype relevance.
We note that the relevance rankings contain different numbers of gene sets for
different diseases, because only gene sets for which evidence/association with the
respective disease has been found are listed in a ranking.
For demonstration, we inspect the relevance rankings for Alzheimer's disease
(ALZ) and breast cancer (BRCA) containing 57 and 142 gene sets, respectively.
```{r malaRankings}
mala.kegg.file <- file.path(data.dir, "malacards", "KEGG.rds")
mala.kegg <- readRDS(mala.kegg.file)
sapply(mala.kegg, nrow)
mala.kegg$ALZ
mala.kegg$BRCA
```
### Mapping between dataset ID and disease code
To obtain the relevance ranking of the respective disease investigated when
assessing results on a specific dataset, a mapping between dataset and
investigated disease is required.
The function `readDataId2diseaseCodeMap` reads such a mapping from a tabular
text file and turns it into a named vector - where the elements correspond to
the disease codes and the names to the dataset IDs.
Here, we read the mapping between GSE ID and disease code for the GEO2KEGG
compendium.
```{r data2dis}
d2d.file <- file.path(data.dir, "malacards", "GseId2Disease.txt")
d2d.map <- readDataId2diseaseCodeMap(d2d.file)
head(d2d.map)
```
### Relevance score of a gene set ranking
To evaluate the phenotype relevance of a gene set ranking obtained from the
application of an enrichment method to an expression dataset, the function
`evalRelevance` assesses whether the ranking accumulates phenotype-relevant gene
sets (i.e. gene sets with high relevance scores) at the top.
Therefore, the function first transforms the ranks from the enrichment analysis
to weights - where the greater the weight of a gene set, the more it is ranked
towards the top of the GSEA ranking.
These weights are then multiplied by the corresponding relevance scores and
summed.
Here, we use `evalRelevance` to assess whether ORA, when applied to the GSE1297
dataset, recovers Alzheimer-relevant KEGG pathways.
```{r evalRelevance}
ea.ranks$ora$GSE1297
obs.score <- evalRelevance(ea.ranks$ora$GSE1297, mala.kegg$ALZ)
obs.score
```
### Random relevance score distribution
To assess the significance of the observed relevance score of an enrichment
method applied to a specific dataset, i.e. to assess how likely it is to
observe a relevance score equal or greater than the one obtained, the function
`compRand` repeatedly applies `evalRelevance` to randomly drawn gene set rankings.
For demonstration, we compute relevance scores for 50 random gene set rankings
and calculate the *p*-value as for a permutation test.
This demonstrates that the relevance score obtained from applying ORA to GSE1297
significantly exceeds random scores.
```{r compRand}
gs.names <- ea.ranks$ora$GSE1297$GENE.SET
gs.ids <- substring(gs.names, 1, 8)
rand.scores <- compRand(mala.kegg$ALZ, gs.ids, perm=50)
summary(rand.scores)
(sum(rand.scores >= obs.score) + 1) / 51
```
### Theoretical optimum
The observed relevance score can be used to compare phenotype relevance of
two or more methods when applied to one particular dataset.
However, as the number of gene sets in the relevance rankings differs between
phenotypes (see above Section *5.3.1 MalaCards disease relevance rankings*),
comparison between datasets is not straightforward as resulting relevance scores
scale differently.
Therefore, the function `compOpt` applies `evalRelevance` to the theoretically
optimal case in which the enrichment analysis ranking is identical to the
relevance score ranking. The ratio between observed and optimal score can then
be used to compare observed scores between datasets.
Here, we compute the optimal score for the Alzheimer relevance ranking, which
indicates that the score observed for ORA, when applied to GSE1297, is about 68%
of the optimal score.
```{r compOpt}
opt.score <- compOpt(mala.kegg$ALZ, gs.ids)
opt.score
round(obs.score / opt.score * 100, digits=2)
```
### Cross-dataset relevance score distribution
Evaluation of phenotype relevance with `evalRelevance` can also be done for
several methods applied across multiple datasets.
This allows to assess whether certain enrichment methods tend to produce
rankings of higher phenotype relevance than other methods when applied to a
compendium of datasets.
As explained in the previous section, observed relevance scores are always
expressed in relation to the respective optimal score.
For demonstration, we use `evalRelevance` to evaluate phenotype relevance of the
gene set rankings produced by ORA and CAMERA when applied to 5 datasets of the
GEO2KEGG compendium.
We can visualize assessment results as before using `bpPlot`, which demonstrates
here that ORA tends to recover more phenotype-relevant gene sets than CAMERA.
```{r evalAll, fig.width=6, fig.height=6}
all.kegg.res <- evalRelevance(ea.ranks, mala.kegg, d2d.map[names(geo2kegg)])
bpPlot(all.kegg.res, what="rel.sets")
```
### User-defined relevance rankings
It is also possible to refine the integrated MalaCards relevance rankings or to
incorporate relevance rankings for additional datasets.
For demonstration, we modify the KEGG relevance ranking for Alzheimer's disease
by providing a random relevance score for each gene set.
```{r}
rel.ranks <- mala.kegg$ALZ[,1:2]
rel.ranks$REL.SCORE <- runif(nrow(rel.ranks), min=1, max=100)
rel.ranks$REL.SCORE <- round(rel.ranks$REL.SCORE, digits = 2)
ind <- order(rel.ranks$REL.SCORE, decreasing = TRUE)
rel.ranks <- rel.ranks[ind,]
rel.ranks
```
We can then compute the aggregated relevance score of the ORA ranking according
to the updated relevance ranking using `evalRelevance` as before.
```{r}
evalRelevance(ea.ranks$ora$GSE1297, rel.ranks)
```
# Advanced
## Caching
Preparing an expression data compendium for benchmarking of enrichment methods
can be time-consuming.
In case of the GEO2KEGG compendium, it requires to summarize probe level
expression on gene level and to subsequently carry out differential expression
analysis for each dataset.
To flexibly save and restore an already processed expression data compendium,
we can use the `cacheResource` function which builds on functionality of the
`r Biocpkg("BiocFileCache")` package.
```{r cacheRes}
cacheResource(geo2kegg, rname="geo2kegg")
```
This adds the selected 5 datasets of the GEO2KEGG compendium (as processed
throughout this vignette) to the cache, and allows to restore it at a later time
via
```{r getRes}
geo2kegg <- loadEData("geo2kegg", cache=TRUE)
names(geo2kegg)
```
Note: to obtain the original unprocessed version of the compendium, set the
`cache` argument of the `loadEData` function to `FALSE`.
To clear the cache (use with care):
```{r clearCache, eval=FALSE}
cache.dir <- rappdirs::user_cache_dir("GSEABenchmarkeR")
bfc <- BiocFileCache::BiocFileCache(cache.dir)
BiocFileCache::removebfc(bfc)
```
## Parallel computation
Leveraging functionality from `r Biocpkg("BiocParallel")`, parallel computation
of the functions `maPreproc`, `runDE`, and
especially `runEA`, when applied to multiple datasets is straightforward.
Internally, these functions call `BiocParallel::bplapply`, which triggers parallel
computation as configured in the first element of `BiocParallel::registered()`.
As a result, parallel computation is implicitly incorporated in the above calls
of these functions when carried out on a multi-core machine.
See the vignette of the `r Biocpkg("BiocParallel")` package for an introduction.
Inspecting
```{r bpRegister}
BiocParallel::registered()
```
shows that the execution uses a `MulticoreParam` per default (on Windows:
a `SnowParam`), where the `bpnworkers` attribute indicates the number of cores
involved in the computation.
To change the execution mode of functions provided in the
`r Biocpkg("GSEABenchmarkeR")` package, accordingly configured computation
parameters of class `BiocParallelParam` can either directly be registered via
`BiocParallel::register`, or supplied with the `parallel` argument of the
respective function.
For demonstration, we configure here a `BiocParallelParam` to display a progress
bar
```{r bpParam}
bp.par <- BiocParallel::registered()[[1]]
BiocParallel::bpprogressbar(bp.par) <- TRUE
```
and supply `runDE` with the updated computation parameter.
```{r runDEBP}
geo2kegg <- runDE(geo2kegg, parallel=bp.par)
```
Users that would like to use distributed computation, on e.g. an institutional
computer cluster, should consult the vignette of the `r Biocpkg("BiocParallel")`
package to similarly configure a `BiocParallelParam` for that purpose.