2.1 A primer on terminology, existing methods & statistical theory

Gene sets, pathways & regulatory networks

Gene sets are simple lists of usually functionally related genes without further specification of relationships between genes.

Pathways can be interpreted as specific gene sets, typically representing a group of genes that work together in a biological process. Pathways are commonly divided in metabolic and signaling pathways. Metabolic pathways such as glycolysis represent biochemical substrate conversions by specific enzymes. Signaling pathways such as the MAPK signaling pathway describe signal transduction cascades from receptor proteins to transcription factors, resulting in activation or inhibition of specific target genes.

Gene regulatory networks describe the interplay and effects of regulatory factors (such as transcription factors and microRNAs) on the expression of their target genes.

Resources

GO and KEGG annotations are most frequently used for the enrichment analysis of functional gene sets. Despite an increasing number of gene set and pathway databases, they are typically the first choice due to their long-standing curation and availability for a wide range of species.

GO: The Gene Ontology (GO) consists of three major sub-ontologies that classify gene products according to molecular function (MF), biological process (BP) and cellular component (CC). Each ontology consists of GO terms that define MFs, BPs or CCs to which specific genes are annotated. The terms are organized in a directed acyclic graph, where edges between the terms represent relationships of different types. They relate the terms according to a parent-child scheme, i.e. parent terms denote more general entities, whereas child terms represent more specific entities.

KEGG: The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a collection of manually drawn pathway maps representing molecular interaction and reaction networks. These pathways cover a wide range of biochemical processes that can be divided in 7 broad categories: metabolism, genetic and environmental information processing, cellular processes, organismal systems, human diseases, and drug development. Metabolism and drug development pathways differ from pathways of the other 5 categories by illustrating reactions between chemical compounds. Pathways of the other 5 categories illustrate molecular interactions between genes and gene products.

Gene set analysis vs. gene set enrichment analysis

The two predominantly used enrichment methods are:

• Overrepresentation analysis (ORA), testing whether a gene set contains disproportional many genes of significant expression change, based on the procedure outlined in the first section
• Gene set enrichment analysis (GSEA), testing whether genes of a gene set accumulate at the top or bottom of the full gene vector ordered by direction and magnitude of expression change Subramanian et al., 2005

However, the term gene set enrichment analysis nowadays subsumes a general strategy implemented by a wide range of methods Huang et al., 2009. Those methods have in common the same goal, although approach and statistical model can vary substantially Goeman and Buehlmann, 2007, Khatri et al., 2012.

To better distinguish from the specific method, some authors use the term gene set analysis to denote the general strategy. However, there is also a specific method from Efron and Tibshirani, 2007 of this name.

Underlying null: competitive vs. self-contained

Goeman and Buehlmann, 2007 classified existing enrichment methods into competitive and self-contained based on the underlying null hypothesis.

• Competitive null hypothesis: the genes in the set of interest are at most as often DE as the genes not in the set,
• Self-contained null hypothesis: no genes in the set of interest are DE.

Although the authors argue that a self-contained null is closer to the actual question of interest, the vast majority of enrichment methods is competitive.

Goeman and Buehlmann further raise several critical issues concerning the 2x2 ORA:

• rather arbitrary classification of genes in DE / not DE
• based on gene sampling, although sampling of subjects is appropriate
• unrealistic independence assumption between genes, resulting in highly anti-conservative p-values

With regard to these statistical concerns, GSEA is considered superior:

• takes all measured genes into account
• subject sampling via permutation of class labels
  
• the incorporated permutation procedure implicitly accounts for correlations between genes

However, the simplicity and general applicability of ORA is unmet by subsequent methods improving on these issues. For instance, GSEA requires the expression data as input, which is not available for gene lists derived from other experiment types. On the other hand, the involved sample permutation procedure has been proven inaccurate and time-consuming Efron and Tibshirani, 2007, Phipson and Smyth, 2010, Larson and Owen, 2015.

Generations: ora, fcs & topology-based

Khatri et al., 2012 have taken a slightly different approach by classifying methods along the timeline of development into three generations:

1. Generation: ORA methods based on the 2x2 contingency table test,
2. Generation: functional class scoring (FCS) methods such as GSEA, which compute gene set (= functional class) scores by summarizing per-gene DE statistics,
3. Generation: topology-based methods, explicitly taking into account interactions between genes as defined in signaling pathways and gene regulatory networks (Geistlinger et al., 2011 for an example).

Although topology-based (also: network-based) methods appear to be most realistic, their straightforward application can be impaired by features that are not-detectable on the transcriptional level (such as protein-protein interactions) and insufficient network knowledge Geistlinger et al., 2013, Bayerlova et al., 2015.

Given the individual benefits and limitations of existing methods, cautious interpretation of results is required to derive valid conclusions. Whereas no single method is best suited for all application scenarios, applying multiple methods can be beneficial. This has been shown to filter out spurious hits of individual methods, thereby reducing the outcome to gene sets accumulating evidence from different methods Geistlinger et al., 2016, Alhamdoosh et al., 2017.

2.2 Data types

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) made adaptations in differential expression analysis and, to some extent, also in gene set enrichment analysis necessary.

Thus, we consider two example datasets - a microarray and a RNA-seq dataset, and discuss similarities and differences of the respective analysis steps.

For microarray data, we consider expression measurements of patients with acute lymphoblastic leukemia Chiaretti et al., 2004. A frequent chromosomal defect found among these patients is a translocation, in which parts of chromosome 9 and 22 swap places. This results in the oncogenic fusion gene BCR/ABL created by positioning the ABL1 gene on chromosome 9 to a part of the BCR gene on chromosome 22.

We load the ALL dataset

library(ALL)
data(ALL)

and select B-cell ALL patients with and without the BCR/ABL fusion, as described previously Gentleman et al., 2005.

ind.bs <- grep("^B", ALL$BT)
ind.mut <- which(ALL$mol.biol %in% c("BCR/ABL", "NEG"))
sset <- intersect(ind.bs, ind.mut)
all.eset <- ALL[, sset]

We can now access the expression values, which are intensity measurements on a log-scale for 12,625 probes (rows) across 79 patients (columns).

dim(all.eset)

## Features  Samples 
##    12625       79

exprs(all.eset)[1:4,1:4]

##              01005    01010    03002    04007
## 1000_at   7.597323 7.479445 7.567593 7.905312
## 1001_at   5.046194 4.932537 4.799294 4.844565
## 1002_f_at 3.900466 4.208155 3.886169 3.416923
## 1003_s_at 5.903856 6.169024 5.860459 5.687997

As we often have more than one probe per gene, we compute gene expression values as the average of the corresponding probe values.

all.eset <- probe.2.gene.eset(all.eset) 
head(featureNames(all.eset))

## [1] "5595" "7075" "1557" "643"  "1843" "4319"

For RNA-seq data, we consider transcriptome profiles of four primary human airway smooth muscle cell lines in two conditions: control and treatment with dexamethasone Himes et al., 2014.

We load the airway dataset

library(airway)
data(airway)

and create an Biobase::ExpressionSet.

air.eset <- as(airway, "ExpressionSet")
annotation(air.eset) <- "hsa"

Note: we are switching here from SummarizedExperiment to ExpressionSet as functionality in the EnrichmentBrowser is explicitly designed for gene-level expression data, for which ExpressionSet can be consistently used.

For further analysis, we only keep genes that are annotated to an ENSEMBL gene ID.

air.eset <- air.eset[grep("^ENSG", rownames(air.eset)), ]
dim(air.eset)

## Features  Samples 
##    63677        8

exprs(air.eset)[1:4,1:4]

##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513
## ENSG00000000003        679        448        873        408
## ENSG00000000005          0          0          0          0
## ENSG00000000419        467        515        621        365
## ENSG00000000457        260        211        263        164

2.3 Differential expression analysis

Normalization of high-throughput expression data is essential to make results within and between experiments comparable. Microarray (intensity measurements) and RNA-seq (read counts) data typically show distinct features that need to be normalized for. As this is beyond the scope of this workshop, we refer to limma for microarray normalization and EDASeq for RNA-seq normalization. See also EnrichmentBrowser::normalize, which wraps commonly used functionality for normalization.

The EnrichmentBrowser incorporates established functionality from the limma package for differential expression analysis. This involves the voom transformation when applied to RNA-seq data. Alternatively, differential expression analysis for RNA-seq data can also be carried out based on the negative binomial distribution with edgeR and DESeq2.

This can be performed using the function EnrichmentBrowser::de.ana and assumes some standardized variable names:

• GROUP defines the sample groups being contrasted,
• BLOCK defines paired samples or sample blocks, as e.g. for batch effects.

For more information on experimental design, see the limma user’s guide, chapter 9.

For the ALL dataset, the GROUP variable indicates whether the BCR-ABL gene fusion is present (1) or not (0).

pData(all.eset)$GROUP <- ifelse(all.eset$mol.biol == "BCR/ABL", 1, 0)
table(pData(all.eset)$GROUP)

## 
##  0  1 
## 42 37

For the airway dataset, it indicates whether the cell lines have been treated with dexamethasone (1) or not (0).

pData(air.eset)$GROUP <- ifelse(colData(airway)$dex == "trt", 1, 0)
table(pData(air.eset)$GROUP)

## 
## 0 1 
## 4 4

Paired samples, or in general sample batches/blocks, can be defined via a BLOCK column in the pData slot. For the airway dataset, the sample blocks correspond to the four different cell lines.

pData(air.eset)$BLOCK <- colData(airway)$cell
table(pData(air.eset)$BLOCK)

## 
## N052611 N061011 N080611  N61311 
##       2       2       2       2

For microarray data, the EnrichmentBrowser::de.ana function carries out differential expression analysis based on functionality from the limma package. Resulting log2 fold changes and t-test derived p-values for each gene are appended to the fData slot.

all.eset <- de.ana(all.eset)
head(fData(all.eset), n=4)

##               FC  ADJ.PVAL limma.STAT
## 5595  0.04296986 0.8996597  0.7346897
## 7075  0.03208350 0.9494365  0.4546963
## 1557 -0.04394014 0.8184362 -1.0658081
## 643  -0.02775435 0.9296705 -0.5674048

Nominal p-values are already corrected for multiple testing (ADJ.PVAL) using the method from Benjamini and Hochberg implemented in stats::p.adjust.

For RNA-seq data, the de.ana function can be used to carry out differential expression analysis between the two groups either based on functionality from limma (that includes the voom transformation), or alternatively, the frequently used edgeR or DESeq2 package. Here, we use the analysis based on edgeR.

air.eset <- de.ana(air.eset, de.method="edgeR")

## Excluding 50740 genes not satisfying min.cpm threshold

head(fData(air.eset), n=4)

##                          FC     ADJ.PVAL  edgeR.STAT
## ENSG00000000003 -0.40494099 9.700883e-05 18.75561361
## ENSG00000000419  0.18305320 1.356563e-01  3.48605567
## ENSG00000000457  0.01455321 9.423545e-01  0.01477962
## ENSG00000000460 -0.13467114 6.537338e-01  0.43278224

2.4 Gene sets

We are now interested in whether pre-defined sets of genes that are known to work together, e.g. as defined in the Gene Ontology or the KEGG pathway annotation, are coordinately differentially expressed. The function get.kegg.genesets downloads all KEGG pathways for a chosen organism (here: Homo sapiens) as gene sets using NCBI Entrez Gene IDs.

kegg.gs <- get.kegg.genesets("hsa")

Analogously, the function get.go.genesets retrieves GO terms of a selected ontology (here: biological process, BP) as defined in the GO.db annotation package.

go.gs <- get.go.genesets(org="hsa", onto="BP", mode="GO.db")

User-defined gene sets can be parsed from GMT file format. Here, we are using this functionality for reading a list of already downloaded KEGG gene sets for Homo sapiens containing NCBI Entrez Gene IDs.

data.dir <- system.file("extdata", package="EnrichmentBrowser")
gmt.file <- file.path(data.dir, "hsa_kegg_gs.gmt")
hsa.gs <- parse.genesets.from.GMT(gmt.file)
length(hsa.gs)

## [1] 39

hsa.gs[1:2]

## $hsa05416_Viral_myocarditis
##  [1] "100509457" "101060835" "1525"      "1604"      "1605"      "1756"      "1981"      "1982"     
##  [9] "25"        "2534"      "27"        "3105"      "3106"      "3107"      "3108"      "3109"     
## [17] "3111"      "3112"      "3113"      "3115"      "3117"      "3118"      "3119"      "3122"     
## [25] "3123"      "3125"      "3126"      "3127"      "3133"      "3134"      "3135"      "3383"     
## [33] "3683"      "3689"      "3908"      "4624"      "4625"      "54205"     "5551"      "5879"     
## [41] "5880"      "5881"      "595"       "60"        "637"       "6442"      "6443"      "6444"     
## [49] "6445"      "71"        "836"       "841"       "842"       "857"       "8672"      "940"      
## [57] "941"       "942"       "958"       "959"      
## 
## $`hsa04622_RIG-I-like_receptor_signaling_pathway`
##  [1] "10010"  "1147"   "1432"   "1540"   "1654"   "23586"  "26007"  "29110"  "338376" "340061"
## [11] "3439"   "3440"   "3441"   "3442"   "3443"   "3444"   "3445"   "3446"   "3447"   "3448"  
## [21] "3449"   "3451"   "3452"   "3456"   "3467"   "3551"   "3576"   "3592"   "3593"   "3627"  
## [31] "3661"   "3665"   "4214"   "4790"   "4792"   "4793"   "5300"   "54941"  "55593"  "5599"  
## [41] "5600"   "5601"   "5602"   "5603"   "56832"  "57506"  "5970"   "6300"   "64135"  "64343" 
## [51] "6885"   "7124"   "7186"   "7187"   "7189"   "7706"   "79132"  "79671"  "80143"  "841"   
## [61] "843"    "8517"   "8717"   "8737"   "8772"   "9140"   "9474"   "9636"   "9641"   "9755"

See also the MSigDb for additional gene set collections.

2.5 GO/KEGG overrepresentation analysis

A variety of gene set analysis methods have been proposed Khatri et al., 2012. The most basic, yet frequently used, method is the over-representation analysis (ORA) with gene sets defined according to GO or KEGG. As outlined in the first section, ORA tests the overlap between DE genes (typically DE p-value < 0.05) and genes in a gene set based on the hypergeometric distribution. Here, we choose a significance level \(\alpha = 0.2\) for demonstration.

ora.all <- sbea(method="ora", eset=all.eset, gs=hsa.gs, perm=0, alpha=0.2)
gs.ranking(ora.all)

## DataFrame with 7 rows and 4 columns
##                                                          GENE.SET  NR.GENES NR.SIG.GENES   P.VALUE
##                                                       <character> <numeric>    <numeric> <numeric>
## 1                hsa05202_Transcriptional_misregulation_in_cancer       153           17    0.0351
## 2 hsa05412_Arrhythmogenic_right_ventricular_cardiomyopathy_(ARVC)        63            8    0.0717
## 3                                                hsa05144_Malaria        45            6    0.0932
## 4                   hsa04670_Leukocyte_transendothelial_migration        94           10    0.1220
## 5                 hsa05100_Bacterial_invasion_of_epithelial_cells        64            7    0.1620
## 6                  hsa04622_RIG-I-like_receptor_signaling_pathway        54            6    0.1780
## 7                  hsa05130_Pathogenic_Escherichia_coli_infection        43            5    0.1840

Such a ranked list is the standard output of most existing enrichment tools. Using the ea.browse function creates a HTML summary from which each gene set can be inspected in more detail.

ea.browse(ora.all)

The resulting summary page includes for each significant gene set

• a gene report, which lists all genes of a set along with fold change and DE \(p\)-value (click on links in column NR.GENES),
• interactive overview plots such as heatmap and volcano plot (column SET.VIEW, supports mouse-over and click-on),
• for KEGG pathways: highlighting of differentially expressed genes on the pathway maps (column PATH.VIEW, supports mouse-over and click-on).

As ORA works on the list of DE genes and not the actual expression values, it can be straightforward applied to RNA-seq data. However, as the gene sets here contain NCBI Entrez gene IDs and the airway dataset contains ENSEMBL gene ids, we first map the airway dataset to Entrez IDs.

air.eset <- map.ids(air.eset, org="hsa", from="ENSEMBL", to="ENTREZID")

## 'select()' returned 1:many mapping between keys and columns

## Excluded 1144 genes without a corresponding to.ID

## Encountered 22 from.IDs with >1 corresponding to.ID (a single to.ID was chosen for each of them)

ora.air <- sbea(method="ora", eset=air.eset, gs=hsa.gs, perm=0)
gs.ranking(ora.air)

## DataFrame with 8 rows and 4 columns
##                                                          GENE.SET  NR.GENES NR.SIG.GENES   P.VALUE
##                                                       <character> <numeric>    <numeric> <numeric>
## 1                                    hsa05206_MicroRNAs_in_cancer       121           71  0.000423
## 2                      hsa05410_Hypertrophic_cardiomyopathy_(HCM)        53           35  0.000765
## 3 hsa05412_Arrhythmogenic_right_ventricular_cardiomyopathy_(ARVC)        48           31  0.002670
## 4                                               hsa05218_Melanoma        50           32  0.002830
## 5                                                 hsa05214_Glioma        47           30  0.004070
## 6                                 hsa04068_FoxO_signaling_pathway       106           58  0.012800
## 7                   hsa04670_Leukocyte_transendothelial_migration        62           35  0.028600
## 8                                                hsa05144_Malaria        21           14  0.029300

Note #1: Young et al., 2010, have reported biased results for ORA on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts. The goseq package and limma::goana implement possibilities to adjust ORA for gene length and abundance bias.

Note #2: Independent of the expression data type under investigation, overlap between gene sets can result in redundant findings. This is well-documented for GO (parent-child structure, Rhee et al., 2008) and KEGG (pathway overlap/crosstalk, Donato et al., 2013). The topGO package (explicitly designed for GO) and mgsa (applicable to arbitrary gene set definitions) implement modifications of ORA to account for such redundancies.

2.6 Functional class scoring & permutation testing

A major limitation of ORA is that it restricts analysis to DE genes, excluding genes not satisfying the chosen significance threshold (typically the vast majority).

This is resolved by gene set enrichment analysis (GSEA), which scores the tendency of gene set members to appear rather at the top or bottom of the ranked list of all measured genes Subramanian et al., 2005. The statistical significance of the enrichment score (ES) of a gene set is assessed via sample permutation, i.e. (1) sample labels (= group assignment) are shuffled, (2) per-gene DE statistics are recomputed, and (3) the enrichment score is recomputed. Repeating this procedure many times allows to determine the empirical distribution of the enrichment score and to compare the observed enrichment score against it. Here, we carry out GSEA with 1000 permutations.

gsea.all <- sbea(method="gsea", eset=all.eset, gs=hsa.gs, perm=1000)  

## Permutations: 1 -- 100

## Permutations: 101 -- 200

## Permutations: 201 -- 300

## Permutations: 301 -- 400

## Permutations: 401 -- 500

## Permutations: 501 -- 600

## Permutations: 601 -- 700

## Permutations: 701 -- 800

## Permutations: 801 -- 900

## Permutations: 901 -- 1000

## Processing ...

gs.ranking(gsea.all)

## DataFrame with 19 rows and 4 columns
##                                                            GENE.SET        ES       NES   P.VALUE
##                                                         <character> <numeric> <numeric> <numeric>
## 1                  hsa05202_Transcriptional_misregulation_in_cancer    -0.530     -1.72   0.00000
## 2                                        hsa04520_Adherens_junction    -0.488     -1.73   0.00000
## 3   hsa05412_Arrhythmogenic_right_ventricular_cardiomyopathy_(ARVC)    -0.511     -1.87   0.00000
## 4                                      hsa05206_MicroRNAs_in_cancer    -0.427     -1.56   0.00201
## 5                                        hsa05416_Viral_myocarditis    -0.630     -1.88   0.00201
## ...                                                             ...       ...       ...       ...
## 15                                     hsa00790_Folate_biosynthesis    -0.734     -1.66    0.0204
## 16                                           hsa05134_Legionellosis    -0.599     -1.56    0.0271
## 17                   hsa05130_Pathogenic_Escherichia_coli_infection    -0.486     -1.51    0.0295
## 18                                               hsa04210_Apoptosis    -0.424     -1.43    0.0428
## 19                                             hsa05131_Shigellosis    -0.479     -1.52    0.0477

As GSEA’s permutation procedure involves re-computation of per-gene DE statistics, adaptations are necessary for RNA-seq. The EnrichmentBrowser implements an accordingly adapted version of GSEA, which allows incorporation of limma/voom, edgeR, or DESeq2 for repeated DE re-computation within GSEA. However, this is computationally intensive (for limma/voom the least, for DESeq2 the most). Note the relatively long running times for only 100 permutations having used edgeR for DE analysis.

gsea.air <- sbea(method="gsea", eset=air.eset, gs=hsa.gs, perm=100)  

While it might be in some cases necessary to apply permutation-based GSEA for RNA-seq data, there are also alternatives avoiding permutation. Among them is ROtAtion gene Set Testing (ROAST), which uses rotation instead of permutation Wu et al., 2010.

roast.air <- sbea(method="roast", eset=air.eset, gs=hsa.gs)
gs.ranking(roast.air)  

## DataFrame with 28 rows and 4 columns
##                                            GENE.SET  NR.GENES       DIR   P.VALUE
##                                         <character> <numeric> <numeric> <numeric>
## 1   hsa05100_Bacterial_invasion_of_epithelial_cells        60         1     0.001
## 2                        hsa05416_Viral_myocarditis        33         1     0.001
## 3               hsa03420_Nucleotide_excision_repair        42        -1     0.001
## 4                          hsa03030_DNA_replication        33        -1     0.001
## 5                        hsa04520_Adherens_junction        63         1     0.002
## ...                                             ...       ...       ...       ...
## 24     hsa04621_NOD-like_receptor_signaling_pathway        40        -1     0.032
## 25          hsa04514_Cell_adhesion_molecules_(CAMs)        60        -1     0.036
## 26              hsa04350_TGF-beta_signaling_pathway        63         1     0.037
## 27                 hsa00561_Glycerolipid_metabolism        39         1     0.040
## 28                  hsa04068_FoxO_signaling_pathway       106         1     0.049

A selection of additional methods is also available:

sbea.methods()

##  [1] "ora"        "safe"       "gsea"       "samgs"      "ebm"        "mgsa"       "gsa"       
##  [8] "padog"      "globaltest" "roast"      "camera"     "gsva"

2.7 Network-based enrichment analysis

Having found gene sets that show enrichment for differential expression, we are now interested whether these findings can be supported by known regulatory interactions.

For example, we want to know whether transcription factors and their target genes are expressed in accordance to the connecting regulations (activation/inhibition). Such information is usually given in a gene regulatory network derived from specific experiments or compiled from the literature (Geistlinger et al., 2013 for an example).

There are well-studied processes and organisms for which comprehensive and well-annotated regulatory networks are available, e.g. the RegulonDB for E. coli and Yeastract for S. cerevisiae.

However, there are also cases where such a network is missing or at least incomplete. A basic workaround is to compile a network from regulations in the KEGG database.

We can download all KEGG pathways of a specified organism (here: Homo sapiens) via

pwys <- download.kegg.pathways("hsa")

For demonstration purposes, we use a selection of already downloaded human KEGG pathways.

pwys <- file.path(data.dir, "hsa_kegg_pwys.zip")
hsa.grn <- compile.grn.from.kegg(pwys)
head(hsa.grn)

##      FROM   TO     TYPE
## [1,] "3569" "3570" "+" 
## [2,] "3458" "3459" "+" 
## [3,] "3458" "3460" "+" 
## [4,] "1950" "1956" "+" 
## [5,] "1950" "2064" "+" 
## [6,] "1950" "3480" "+"

Signaling pathway impact analysis (SPIA) is a network-based enrichment analysis method, which is explicitly designed for KEGG signaling pathways Tarca et al., 2009. The method evaluates whether expression changes are propagated across the pathway topology in combination with ORA.

spia.all <- nbea(method="spia", eset=all.eset, gs=hsa.gs, grn=hsa.grn, alpha=0.2)
gs.ranking(spia.all)

More generally applicable is gene graph enrichment analysis (GGEA), which evaluates consistency of interactions in a given gene regulatory network with the observed expression data Geistlinger et al., 2011.

ggea.all <- nbea(method="ggea", eset=all.eset, gs=hsa.gs, grn=hsa.grn)
gs.ranking(ggea.all)

## DataFrame with 4 rows and 5 columns
##                           GENE.SET   NR.RELS RAW.SCORE NORM.SCORE   P.VALUE
##                        <character> <numeric> <numeric>  <numeric> <numeric>
## 1       hsa05416_Viral_myocarditis         7      3.30      0.471   0.00699
## 2 hsa04390_Hippo_signaling_pathway        56     20.60      0.368   0.01800
## 3       hsa04520_Adherens_junction         9      3.88      0.431   0.02300
## 4               hsa04210_Apoptosis        19      7.45      0.392   0.02600

A selection of additional network-based methods is also available:

nbea.methods()

##  [1] "ggea"        "spia"        "nea"         "pathnet"     "netgsa"      "degraph"    
##  [7] "topologygsa" "ganpa"       "cepa"        "lego"

Note #1: As network-based enrichment methods typically do not involve sample permutation but rather network permutation, thus avoiding DE re-computation, they can likewise be applied to RNA-seq data.

Note #2: Given the various enrichment methods with individual benefits and limitations, combining multiple methods can be beneficial, e.g. combined application of a set-based and a network-based method. This has been shown to filter out spurious hits of individual methods and to reduce the outcome to gene sets accumulating evidence from different methods Geistlinger et al., 2016, Alhamdoosh et al., 2017.

The function comb.ea.results implements the straightforward combination of results, thereby facilitating seamless comparison of results across methods. For demonstration, we use the ORA and GSEA results for the ALL dataset from the previous section:

res.list <- list(ora.all, gsea.all)
comb.res <- comb.ea.results(res.list)
gs.ranking(comb.res)

## DataFrame with 20 rows and 6 columns
##                                                            GENE.SET  ORA.RANK GSEA.RANK MEAN.RANK
##                                                         <character> <numeric> <numeric> <numeric>
## 1                  hsa05202_Transcriptional_misregulation_in_cancer       2.6       7.7       5.1
## 2   hsa05412_Arrhythmogenic_right_ventricular_cardiomyopathy_(ARVC)       5.1       7.7       6.4
## 3                     hsa04670_Leukocyte_transendothelial_migration      10.3      17.9      14.1
## 4                                        hsa04520_Adherens_junction      20.5       7.7      14.1
## 5                                      hsa05206_MicroRNAs_in_cancer      23.1      12.8      17.9
## ...                                                             ...       ...       ...       ...
## 16                       hsa05410_Hypertrophic_cardiomyopathy_(HCM)      41.0      51.3      46.2
## 17                                     hsa00790_Folate_biosynthesis      53.8      38.5      46.2
## 18                     hsa04621_NOD-like_receptor_signaling_pathway      61.5      33.3      47.4
## 19                              hsa04350_TGF-beta_signaling_pathway      38.5      59.0      48.7
## 20                         hsa05150_Staphylococcus_aureus_infection      79.5      25.6      52.6
##      ORA.PVAL GSEA.PVAL
##     <numeric> <numeric>
## 1      0.0351   0.00000
## 2      0.0717   0.00000
## 3      0.1220   0.00765
## 4      0.2010   0.00000
## 5      0.2250   0.00201
## ...       ...       ...
## 16      0.403    0.0565
## 17      0.556    0.0204
## 18      0.631    0.0198
## 19      0.390    0.0982
## 20      0.851    0.0120