In lianos/multiGSEA: A unified interface to a plethora of gene set enrichment analysis methods

knitr::opts_chunk$set(
  echo=TRUE, warning=FALSE, message=FALSE, error=FALSE) #, dpi=150)

# TODO: 
# 1. update text to favor overrepresentation analysis method `ora`
# 2. describe how its data.frame input works (ora)
# 3. provide a brief comparison to its performance vs goseq, with and without 
#    accounting for bias). You can take the code for that from the 
#    enrichtest/goseq PAC unit test.

Overview

The multiGSEA 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 multiGSEA 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 MultiGSEAResult object which stores the results of each GSEA method internally, allowing for easy query and retrieval.
• A sister multiGSEA.shiny package provides an explore function, which is invoked on MultiGSEAREsult objects returned from a call to multiGSEA. 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 MultiGSEAResult *.rds file.
• An "over representation analysis" method ora() which wraps the biased enrichment functionality found within limma::kegga and generalizes it to work against data.frame inputs with arbitrary genesets.
• The 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.
• The GeneSetDb class: a new class to store collections of gene sets. This provides a bit more functionality than the base GeneSetCollection class defined in the GSEABase.
• A convenience getMSigGeneSetDb() function, which enables the user to fetch arbitraty subsets of the MSigDB collections mapped to different organisms and identifier types.

The [BiocSet][] package was introduced in Bioconductor version 3.10, which should supersede the `GeneSetDb` class introduced here. multiGSEA will endeavor to transition to using that class to store gene set collections in the near future.

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

• Expression data (an EList, DGEList, or expression matrix)
• A design matrix
• A contrast vector/matrix (if your design and comparison require it)

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 multiGSEA call that can orchestrate multiple analyses of any type.

All GSEA methods require the use of a GeneSetDb object, which is a new class provided by this package that holds geneset collection information which is stored primarily as tables / data.frames.

Currently supported gene set enrichment methods include:

dplyr::select(multiGSEA::multiGSEA_methods(),
              method, test_type, package, comment)

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

The multiGSEA 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.

Standard Workflow

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

library(multiGSEA)
library(magrittr)
library(reshape2)
library(dplyr)
library(ggplot2)
library(ComplexHeatmap)
library(circlize)
library(edgeR)

theme_set(theme_bw())

Internally, mulitGSEA leverages the [data.table](https://CRAN.R-project.org/package=data.table) package for fast indexing and manipulation over data.frames. All functions that return these objects back to the user have an `as.dt` argument, which is set to `FALSE` by default. This means that unless the end user wants to work with `data.table` objects, they will always be returned a `data.frame`.

Data Setup

multiGSEA 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 Entrez IDs for gene gene identifiers.

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)

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](http://bioconductor.org/packages/release/data/experiment/vignettes/airway/inst/doc/airway.html) for full details (among others).

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:

vm$design

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

(cm <- makeContrasts(BvH=Basal - Her2, levels=vm$design))

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) %>%
  contrasts.fit(cm) %>%
  eBayes
tt <- topTable(fit, 'BvH', n=Inf, sort.by='none')

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 multiGSEA 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.

Gene Sets to Test

The multiGSEA package provides a GeneSetDb class, which is used instead of a GSEABase::GeneSetCollection to house collections of genesets. You can easily create a GeneSetDb object from a custom set of gene sets you have on hand, but the package also provides convenience wrappers to generate GeneSetDb objects from popular gene annotation resources such as MSigDB, PANTHER, etc.

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

gdb <- getMSigGeneSetDb(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 %>%
  select(1:4)

For more details on creating and manipulating GeneSetDb objects, please jump the to The GeneSetDb Class section.

Running multiGSEA

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

The call to multiGSEA will perform these analyses and return a MultiGSEAResult object which you can then use for downstream analysis.

mg <- multiGSEA(
  gdb, vm, vm$design, cm[, 'BvH'],
  methods = c('camera', 'fry', 'ora'),
  # 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 multiGSEA call shortly ...

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 multiGSEA, a limma-based lmFit %>% (eBayes|treat) %>% (topTable|topTreat) pipeline was run. If a DGEList was passed, then multiGSEA 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 MultiGSEAResult object:

lfc <- logFC(mg)
lfc %>%
  select(symbol, entrez_id, logFC, t, pval, padj) %>%
  head

We can confirm that the statistics generated internally in multiGSEA 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)

The internally performed differential expression analysis within the mulitGSEA 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 multiGSEA'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 multiGSEA function) for more information.

Explicit GSEA

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

mg

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 MultiGSEAResult object simply calls the tabulateResults function, which you can call directly with the value of max.p that you might find more appropriate.

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 MultiGSEAResult object:

resultNames(mg)

Note that when running a "goseq" analysis, multiGSEA will (by default) run it three different ways. By running an enrichment analysis on all differentially expressed genes, 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) %>% 
  head

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, 'H', 'HALLMARK_WNT_BETA_CATENIN_SIGNALING') %>% 
  select(symbol, entrez_id, logFC, pval, padj) %>%
  head

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.

Plotting

multiGSEA 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.

Boxplot

iplot(mg, 'H', 'HALLMARK_WNT_BETA_CATENIN_SIGNALING',
      type='boxplot', value='logFC')

Density

iplot(mg, 'H', 'HALLMARK_WNT_BETA_CATENIN_SIGNALING',
      type='density', value='logFC')

Interactive Exploration

A sister multiGSEA.shiny package is available that can be used to interactively explore MultiGSEAResult objects. The application can be invoked as follows:

library("multiGSEA.shiny")
explore(mg)

Please refer to the "multiGSEA-shiny"" vignette in the multiGSEA.shiny package for documentation on the application's use.

Singe Sample Gene Set Scoring

It can be both convenient and effective to transform a gene-by-sample expression matrix to a geneset-by-sample expression matrix. By doing so, so we can quickly identify biological processes that are up/down regulated (loosely speaking) in each sample.

We can generate single sample gene set scores using the gene sets defined in a GeneSetDb using the scoreSingleSamples function. This function takes a GeneSetDb, an expression container, and a methods argument, which is analagous to the methods argument in the multiGSEA call: it defines all of the scoring methos the user wants to apply to each sample.

Let's pick a few gene sets to score our samples with for this exercise. We'll take the significant hallmark gene sets, or any other significant gene set that has a large (on average) log fold change beteen conditions.

sig.res <- cam.res %>% 
  filter(padj < 0.05 & (collection == 'H' | abs(mean.logFC) >= 2))
gdb.sub <- gdb[geneSets(gdb)$name %in% sig.res$name]

Refer to the [Subsetting a GeneSetDb](#subsetting-a-genesetdb) section to learn how to subset a `GeneSetDb` object to create a derivative object with fewer gene sets.

Recall that the GSEA analysis we performed was perfomed between the Basal and Her2 subtypes, so we will use an expression matrix that only has the samples from those two groups.

vm.bh <- vm[, vm$targets$PAM50subtype %in% c("Basal", "Her2")]

Generating Single Sample Gene Set Scores

Once we have a GeneSetDb object that contains all of the gene sets we wish to use to create single sample gene set scores, we can use the scoreSingleSamples function to produce these scores using a variety of algorithmes, which the user species using the methods parameter.

The scoreSingleSamples function will return a long data.frame with length(methods) * ncol(exprs) rows. Each row represents the score for the given sample using the specified method. You can subset against the method column to extract all of the single sample scores for a given method.

scores <- scoreSingleSamples(gdb.sub, vm.bh,
                             methods=c('ewm', 'ssgsea', 'zscore'),
                             ssgsea.norm=TRUE,
                             unscale=FALSE, uncenter=FALSE)

We can see how the scores from different methods compare to each other with using a little reshape2 mojo and the multiGSEA::corplot function.

sm <- acast(scores, name + sample_id ~ method, value.var="score")
corplot(sm, cluster=TRUE)

It is, perhaps, interesting to compare how the ewm method scores change when we choose not to "uncenter" and "unscale" them:

ewmu <- scoreSingleSamples(gdb.sub, vm.bh,
                           methods=c('ewm'),
                           unscale=TRUE, uncenter=TRUE) %>% 
  mutate(method='ewm_unscale')
scores.all <- bind_rows(scores, ewmu)
sma <- acast(scores.all, name + sample_id ~ method, value.var="score")
corplot(sma, cluster=TRUE)

Further exposition on the "ewm" (eigenWeightedMean) scoring method can be found in the ?eigenWeightedMean function.

Similarity of ewm:ssGSEA and ewz/zscore to gsva

TODO: explore this more deeply, but the examples in ?scoreSingleSamples suggests that to be the case in that example.

Visualizing Single Sample Gene Set Scores

The "long" data.frame nature of the results produced by scoreSingleSamples makes it convenient to use with graphing libraries like ggplot2 so that we can create arbitrary visualizations. Creating boxplots for gene sets per subtype is an easy way to explore these results.

Let's annotate each row in scores.all with the subtype annotation and observe how these methods score each sample for a few gene sets.

all.scores <- scores.all %>% 
  inner_join(select(vm.bh$targets, sample_id=Sample_ID, subtype=PAM50subtype),
             by = "sample_id")

some.scores <- all.scores %>% 
  filter(name %in% head(unique(all.scores$name), 5))

ggplot(some.scores, aes(subtype, score)) +
  geom_boxplot(outlier.shape=NA) +
  geom_jitter(width=0.25) +
  facet_grid(name ~ method)

Gene Set Based Heatmap with mgheatmap

We often want to create expression based heatmaps that highlight the behavior of gene sets across our samples. The mgheatmap function uses the ComplexHeatmap package to create two different types of heatmaps:

1. Gene based heatmaps, that split the genes (rows) based on their genesets
2. Single sample gene set based heatmaps, optionally split by gene set collection.

The mgheatmap function has a set of arguments that customize how the heatmap is to be created (gene level vs. gene set level, whether to split it, etcv) and will also use the ... argument to pass any parameters down to the inner ComplexHeatmap::Heatmap function call and customize its behavior. The mgheatmap function returns a ComplexHeatmap,Heatmap object for plotting or combining with other ComplexHeatmap heatmaps or annotations in order to create arbitrarily complex/informative heatmap figures.

Gene level based heatmap (from genesets)

You can plot a heatmap of the genes from a predefined set of gene sets by providing the gene sets you want to visualize in a GeneSetDb object.

We'll create a new GeneSetDb object using the first two gene sets in gdb.sub and draw a heatmap of their expression.

gs.sub <- geneSets(gdb.sub)
gdb.2 <- gdb.sub[geneSets(gdb.sub)$name %in% head(gs.sub$name, 2)]

col.anno <- HeatmapAnnotation(
  df = vm.bh$targets[, 'PAM50subtype', drop = FALSE],
  col = list(PAM50subtype = c(Basal = "gray", Her2 = "black")))

mgheatmap(vm.bh, gdb.2, aggregate.by = "none", split = TRUE,
          show_row_names = FALSE, show_column_names = FALSE,
          recenter = TRUE, top_annotation = col.anno, zlim = c(-3, 3))

Gene set-based heatmap

You can often get a higher information:ink ratio by plotting heatmaps based on single sample gene set scores as opposed to the genes that make up a geneset.

Let's see what the simple 2-geneset version of the heatmap above looks like:

mgheatmap(vm.bh, gdb.2, aggregate.by = "ewm", split = FALSE,
          show_row_names = TRUE, show_column_names = FALSE,
          top_annotation = col.anno)

Plotted in this way, we can now show the activity of a greater number of genesets

mgheatmap(vm.bh, gdb.sub, 
          aggregate.by = 'ewm', split=TRUE, recenter = TRUE,
          show_row_names=TRUE, show_column_names=FALSE,
          top_annotation=col.anno, zlim = c(-2.5, 2.5))

The GeneSetDb Class

The GeneSetDb class is a new container to store collections of genesets which provides different types of functionality than is found in GSEABase::GeneSetCollection objects.

We can, for instance, identify all 84 gene sets that have genes "10014" and "1454" (entrez ids) as members (HDAC5 and CSNK1E, respectively).

gdb.sub <- subsetByFeatures(gdb, c('10014', '1454'))
nrow(gdb); nrow(gdb.sub)

The GeneSetDb object uses the data.table package internally for fast lookup. The code will be optimized in the future to be even more performant. Internally the collection of gene set information is minimally stored as a three-column data.table in "long form", which has the following columns:

• collection
• name
• feature_id

More columns can be added to the internal data.table (a "symbol" column, for instance), but those are the only three you need. To see what we are talking about, exactly, you can call the as.data.frame function on a GeneSetDb object:

as.data.frame(gdb)[c(1:5, 201:205),]

The (collection,name) tuple is the primary key of a gene set. The feature_id column stores gene identifiers. For the time being, it will be most natural for these IDs to simply be ensembl gene identifiers (or entrez ids) as many of the annotation databases use these identifiers, as well.

Building a GeneSetDb

The multiGSEA package provides convenience funcitons to fetch genesets from many sources and convert them into a GeneSetDb object. The two most useful sources may be:

• MSigDB via getMSigGeneSetDb(...). Although the core multiGSEA package provides the getter function for these genesets, the user needs to install the msigdb.data package to include the data and utility functions for the identifier and species mapping.
• PANTHER (pathways and GOSLIM) via getPantherGeneSetDb()

You can create a custom GeneSetDb via the GeneSetDb constructor, which accpets the following inputs:

1. A data.frame of geneset membership. This requires collection, name, and feature_id columns. Reference the output of as.data.frame(gdb) shown above.
2. A named list of gene identifier vectors that represent genesets for a single collection
3. A named list of (2)-like lists. The top level names are the names of the different collections, and each sublist represents the genesets in that collection.

Two GeneSetDb objects can be combined using the append function. For now it is your responsibility to ensure that the two GeneSetDb objedts are "reasonably conformable", ie. they use the same types of gene identifiers, and are referencing the same species, etc.

msigdb <- getMSigGeneSetDb('H', 'human')
goslimdb <- getPantherGeneSetDb('goslim', 'human')
gdb.uber <- combine(msigdb, goslimdb)

See the help and examples in ?GeneSetDb for more information.

For some reason the `PANTHER.db` package needs to be installed in a user-writable package location for this to work properly. If you see an error that speaks to using "rsqlite to write to a readonly database", you will have to install `PANTHER.db` in a user-writable directory using `BiocInstaller::biocLite("PANTHER.db")`

Subsetting a GeneSetDb

The subsetting functionality for a GeneSetDb isn't quite as fluent as I would like to be. Improvements will be provided in a future release.

Let's say we wanted to include only the genesets where their "subcategory" was "[CP:PID][msigdbpid]", we would like to do it like this:

gdb.sub <- subset(gdb, subcategory == "CP:PID")

Unfortunatey we aren't there yet. Currently have to first define a logical vector over the data.frame you get from geneSets(gdb). Then use that vector to remove gene sets using the GeneSetDb,"[" function:

keep <- geneSets(gdb)$subcategory == "CP:PID"
gdb.sub <- gdb[keep]

To reiterate, length(keep) must equal nrow(geneSets(gdb)).

We have also showed you above how to create a subsetted GeneSetDb by keeping only the gene sets that have certain features in them using the subsetByFeatures function.

Packaged MSigDb Gene Sets

TODO: This needs update to reflect our use of msigdbr within the msigdb.data package.

The multiGSEA package provides the getMSigGeneSetDb utility function to provide easy access to the MSigDB gene set collection for users. The data for this function is stored in data package that need to be downloaded first before use. The GeneSetDb.MSigDB.Hsapiens.v*, and GeneSetDb.MSigDB.Mmusculus.v* packages provide these genesets for human and mouse (using entrez identifiers). The latest version of these packages (v61) provides version 6.1 of the MSigDB gene sets, which is the latest version as of this relase.

The MSigDB GeneSetDb objects include an organism column in their geneSets(gdb) table. Note that this column indicates which organism that the geneset was defined in. The identifiers in the genesets provided in the GeneSetDb are all either "human" or "mouse" identifiers, depending on the value of the species argument in your getMSigGeneSetDb call. The example above shows how to subset the GeneSetDb based on the organism column so that the resulting GeneSetDb only has genesets defined in mouse. The user can also set the species.specific parameter to TRUE so that we only return gene sets defined in the organism ('human' or 'mouse') that the user is retrieving the genesets from. For instance, the call below only retrieves gene sets in the "c2" and "c5" collections that were only defined in mouse.

# This no longer is included
mgdb <- getMSigGeneSetDb(c('C2' 'C7'), species='mouse',
                         species.specific=TRUE)

Active vs Inactive Gene Sets

A GeneSetDb is used to hold "the universe" of genes that belong to different gene sets across different collections. Depending on the assay performed to measure these genes, the set of genes you observe in your study will likely be a subset of the genes in the GeneSetDb. As such, prior to using a GeneSetDb for GSEA, it must be "conformed" to a target object that will be used for the input to the GESA (either a matrix of expression, or a pre ranked vector of statistics). This step will index into the target expression object and identify which rows of the object correspond to which genes in the GeneSetDb.

"Conformation" happens automatically within the multiGSEA call, but we call it explicitly below to outline its functionality. The command below conforms the GeneSetDb to our target "voomed" EList, and deactivates gene sets (i.e. removes them from downstream GSEA) that have less than 10 or more than 100 genes that were found in vm:

gdbc <- conform(gdb, vm, min.gs.size=10, max.gs.size=100)
head(geneSets(gdbc, active.only=FALSE))

We can see that, only 23 of the 26 genes in the (C2,ABBUD_LIF_SIGNALING_1_DN) were found in the rows of vm, and the (C2,ABBUD_LIF_SIGNALING_2_DN) was "deactivated." Deactivated (active == FALSE) gene sets will be ignored during downstream analyses. This gene set was deactivated because it only has five "conformed" genes, but the minimum geneset size we wanted to consider (min.gs.size) was set to ten in our call to conform.

Accessing members of a gene set

The geneSet and featureIds functions allow the user to identify the genes found in a gene set. Both of these functions take an active.only argument, which is TRUE by default. This specifies that only the genes that have been successfully conformed to a gene set should be the ones that are returned.

For instance, we can identify which genes belong to the (C2,ABBUD_LIF_SIGNALING_1_DN), and which three were not found in vm like so:

missed <- setdiff(
  featureIds(gdbc, 'C2', 'ABBUD_LIF_SIGNALING_1_DN', active.only=FALSE),
  featureIds(gdbc, 'C2', 'ABBUD_LIF_SIGNALING_1_DN', active.only=TRUE))
missed

or we can use the geneSet function to return a data.frame of these results:

gdbc %>% 
  geneSet('C2', 'ABBUD_LIF_SIGNALING_1_DN', active.only = FALSE) %>%
  subset(feature_id %in% missed)

Mapping of gene set featureIds to target expression containers

It may be that the IDs used in a gene set collection are different from the ones used as the rownames of your expression container. For instance, the IDs used for a given gene set collection in the GeneSetDb might be Ensembl gene identifiers, but the rownames of the expression object might be Entrez ID. This is where the mapping parameter becomes useful.

The GeneSetDb class has a concept of an internal featureIdMap to accommodate these scenarios, which would allow for a non-destructive mapping of the original IDs to a new "ID space" (entrez to ensembl, for instance).

This functionality is not ready for this release, but it's just a note to keep the user aware of some future development of the package. For the time being, the user is required to manually map the feautreIds in their expression matrix to be concordant with the ones found in the GeneSetDb.

In the meantime, a rename_rows convenience function is provided here to easily rename the rows of our expression container to different values. For instance, to rename this is how you might rename the rows of your assay container to use symbols:

vm <- exampleExpressionSet()
vms <- rename_rows(vm, "symbol")
head(cbind(rownames(vm), rownames(vms)))

We grabbed the symbol column from vm$genes and "smartly" renamed the rows of vm with the values there. Refer to the ?rename_rows man page for more details. This, of course, still requires you to manually fetch and map identifiers, but still ...

Customizing Analyses

The internal differential expression analysis as well the gene set enrichment analyses can be customized by passing parameters through the ... in the multiGSEA function.

Custom Differential Expression

The internal differential expression pipeline, exported via the calculateIndividualLogFC function allows the end user to configure an "arbitrarily complex" differential expression analysis using either edgeR's quasilikelihood framework (if the input is a DGEList) or a direct limma analysis (with a pre-voomed EList, expression matrix, or whatever).

User's should refer to the ?calculateIndividualLogFC help page to see which parameters are exposed for a differential expression analysis and configure them accordingly. When calling multiGSEA use these same parameters in the call and they will be provided to calculateIndividualLogFC.

For instance, if you wanted to use limma's "treat" functionality to specify a minimal log fold change threshold for statistical significance, you would do so as follows:

mg <- multiGSEA(gdb, vm, vm$design, cm[, 'BvH'],
                methods=c('goseq'),
                treat.lfc=log2(1.5),
                ## feature length vector required for goseq
                feature.bias=setNames(vm$genes$size, rownames(vm)))

Using the internal treat functionality would really only affect enrichment tests that first threshold the genes in your experiment as "significant" or not, like goseq and not tests like camera.

Custom GSEA

The GSEA methods that are wrapped by multiGSEA all take the same parameters that are defined by their implementation. Simply pass these parameters down via the ... in the multiGSEA call.

For instance, you can read ?camera to find that the camera method accepts an inter.gene.cor parameter, and ?roast will tell you that you can specify the number of rotations used via the nrot parameter.

mgx <- multiGSEA(gdb, vm, vm$design, cm[, 'BvH'],
                 methods=c('camera', 'roast'),
                 inter.gene.cor=0.04, nrot=500)

Developing multiGSEA

Adding new GSEA methods

This section needs more explanation

Suppose we wanted to add a new GSEA method named superGSEA to the methods that the multiGSEA can delegate to via its methods argument, we need to add the following internal multiGSEA methods.

• validate.inputs.superGSEA: Choose either .validate.inputs.full.design or .validate.inputs.preranked depending on whether or not the method needs:
  • a full expression,design,contrast tripple like roast, romer, or camera, use .validate.inputs.full.design;
  • or a preranked set of stats vector like cameraPR or fgsea, use .validate.inputs.preranked
• validate.x.superGSEA: this is redundant, just put validate.X. This will be fixed in a future release
• do.superGSEA: this method should take the parameters listed below, and return the result of the superGSEA call unmodified from its original form that the superGSEA method that is being wrapped returns it. The parameters of the do.superGSEA function are:
  • gsd: the pre-conformed GeneSetDb
  • x: the expression matrix or pre-ranked stats vector
  • design: the design matrix. If the method uses a pre-ranked stats vector just ignore this argument in the do.superGSEA function body.
  • contrast: the contrast to test. If the method uses a pre-ranked stats vector just ignore this argument in the do.superGSEA function body.
  • gsd.idxs: this will be a list of gene sets. The names are the collection;;name tuples pasted together, and the values are integers indices into the rows of x for each gene in the gene set. Genes in the gene set that are not in x have already been removed.
  • any custom parameters you want to set the defaults for in superGSEA when the default values are not what you want.
  • ...: any other formal arguments defined in superGSEA will be passed into here, and it will be your responsibility to extract them and pass them down into the superGSEA call.
• mgres.superGSEA: The function takes the output from do.superGSEA and turns it into a data.table that minimally has collection, name, pval, and padj columns.

Look to the implementation in the do.camera.R file for a reference.

lianos/multiGSEA documentation built on Nov. 17, 2020, 1:26 p.m.