Two subtle problems with over-representation analysis
Mark Ziemann1,2*, Anusuiya Bora1,2
Affiliations
Burnet Institute, Melbourne, Australia.
Deakin University, School of Life and Environmental Sciences, Geelong, Australia.
(*) Corresponding author: mark.ziemann@burnet.edu.au Burnet Institute, 85 Commercial Rd Melbourne VIC 3004, Australia.
Abstract
Background: Over-representation analysis (ORA) is used widely to assess the enrichment of functional categories in a gene list compared to a background list. ORA is therefore a critical method in the interpretation of ’omics data, relating gene lists to biological functions and themes. Although ORA is hugely popular, we and others have noticed two potentially undesired behaviours of some ORA tools. The first one we call the “background problem,” because it involves the software eliminating large numbers of genes from the background list if they are not annotated as belonging to any category. The second one we call the “false discovery rate problem,” because some tools underestimate the true number of parallel tests conducted.
Results: Here we demonstrate the impact of these issues on several real RNA-seq datasets and use simulated RNA-seq data to quantify the impact of these problems. We show that the severity of these problems depends on the gene set library, the number of genes in the list, and the degree of noise in the dataset.
Conclusion: These problems can be mitigated by changing packages/websites for ORA or by changing to another approach such as functional class scoring.
Keywords
Bioinformatics
Pathway analysis
Gene set enrichment analysis
Research rigour
Background
Over-representation analysis (ORA) is a type of functional enrichment analysis (FEA) that involves the summarisation of omics data to reflect biological differences. FEA has become one of the most popular methods in bioinformatics, collectively accumulating 131,332 citations as of late 2019 [1].
ORA involves the selection of genes of interest, followed by a test to ascertain whether certain functional categories are over-represented in the selected genes. Since its initial development in 1999 [2], ORA has proliferated widely, becoming part of many bioinformatics websites and software packages. The most popular ORA website is DAVID [3–5], followed by Ingenuity Pathway Analysis, a commercial software package provided by QIAGEN Inc. ORA has been implemented into R/Bioconductor packages including clusterProfiler
[6, 7], limma
[8] and GOseq
[9].
Another popular approach to FEA is functional class scoring (FCS). In FCS, all detected genes are ranked by their degree of differential expression followed by a statistical test for enrichment of gene categories at either extreme of the ranked list [10].
According to a study that used simulated differential expression data, FCS accuracy was superior to ORA over a range of experimental designs [11]. Despite this performance difference, a survey of articles published in 2019 showed that seven of the top eight most popular tools were based on ORA [12]. The discrepancy between performance and popularity is likely due to the relative ease of conducting ORA as compared to FCS. A minimal ORA might involve pasting a list of gene symbols into a text box on a website, with results appearing nearly instantly. On the other hand, FCS typically involves installing specific software and dealing with a dataset representing every detected gene (~20,000 rows) and ensuring that the file format from upstream tools is compatible with FCS packages.
Despite the popularity of ORA, there are concerns that this technique is being misused. Due to technological and biological reasons, not all genes are detected in transcriptomic studies, meaning that some genes will more readily appear as differentially expressed. To address this bias, a background list of detected genes is required for comparison to the list of selected genes (foreground list) [13]. Failure to use a background gene list leads to dramatic changes in enrichment results, rendering them unreliable [12, 14]. Unfortunately, use of an appropriate background list is reported in only a small fraction (~4%) of peer-reviewed articles describing enrichment analysis [12]. Correction of p-values for multiple testing is also crucial in controlling the false positive rate as enrichment analysis typically involves thousands of parallel tests [13]. But appropriate p-value correction is reported in only ~54% of studies [12]. These issues have resulted in the development of best practices for end users and stronger reporting standards [15].
The focus of this work is to raise awareness to two existing problems in ORA implementations of some popular enrichment tools. The first problem is that genes belonging to the background list are removed from the analysis entirely if they are not annotated as belonging to any functional categories. This results in the background list appearing smaller than it should, leading to an underestimation of fold enrichment scores and significance values. In principle, this problem would differentially affect analyses involving different gene set libraries, with smaller libraries such as Kyoto Encyclopedia of Genes and Genomes (KEGG) [16], which describes 2,727 genes affected more severely as compared to larger libraries like Gene Ontology [17] which describes 19,428 genes (figures obtained from MSigDB June 2024 [18]).
The second problem is that adjusting p-values for multiple testing, also known as false-discovery rate correction, is sometimes implemented improperly. Determining which pathways should be subject to false discovery rate correction depends on how a detection threshold is set. Ideally, a pathway should be considered detected based on the presence of a predetermined number of member genes in the background set. However, some tools use the foreground set to define a pathway as detected. This is problematic, as the act of trying to find an intersection between a foreground gene list and a pathway is a test, even if there are no common genes. This results in pathways being discarded from the analysis if no common genes are found in the foreground, despite many genes being present in the background list. This artificially reduces the number of tests conducted and makes FDR values appear lower than they should, potentially raising the rate of false positives. We expect this problem to have a more severe effect when the foreground list is small.
Here, we define the effect of these two problems on real RNA-seq based enrichment analysis results and demonstrate how these two issues impact performance using simulated expression data. Finally, we classify popular tools regarding these two problems and provide general recommendations for end users.
Methods
Expression data preparation
To quantify the effect of these two problems on RNA-seq based enrichment analysis, data from seven transcriptomic studies were downloaded from DEE2 using the getDEE2 R/Bioconductor package [19]. Raw data for these seven studies are available from NCBI Sequence Read Archive under the accession numbers SRP128998, SRP038101, SRP037718, SRP096177, SRP097759, SRP253951 and SRP068733 [20–25]. For each dataset, kallisto (v0.43.1) [26] transcript counts were aggregated to the gene level. Genes with fewer than 10 reads per sample on average were removed from downstream analysis. Remaining genes passing this selection were included in the background gene set. Differential expression analysis was conducted with DESeq2 v1.44.0 [27], with genes identified by their Ensembl identifiers. Gene symbols were then fetched using biomaRt v2.60.0, based on Ensembl version 112 [28].
Examining the background problem using real transcriptome data
Human gene sets were obtained from MSigDB v2023.2 [18], and included KEGG Medicus [16], Reactome [29], Wikipathways [30], MicroRNA targets from miRdb [31], transcription factor targets from GTRD [32], Gene Ontology (GO) [17], Human Phenotype Ontology (HPO) [33], Cell markers and Hallmark pathways [18]. To demonstrate the background problem, differentially expressed genes with FDR<0.05 were subset, with separate lists for up and downregulated genes. If fewer than 200 genes met this criterion, then the 200 genes with the smallest p-values were selected, as suggested by a previous report [34]. These gene lists were subjected to ORA using clusterProfiler’s enricher
function v4.12.0, using a minimum set size of 5 and a no maximum set size.
To systematically assess the impact of the background problem, we conducted ORA with clusterProfiler using the different gene set libraries described above for the seven RNA-seq datasets. ClusterProfiler was used with the same parameters as above, and we devised a simple workaround to the background problem, which is to append the entire background list as a gene set to the library. This forces clusterProfiler to retain all background genes. The ORA analyses were then filtered for significant sets using an FDR threshold of 0.05. Jaccard index was used to compare original and corrected analyses.
Examining the FDR problem
For this test, 2000 of the top up and down regulated genes were selected from each differential expression dataset based on p-value. A workaround for the FDR problem was devised by filtering the gene sets by the presence of at least two genes in the background list prior to running clusterProfiler. After running clusterProfiler using the parameters above, the number of gene sets in the results was quantified, and compared to the number that should have been reported using the two gene limit. The difference between these numbers (n) represents the number of missing gene sets. To account for these tests, n values of 1 were appended to the p-values obtained by clusterProfiler, followed by FDR correction in R, to obtain the properly corrected FDR values. As above, analysis of seven datasets was done with nine gene set libraries, the significant sets were selected and the Jaccard statistics were collected for each run to compare original and corrected analysis.
As we predicted that this effect could be more severe for smaller foreground lists, we performed parallel analysis with foreground gene lists with sizes varying between 125 and 2,000 genes, as ranked by FDR value.
Quantifying the effect of these problems on precision and recall using simulated expression data
RNA-seq counts for SRA run accession ERR2539161 were obtained from DEE2 [35]. This dataset was selected due to its high number of assigned reads (367M). A library of 200 random gene sets were generated, each containing 30 genes. To generate differential expression profiles, 5% of gene sets were selected to be differentially expressed, with equal numbers of up- and down-regulated pathways. Pseudosamples were generated by first downsampling the counts using the thincounts
function of the edgeR
package v4.2.0 [36] to 20M reads, followed by the addition of some extra variation drawn from a normal distribution with mean of 1 and standard deviation varying between 0.1 and 0.6, this value selected randomly for each gene. Three “control” and five “case” pseudosamples were generated. For the genes selected to be differentially expressed, expression values were multiplied by a log fold change in the case samples of 0.5 for upregulated genes and -0.5 for downregulated genes. Genes that were selected to be both up and downregulated were left unchanged. With these expression profiles, DESeq2 was used for differential analysis, and the results were subjected to clusterProfiler which suffers from both problems, clusterProfiler with the background problem fix, clusterProfiler with the FDR problem fix, clusterProfiler with the fixes to background and FDR problems, an ORA function called fora
from fgsea
bioconductor package v1.30.0 that does not suffer these problems, and fgsea
a fast implementation of the preranked gene set enrichment analysis (GSEA) algorithm [37]. For the ORA functions above, significantly differentially expressed genes with FDR 0.05 were selected for the foreground gene set. If fewer than 200 genes were detected as significant, then 200 genes with the smallest p-values were used for ORA, like a previous report [34]. For fgsea
, the DESeq2 test statistic was used for scoring differential expression. After enrichment analysis with these approaches, the significant up and down regulated sets were selected and compared to the ground truth to calculate precision, recall and F1 score. At each value of added variance, 1000 replications were conducted. No set seed was used, but results were similar when run on different computer systems. All analyses were undertaken with R v4.4.0 in a Docker container.
Comparing ORA and FCS methods by downsampling real RNA-seq data
To understand the sensitivity and false positive rate of ORA and FCS methods, we employed a downsampling approach to a previously published dataset derived from 37 lung cancer patients, having tumor and normal tissues for each patient [38]. These RNA-seq based gene expression counts were downloaded from NCBI GEO under accession number GSE158420. After reading them into R, we noticed some gene names converted to dates [39], so we used the HGNChelper package to fix them [40]. DESeq2 was then used for differential expression analysis correcting for patient-of-origin to identify genes differentially expressed in tumor compared to normal tissue. FCS was conducted using fgsea using the DESeq2 test statistic to score each gene. ORA was conducted using fora, using genes with differential expression FDR<0.05. Separate tests were conducted for up- and down-regulated genes. For the ORA background, genes with a mean of 10 or more reads per sample on average across the comparison were included. Reactome pathways were used. For downsampling, sample size was varied between 2 and 30 patients, selected pseudorandomly with a set seed, followed by parallel analysis with fgsea and fora pipelines, using a significance filter of FDR<0.05. We then calculated the number of pathways identified in downsampled data that were consistent with the full dataset. The number of consistent pathways gives an indication of sensitivity of each method. Then, we calculated the number of pathways identified in the downsampled data that were inconsistent with the full dataset. The proportion of significant pathways classified as inconsistent provides an estimate of the real false discovery rate. This was repeated 50 times at each sample size, and the median values were reported.
Results
Understanding the background problem
A workaround to eliminate the background problem was developed and we used it to compare the results of ORA of seven datasets with the original and corrected methods for nine different gene set libraries. The mean number of statistically significant results was calculated for each of the gene set libraries before and after correcting the problem (Figure 1A). Each of the gene set libraries showed an uplift when correcting the background problem. The smallest increase was observed in Cellmarkers gene library (11.8%), while the largest increase was observed for HPO (224%). The similarity between results, as quantified with the Jaccard similarity, was highest for Cellmarkers and GTRD at 0.76 and 0.77 respectively, and the lowest similarity was recorded for HPO and KEGG with 0.20 and 0.48 respectively. Observed Jaccard scores correlated with the number of genes annotated to one or more functional categories in the gene set library (Figure 1B). This result suggests that ORA with smaller gene set libraries like HPO and KEGG are impacted more severely.
Figure 1A. Effect of correcting the background problem. Figures represent the mean number of statistically significant sets (FDR<0.05) over seven independent datasets detected with clusterProfiler with and without correcting the background problem. Separate tests for up- and down-regulated gene lists were conducted. The percent increase in significant sets are shown.
Figure 1B. Impact of the background problem is worse for smaller gene set libraries. Mean Jaccard similarity index was calculated for original and corrected ORA for seven independent datasets. The number of genes represented in one or more functional sets is shown on the x-axis, and the mean Jaccard index is shown on the y-axis.
Understanding the FDR problem
A workaround for the FDR problem was developed, then original and corrected analyses were conducted for the seven datasets and nine gene set libraries. Results showed no major differences in the number of significant sets between original and corrected analyses (Figure 2A).
We posited that the foreground gene list length might exacerbate the problem, so we performed parallel analysis with foreground lists of length between 125 and 2,000 genes (Figure 2B). A clear trend was observed for GO, Reactome, Cellmarkers and TFT GTRD, with Jaccard similarity being significantly reduced when foreground gene lists were shorter than 500 genes. Interestingly, the drop in Jaccard similarity was more severe for GO and Reactome sets as compared to Cellmarkers and TFT GTRD sets and may be associated with the size of the gene sets in the library. Median gene set size for these are GO:18, Reactome:11, TFT GTRD:287 and Cellmarkers:115. These results demonstrate that the FDR problem is more severe when dealing with shorter foreground lists, and this is exacerbated when functional categories in the gene set library are also small.
Figure 2A. Effect of correcting the FDR problem. Figures represent the mean number of significant sets with and without correction of the FDR problem over seven independent RNA-seq datasets. Foreground list consisted of 2000 genes with the smallest p-values
Figure 2B. Impact of the FDR problem is worse for shorter foreground gene lists. Mean Jaccard index between original and corrected analysis were calculated for seven RNA-seq datasets and nine gene set libraries. Only GO, Reactome, Cellmarkers and TFT GTRD are shown.
Impact of these two problems on accuracy determined with simulated data
Simulated differential expression profiles with a priori changes to predetermined pathway genes were analysed in parallel by clusterProfiler default and with mitigations for problems described above. In addition, these data underwent ORA with fora which does not suffer from these two problems. For comparison, fgsea was used to compare the performance of these two ORA methods against a prototypical FCS method. These methods were tested with different levels of introduced variation and the mean precision, recall and F1 scores were recorded (Figure 3A). Mitigating the background problem improved recall (by 1.8%), however it caused a similar-sized decrease in precision (1.6% lower). Addressing the FDR problem improved precision by 9.9%, but the recall dropped by 7.7%. Addressing both problems improved precision by 8.2% but decreased recall by 5.4%. Under these conditions, the FDR problem affected results more severely than the background problem. Results obtained for clusterProfiler after implementing the two mitigations were identical to fora to three significant figures. fgsea recorded precision of 0.940, similar to the ORA methods, but recall was far superior, with fgsea scoring 22.9% higher than default clusterProfiler and 28.3% higher than fora. As a result, fgsea’s overall accuracy as defined by F1 index was 21.0% higher than fora and 18.6% higher than default clusterProfiler.
When looking at the accuracy profile at different levels of added noise (Figure 3B), ORA techniques show high precision when the added noise is low (<0.2), but precision drops at higher levels of added noise. The rate of drop in precision was more severe for approaches that suffer from the FDR problem. The precision of fgsea was relatively stable despite high levels of added noise, ranging from 0.95 to 0.92. The recall of fgsea was superior to all ORA approaches across all levels of added noise.
Figure 3A. Impact of correcting these two problems on accuracy of simulated expression data. Precision, recall and F1 scores are the mean of trials at a range of added noise. CP; clusterProfiler, CP BG fix; clusterProfiler with a workaround to the background problem, CP FDR fix; clusterProfiler with a workaround to the FDR problem, CD BG and FDR fix; clusterProfiler with mitigations for both background and FDR problems, fora; ORA function of the fgsea package, fgsea; an FCS method similar to GSEA.
Figure 3B. Accuracy of clusterProfiler, fora and fgsea with simulated data with various levels of added noise. Precision, recall and F1 score are shown with the addition of different amounts of random noise to expression levels. The purple line for clusterProfiler with mitigations for both problems is not visible as the performance was identical to fora.
Exploring sensitivity in real RNA-seq data
To understand whether FCS has better sensitivity than ORA in a real setting, we conducted a systematic downsampling of a gene expression dataset cosisting of normal and tumor samples of 37 lung cancer patients. In the full set of patient samples, fora (ORA) yielded 375 significantly differentially regulated Reactomes, while fgsea (FCS) gave 470 (FDR<0.05). Fgsea yielded consistently more significant pathways than fora in downsampled data (Figure 4A). Fgsea also identified a larger proportion of pathways in downsampled data compared to fora. For example at a sample size of 5, fgsea identified 74.5% of pathways from the full set as compared to 41.5% for fora, confirming better recall for FCS over ORA. It was noted that fgsea identified more pathways in downsampled data that were inconsistent with the full dataset as compared to fora (Figure 4B). However, inconsistent pathways were a smaller proportion of all findings for fgsea as compared to fora (Figure 4C), indicating an overall lower false positive rate for FCS compared to ORA methods.
Figure 4. Downsampling confirms higher sensitivity for FCS over ORA. (A) The number of significant pathways identified in downsampled data that are consistent with findings in the full dataset, using FCS (fgsea) and ORA (fora) methods. The blue and red horizontal lines indicate the number of significant pathways identified in the full set of 27 patients using fgsea and fora respectively. (B) The number of significant pathways identified in downsampled data that are not consistent the full dataset. (C) The proportion of pathways in downsampled data that are inconsistent with the full dataset. All values are the median of 50 simulations.
Discussion
The background problem came to our attention when a doctoral student in our department was puzzled by results obtained from clusterProfiler when using a custom set of pathways. This exact issue was the subject of several posts to online bioinformatics forums dating to ~2020 where confused users were wondering why the number of background genes changed so much. The results shown here indicate the background problem affects smaller gene set libraries like KEGG and HPO more severely as compared to Cellmarkers and GO libraries which describe a larger proportion of all genes. ClusterProfiler maintainers recently provided an option to prevent the loss of unannotated genes from the background options(enrichment_force_universe = TRUE)
, however this feature isn’t yet described in the official documentation.
The FDR problem is potentially more serious in that it could lead to increasing false positives under specific conditions; when the foreground list is short, and the median pathway size is relatively small. These false positives arise due to an underestimate of the actual number of tests conducted. In extreme cases, correction of this problem can lead to around half of the statistically significant results changing (see Figure 2B).
So end-users can explore the effect of these two problems on their work, we have developed an R/Shiny web application (see “Availability of data and materials” below). Users can upload their foreground and background lists, select their preferred gene set library and compare results from default and patched analyses.
ClusterProfiler is not the only tool to suffer from these problems, they are widespread (Table 1). Still, there are some tools free of such problems. For example, fora is a good option for R-based workflows, however the function does not calculate the fold enrichment score, so users will need to do it themselves [37]. ShinyGO is an excellent alternative for web based analysis [41]. It doesn’t suffer from the two issues described here, it is easy to use, and has a high degree of reproducibility thanks to its stand-alone versioned docker container option for local execution. One problem that all of these tools listed in Table 1 have is that a background list is optional. Sometimes this feature is hidden under “advanced options.” Given the strong biological and technical biases that are present in contemporary omics such as including single cell and spatial transcriptomics, it should be a mandatory step for users to provide a custom background list.
Tool | Version | Type | Provides FDR values | Provides enrichment score | Proper background handling | Proper FDR | Reference |
---|---|---|---|---|---|---|---|
DAVID | v2023q4 | Web | Yes | Yes | No | Yes | [5] |
Panther | v19 | Web | Yes | Yes | Yes | Yes | [42] |
Enrichr | June 2023 | Web | Yes | Yes | Yes | No | [43] |
KOBAS-i | KOBAS 3.0 | Web | Yes | No | Yes | No | [44] |
WebGestalt | 2024 | Web | Yes | Yes | No | Yes | [45] |
g:Profiler | Feb 2024 | Web | Yes | No | No | Yes | [46] |
STRING-DB | v12.0 | Web | Yes | Yes | Yes | Yes | [47] |
ShinyGO | v0.80 | Web | Yes | Yes | Yes | Yes | [41] |
BinGO | v3.0.5 | Cytoscape | Yes | No | No | No | [48] |
ClueGO | v2.5.10 | Cytoscape | Yes | No | No | No | [49] |
clusterProfiler | v4.12.0 | R package | Yes | No | No | No | [7] |
topGO | v2.56.0 | R package | No | No | No | Yes+ | [50] |
GOseq | v1.56.0 | R package | No | No | No* | Yes+ | [9] |
goana/kegga | limma v3.60.3 | R package | No | No | No** | Yes+ | [8] |
fora | fgsea v1.30.0 | R package | Yes | No^ | Yes | Yes | [37] |
The purpose of this work was to thoroughly characterise two subtle problems with ORA, but there are some bigger problems noted. Firstly, all variations of ORA were less accurate as compared to FCS in this simulation work, confirming a previous report [11]. Specifically, precision of ORA deteriorated with greater levels of added noise, while FCS was relatively stable. This means researchers can be more confident with FCS results, even if the data are noisy. The higher recall of FCS demonstrated in simulated differential expression profiles and confirmed with downsampled real cancer RNA-seq data means that researchers will have more power to identify pathways relevant to the biological processes being studied. Secondly, many popular ORA tools listed in Table 1 do not provide enrichment scores in their default output tables. Enrichment scores are a surrogate for effect size, and without them, users may over-rely on significance scores when interpreting enrichment results. Lastly, ORA results can vary dramatically depending on the number of genes selected in the foreground, a mostly arbitrary choice.
These issues together with previously mentioned pitfalls such as lack of background correction and methodological misreporting [12, 14], suggest that if the data can be scored/ranked, methods like FCS are recommended over ORA. ORA is better suited to other cases where genes are classified into binary groups. For example when investigating the over-represented categories in gene expression modules [51], or with variant-harbouring genes associated with a disease [52].
Declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Availability of data and materials
Publicly available data were obtained from Digital Expression Explorer 2 (dee2.io) and NCBI GEO using the accession numbers mentioned in the Methods section.
Reproduction of the results here can be achieved using the GitHub code repository GitHub (https://github.com/markziemann/background) and Docker image (https://hub.docker.com/r/mziemann/background), which have both been deposited to Zenodo for long-term preservation (https://zenodo.org/record/TODOXXX).
The R/Shiny tool for comparing ORA with and without the two problems is currently available at https://oratool.ziemann-lab.net/ and available as a downloadable docker container for local execution (https://hub.docker.com/repository/docker/mziemann/background_app/general).
Competing Interests
The authors declare that they have no competing interests.
Funding
This work was supported by Burnet Institute and Deakin University. Authors declare that no grants were involved in supporting this work.
Acknowledgements
We thank Ms Megan Soria and Ms Kaumadi Wijesooriya (Deakin University) for discussions on the concept and manuscript. This research was supported by use of the Nectar Research Cloud, a collaborative Australian research platform supported by the NCRIS-funded Australian Research Data Commons (ARDC). The authors gratefully acknowledge the contribution to this work of the Victorian Operational Infrastructure Support Program received by the Burnet Institute.