Assessment of Gene Set Enrichment Analysis using curated RNA-seq-based benchmarks

Abstract

Pathway enrichment analysis is a ubiquitous computational biology method to interpret a list of genes (typically derived from the association of large-scale omics data with phenotypes of interest) in terms of higher-level, predefined gene sets that share biological function, chromosomal location, or other common features. Among many tools developed so far, Gene Set Enrichment Analysis (GSEA) stands out as one of the pioneering and most widely used methods. Although originally developed for microarray data, GSEA is nowadays extensively utilized for RNA-seq data analysis. Here, we quantitatively assessed the performance of a variety of GSEA modalities and provide guidance in the practical use of GSEA in RNA-seq experiments. We leveraged harmonized RNA-seq datasets available from The Cancer Genome Atlas (TCGA) in combination with large, curated pathway collections from the Molecular Signatures Database to obtain cancer-type-specific target pathway lists across multiple cancer types. We carried out a detailed analysis of GSEA performance using both gene-set and phenotype permutations combined with four different choices for the Kolmogorov-Smirnov enrichment statistic. Based on our benchmarks, we conclude that the classic/unweighted gene-set permutation approach offered comparable or better sensitivity-vs-specificity tradeoffs across cancer types compared with other, more complex and computationally intensive permutation methods. Finally, we analyzed other large cohorts for thyroid cancer and hepatocellular carcinoma. We utilized a new consensus metric, the Enrichment Evidence Score (EES), which showed a remarkable agreement between pathways identified in TCGA and those from other sources, despite differences in cancer etiology. This finding suggests an EES-based strategy to identify a core set of pathways that may be complemented by an expanded set of pathways for downstream exploratory analysis. This work fills the existing gap in current guidelines and benchmarks for the use of GSEA with RNA-seq data and provides a framework to enable detailed benchmarking of other RNA-seq-based pathway analysis tools.

Fig 1. Cancer-type-specific pathways.

Procedure to identify lists of preselected, target, and positive-control pathways using TCGA-BRCA (breast cancer) as example. See S2 Table for full details.

Fig 2. Significant TCGA-BRCA positive control pathways across different weight parameter choices.

(a) Gene-set permutation with p—value < 0.05. (b) Gene-set permutation with p—value < 0.01. (c) Phenotype permutation with p—value < 0.05. (d) Phenotype permutation with p—value < 0.01. GSEA enrichment statistics: classic (“cl”), weight parameter p = 1 (“p1”), weight parameter p = 1.5 (“p1.5”), and weight parameter p = 2 (“p2”).

Fig 3. Significant TCGA-BRCA positive control pathways using gene-set (“gs”) and phenotype (“ph”) permutation approaches for different enrichment statistics.

The significance criterion was p—value < 0.05. (a) Classic (unweighted). (b) Weight parameter p = 1. (c) Weight parameter p = 1.5. (d) Weight parameter p = 2.

Fig 4. ROC curves for different GSEA and ORA approaches using 72 TCGA-BRCA positive control pathways and 1000 randomized negative controls.

(a) Gene-set permutation GSEA. (b) Phenotype permutation GSEA. (c) Signed ORA. (d) Unsigned ORA. GSEA approaches used different enrichment statistics, as indicated. ORA approaches used Bonferroni and Benjamini-Hochberg (B-H) adjusted q-values as different inclusion criteria to select differentially expressed genes, as indicated.

Fig 5. AUC across TCGA projects for different GSEA and ORA approaches using cancer-type-specific positive control pathways and 1000 randomized negative controls.

(a) Gene-set permutation GSEA. (b) Phenotype permutation GSEA. (c) Signed ORA. (d) Unsigned ORA. GSEA approaches used different enrichment statistics, as indicated. ORA approaches used Bonferroni and Benjamini-Hochberg (B-H) adjusted q-values as different inclusion criteria to select differentially expressed genes, as indicated.

Fig 6. Pathway-level enrichment evidence scores for thyroid cancer and hepatocellular carcinoma cohorts.

(a) Comparison between significant target pathways in REBC-THYR vs TCGA-THCA. (b) Contingency table of grouped EES intervals in REBC-THYR vs TCGA-THCA (Fisher’s exact test p—value = 2.5 × 10⁻⁶). (c) Comparison between significant target pathways in MO-HCC vs TCGA-LIHC. (d) Contingency table of grouped EES intervals in MO-HCC vs TCGA-LIHC (Fisher’s exact test p—value = 4.4 × 10⁻¹⁵). Distance to the diagonal is represented with increasingly darker shades of blue.

Fig 7. Gene-level enrichment evidence scores for thyroid cancer cohorts.

(a) Comparison between EES leading-edge genes in REBC-THYR vs TCGA-THCA for genes in LUI_THYROID_CANCER_CLUSTER_1, previously identified as a high-consensus tumor-enriched pathway in thyroid cancer. (b) Contingency table of grouped EES intervals for the same case as in panel (a) (Fisher’s exact test p—value = 4.6 × 10⁻⁸). (c) Network representation showing seven core thyroid cancer pathways and high-consensus leading-edge genes (|EES| ≥ 3 in at least one of the cohorts). Only genes connected to two or more pathways are shown.