Set-based enrichment methods are an integral part of the analysis of high-throughput expression data, such as those originating from transcriptomics or proteomics experiments. Enrichment methods allow the identification of molecular pathways, Gene Ontology (GO) terms, and other gene sets that might play a role in the disease of interest. Most enrichment methods are subsequently conducted for differential expression analysis; that is, they rely on the ranking of genes after comparing two groups of samples. Statistical tests are used to determine whether the genes of a particular set are disproportionately highly enriched among the differentially expressed genes (DEGs) (Beissbarth and Speed, 2004; Subramanian et al., 2005; Ackermann and Strimmer, 2009). This contrasts with self-contained gene set tests, which are based on subsets of expression data related to a particular gene set (Goeman et al., 2004; Hummel et al., 2008; Jung et al., 2011; Bayerlová et al., 2015).

Set information for enrichment analyses is usually obtained from public databases, for example, on molecular pathways or GO terms. The most commonly used databases are the “Reactome pathway knowledgebase” (Fabregat et al., 2016) (https://reactome.org/), “Kyoto Encyclopedia of Genes and Genomes” (Kanehisa and Goto, 2000) (KEGG; https://www.genome.jp/kegg/), “WikiPathways” database (Kelder et al., 2012) (https://www.wikipathways.org), “GO” knowledgebase (Consortium, 2004) (http://geneontology.org/), and the “Molecular Signatures Database” (Liberzon et al., 2011) (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb). In the GO database, a particular GO term comprises a set of genes that can be assigned to a biological process (BP), molecular function (MF), or cellular component (CC).

The contents of the databases are curated either automatically by computer algorithms or manually by experts (Consortium, 2004). Specifically, WikiPathways provides community-based curation by registered contributors (Kelder et al., 2012). An example of a database where curation is done both ways, manually and computationally, is the MSigDB. Furthermore, pathway membership can be experimentally validated or predicted computationally. However, none of the modes of curation can prevent the uncertainties remaining regarding the membership of individual genes to a particular pathway (Gillis and Pavlidis, 2013). This is important because all enrichment analyses rely on the correctness of the database information, and the results of enrichment analyses would change if features of a set are removed or added. This can especially happen when the database information is retrieved at different times. For example, the GO database contained 42,442 terms classified as valid and 5,287 classified as obsolete in January 2024. Two months before, only 4,889 terms were classified as obsolete, meaning that nearly 400 terms would have to be reconsidered when a Gene-set enrichment analysis (GSEA) is performed after January 2024. In addition, the WikiPathways database reports roughly between 100 and 700 edits per month. Furthermore, in the KEGG database, complete pathways can be merged, leading to a large number of changes. For example, the KEGG pathway map00471 has been deleted and then added to the KEGG pathway map00470 (“D-amino acid metabolism”).

In this work, we present bootGSEA as a novel bootstrap approach to repeatedly sample subsets of pathways or other gene sets to study whether a result remains significant when the set composition is changed. The ranking lists of the gene sets obtained from each bootstrap replicate were aggregated using a score that can be used for a new ranking list. The analyst can then compare the original ranking with the bootstrap-based ranking list to study whether the association of a pathway or GO term with the disease gains or loses evidence. A similar approach was proposed by Schmid et al. (2016), who generated a robustness score for each gene set using random subsets of gene sets. In contrast to their approach, our method results in a new ranking of gene sets that can be helpful in aggregating findings from different independent studies or different omics levels. Thus, our approach for multi-omics follows the idea of aggregating the different omics levels after performing primary analysis on the individual levels first. This way of multi-omics analysis has also been implemented in other studies. For example, Wang et al. (2014) fused networks that were first derived on individual omics levels, Xiong et al. (2012) integrated genetic and transcriptomics results in a joint score, and Kang et al. (2022) fitted neural networks from individual omics levels and merged them into a joint model.

We also demonstrated that this method is useful for aggregating the results of enrichment analysis from different omics domains in the same experiment. We applied our new bootstrap pipeline to a single-omics scenario (transcriptomics only) comprising six independent renal cancer datasets. This example was used to show how our bootstrap pipeline can help study the robustness of GSEA when comparing results from multiple independent datasets from studies on a similar research question. In addition, we analyzed our multi-omics kidney data (transcriptomics and proteomics) from our research consortium on spinal muscular atrophy (SMA). The data were obtained from a SMA mouse model to demonstrate the usefulness of our approach when comparing GSEA results between different omics levels. SMA is a monogenic disease caused by the mutation or deletion of the survival motor neuron 1 (SMN1) gene (Lefebvre et al., 1995). The disease is characterized by the degeneration of motoneurons, with the subsequent atrophy of skeletal muscles to muscular atrophy since the SMN affects all tissues, which also include non-skeletal muscles. Moreover, SMA is a multi-system disorder that also affects peripheral organs, such as the kidney (Allardyce et al., 2020). Three treatment methods are available, all increasing SMN expression. The SMN is expressed ubiquitously and has several important cellular functions, including snRNP assembly, R-loop resolution, and regulation of the actin cytoskeleton and translation (Hensel et al., 2020). Therefore, SMA is a highly complex disease with expected dysregulations in pathways in several cell types and on several molecular levels.

In this section, we describe the analysis pipeline, including the approach for the bootstrap step used to repeatedly analyze different random subsets of the data. Furthermore, the rank aggregation step and examples of transcriptomics and proteomics data are presented. The complete workflow is shown in Figure 1. All the analyses were implemented in the R programming environment [www.r-project.org, version 4.2 (R Core Team, 2022)].

The bootGSEA pipeline used the results from the differential expression analysis as input, and GSEA was performed separately for the following types of gene sets: GO terms (BP, MF, and CC), KEGG, Reactome, and WikiPathways databases. The enrichment analysis performed in this pipeline provided two ranking lists of gene sets as outputs: one list from the original analysis of the complete data and one aggregated list from the analysis of the bootstrap samples. The original analysis included the entire list of genes or proteins from differential expression analysis. The bootstrap analysis involved taking the subsampled lists of genes or proteins for enrichment analysis, providing B additional ranking lists of results for enrichment analysis. To further determine the robustness of the bootstrap analysis, we aggregated the B lists to make them comparable to the ranking of the original analysis. All analyses were initially performed at a single-omics level.

For multi-omics analysis, rank scores from single-omics levels (transcriptomics and proteomics) were further integrated by rank aggregation to obtain an integrated score for the pathways or GO terms retrieved.In the following sections, we describe the results of the differential expression analysis, original findings of the GSEA, bootstrap enrichment analysis, and aggregated results from the six cancer datasets and the two omics levels from the SMA data. First, the results for the renal cancer data are shown, followed by the results for the multi-omics data from SMA mouse kidney samples.To compare the ranking obtained from the original and bootstrap GSEA, we mainly described rank gains and losses of individual GO terms or pathways. We avoided using the correlation coefficient since a correlation of, for example, 0.90 sounds high but can still include large rank differences. Only for the comparison between transcriptomics and proteomics in the SMA example did we use Kendall’s τ to describe the advantage of the bootstrap approach.

Set-based enrichment tests are a necessary part of omics data analyses to better understand the biological meaning of differentially expressed genes, proteins, or other molecules. However, because the set compositions provided in the databases are subject to uncertainties, incorrect pathways can emerge from the analysis and lead to biological interpretations in the wrong direction. We observed in the two data examples of this work that pathways or terms obtained from several independent datasets of the same omics domain or different omics domains can have only a moderate overlap. For example, in the SMA dataset, some pathways or GO terms were selected from either transcriptomics or proteomics data analysis. Here, cellular component terms ordered by their integrated score comprise of complex and subunits of complex terms that are indirectly related to neuromuscular diseases. The SMN is a part of the SMN complex and interferes functionally with several other complexes, such as the cytoplasmic ubiquitin ligase complex (GO:0000153) or ER ubiquitin ligase complex (GO:0000835) (Chaytow et al., 2018), which are highly ranked when ordered by their integrated score but were ranked much lower using the original ranking (Table 4). Furthermore, the top GO-BP term (Supplementary Table S2 (BP)) with the highest integrated rank (GO:0000245, spliceosomal complex assembly) has been associated with SMA (Price et al., 2018) but was ranked much lower with the standard GSEA of the single-omics analyses (transcriptomics rank: 1,596; proteomics rank: 81). Furthermore, vitamins play, in general, an important role in neurodegenerative disorders. Lack of water-soluble vitamins (GO:0006767, water-soluble vitamin metabolic process, integrated rank: 6) (Supplementary Table S2 (BP)) can lead to neurological diseases (Rai et al., 2021). In particular, vitamin B6 (GO:0042816, vitamin B6 metabolic process, integrated rank: 2) is necessary for the production of various neurotransmitters such as serotonin, dopamine, and γ-aminobutyric acid (GABA). Deficiencies in vitamin B6 have been linked to depression and impaired brain function, such as epilepsy (Hellmann and Mooney, 2010). This suggests that it might be beneficial to consider providing sufficient supplementation of nutrients involved in maintaining an optimal methylation state, including folic acid, vitamin B12, and vitamin B6, for individuals with SMA (Fitzgerald and McArdle, 1941; Friesen et al., 2001). Therefore, our new pipeline provides not only robust terms but also biologically relevant terms when ordered by their integrated score.

Several attempts have been made to obtain a more robust enrichment analysis, for example, by integrating information about pathway topology (Draghici et al., 2007; Glaab et al., 2010; Massa et al., 2010) or GO hierarchy (Alexa et al., 2006) into the algorithms. Other approaches use sample permutations (Subramanian et al., 2005) or comparisons with results for random gene sets (Kim and Volsky, 2005) to account for uncertainties. The GAGE method (Luo et al., 2009) improves the robustness of GSEA by treating curated gene sets as either a pathway or an experimentally derived differential expression set. A method to evaluate the contribution of individual features to the significance of a gene set was presented as well by Jung et al. (2011).

Here, we present a new approach that combines bootstrap analysis at the gene set level with rank aggregation. This approach accounts for the uncertainty in set compositions by repeatedly analyzing subsets of each gene set. The percentage of genes or proteins to be selected for bootstrap can be chosen by the user of our R package. In case the user assumes much uncertainty in the database, a lower percentage should be taken. A major advantage of this approach is that it can be easily combined with other GSEA approaches. Exemplarily, we demonstrated the combination of our approach with pathway and GO term enrichment analyses implemented using the R packages clusterProfiler and topGO. We showed that overlaps of the detected GO terms between independent datasets were larger when using the bootstrap approach instead of the ranking from the standard analysis (Figure 3; Supplementary Figure S36). Similarly, we showed a higher rank correlation for the detected GO terms from transcriptomics and proteomics in the SMA data. Although increased overlaps were not observed for all types of sets, overlaps obtained using the bootstrap analysis were in no case smaller than the overlaps from the standard analysis.

In contrast to other approaches to account for uncertainty in GSEA, the rank aggregation step of our pipeline also allows the combination of results from multiple datasets or omics domains. This could also improve the stability and reproducibility of the findings.

In contrast to the approach proposed by Schmid et al. (2016), we derived not only a measure for the robustness of the result for each set but also provided a new ranking of sets. The size of either a gain or loss can be used as a measure of robustness. While we used the 2.5% and 97.5% quantiles, the user is free to use other thresholds to identify sets with extreme gains or losses. Nevertheless, thresholds are useful to account for different numbers of sets in an analysis. When having an overall large number of sets, the values of the quantiles will be larger as well, meaning that larger gains or losses are allowed before flagging a gain or loss as extreme. A disadvantage of the current rank aggregation approach is that the new ranking is based on a score and not on a p-value. Therefore, it is a bit more difficult to specify a threshold for the selected sets.

To conclude, set-based analyses now have a long history of omics data analysis to facilitate the biological interpretation of selected features from differential expression analysis. However, the large number of different computational GSEA methods presented in the last two decades and the huge databases with pathway annotations provide an unmanageable number of possible results, and analysts may be conventional in their biological interpretations. Moreover, some entries in the databases may be less supported by experimental findings or by the literature than other entries. In this regard, our bootstrap approach can help separate less robust findings from more robust findings. The rank aggregation step can additionally help combine gene set results from multiple datasets of the same or different omics levels. In particular, a GO term or pathway is only highly ranked by the integrated score if there is evidence for the importance of a term or pathway from different omics levels. We demonstrated the use of our approach in an example with transcriptomics and proteomics data, but it could be extended by GSEA from other omics domains, such as metabolomics (Mahajan et al., 2024). The rank aggregation step also supports the idea of research synthesis, that is, integrating findings from different studies or data sources to obtain a higher level of scientific evidence. Our new pipeline bootGSEA is universal as it can be combined with the most common GSEA methods. However, when using “topGO” for GO analysis, which works in the sense of overrepresentation analysis, users must keep in mind that the results depend on thresholds for differentially expressed features.

As a future extension of our approach, we also consider to not only remove features of pathways but to also move features between pathways, which is also an action we observed in databases. This can, however, only be done using biological information about whether a pathway feature makes biological sense in a certain pathway.