For each dataset, we created lists with 500, 200, 100, and 50 genes. Furthermore, Entrez ID gene lists were created for all lists by converting their gene symbols using the biomaRt library in R (Durinck et al., 2005; Durinck et al., 2009). Both lists, with gene symbols and Entrez IDs, were used because some tools either accepted or displayed imprecise results for one of the identifiers. Our datasets can be acquired from [https://github.com/LARA-Lab-Aging].

The random lists were generated by selecting 500, 200, 100, and 50 protein-coding Homo sapiens genes at random. This process was repeated 5 times to produce 5 different lists for each size.

The Hallmark lists were built from the gene sets available in the MSigDB human hallmarks collection (Subramanian et al., 2005; Liberzon et al., 2011; 2015). To create the 500 genes Hallmark list, we combined the HALLMARK HYPOXIA, HALLMARK DNA REPAIR, and HALLMARK UV RESPONSE UP gene sets. The 200-gene list contained only genes in the hypoxia set. The 100 and 50-gene lists were built by downsampling the 200-hypoxia dataset.

The GOBP lists were generated by combining gene sets from different GO Biological Process datasets in the MSigDB Ontologies collection. To build the lists, we aimed to select ontologies with minimal overlap. The ontology gene sets that compose each gene list are gathered in Supplementary Table S1. Additionally, all lists had their duplicates removed.

The Contextual dataset was obtained from the lung cancer GSE18842 dataset available in the Gene Expression Omnibus (GEO) [https://www.ncbi.nlm.nih.gov/geo/], which comprises 46 cancer and 45 control samples (Sanchez-Palencia et al., 2011). Firstly, the dataset was imported into Gene Expression Analysis Platform (GEAP) software, which employs underlying R-based tools through a graphical user interface to perform microarray data analyses (Nunes et al., 2022). Quality analysis was conducted to filter out low-quality data within GEAP, which utilizes the arrayQualityMetrics package internally. Samples that failed the quality metrics of at least two of the three metrics analyzed by arrayQualityMetrics would be discarded before the Differential Gene Expression (DGE) analysis; however, no samples had to be discarded. Next, DGE analysis was conducted through the “Comparison Between Two Groups” tab in GEAP, which employs the limma package. The parameters used were the default eBayes method and the False Discovery Rate (FDR) correction method (Benjamini and Hochberg, 1995). Finally, the results were filtered for logFC >1 and p-value <0.05, yielding a table with 3,222 DEG (overexpressed = 1,386; underexpressed = 1836). The top 500 DEGs were selected to build the 500 contextual gene list. The same downsampling process as described before was used to create the smaller lists.

Tools were selected based on the following criteria: (i) tools widely used in the scientific community (≥500 citations on Google Scholar), (ii) tools that have been updated in the past 5 years or allow the usage of up-to-date GO annotations and ontology files, and (iii) tools that are functioning as described in their available documentation. For the FEA, we attempted to maintain the parameters as similar as possible and used the default options when no equivalents were available. All analyses were conducted using human GO annotations. The whole human genome was used as background for the Random, Hallmark, and GOBP sets, whereas the Affymetrix table from GSE18842 served as background for the Contextual dataset. In general, either the FDR or the Benjamini-Hochberg (BH) correction method was used, and all ontologies with fewer than 2 genes were removed from the results. Enrichment results for the Hallmark, Contextual, and GOBP datasets are compiled in Supplementary Table S2. All raw outputs are available at [https://github.com/LARA-Lab-Aging]. Additionally, gene-mapping success rates are provided in Supplementary Table S5.

GO term biological specificity was assessed based on the term’s annotation size and its depth in the ontology structure. Both properties have been widely employed to root the evaluation of term specificity (Lewin and Grieve, 2006; Louie et al., 2010; Tomczak et al., 2018). Additionally, the link between them and biological specificity is highly intuitive: a term tends to be more general as more genes are associated with it, and, because the ontology is hierarchical, terms deeper in the hierarchy tend to be more biologically precise. Metrics were retrieved using the GOATOOLS Python library (Klopfenstein et al., 2018).

Semantic similarity (SS) analyses were conducted using GOATOOLS, which implements multiple methods to determine semantic similarity between two GO terms. We selected Wang’s method (Wang et al., 2007) to conduct the SS analyses, as it relies solely on the GO Directed Acyclic Graph structure to define SS and attempts to translate the similarity of two ontologies into biological meaning. The GO’s relationships ‘is_a’ and ‘part_of’ with edge scores of 0.8 and 0.6, respectively, were used to determine Wang’s SS score.

To calculate the defined metrics, we used the enrichment results obtained with the GOBP dataset. True positives (TP) were defined as statistically significant ontologies (adjusted p-value <0.05) that had a Wang’s semantic similarity (SS) score of at least 0.7 with at least one of the target ontologies for that input list. Conversely, false negatives (FN) were ontologies with a Wang’s SS score of at least 0.7, but that were not statistically significant. Similarly, false positives (FP) and true negatives (TN) were ontologies with low maximum semantic similarity (Wang’s SS score <0.3) compared to the original target pathways, which were statistically significant or not, respectively. These Wang’s SS score thresholds were selected to center the analysis on the ontologies that are very similar (Wang’s SS score >0.7) to one of the targets and are, therefore, expected in the results; and on ontologies that are very dissimilar (Wang’s SS score), thus are unexpected.

To group the enriched ontologies into functionally similar clusters, we built interaction networks for all FEA results, in which the edge score between two GO terms was their respective Wang’s SS scores. The networks were then pruned with an edge-weight cutoff of 0.5 to remove edges connecting dissimilar ontologies. The resulting ontology similarity networks were clustered using Markov Clustering (MCL) (Van Dongen, 2008), one of the most robust algorithms for clustering biological data (Brohée and van Helden, 2006; Satuluri et al., 2010; Lim et al., 2019). Of the four Markov Clustering inflation values tested – 1.5, 2.0, 3.5, and 5.0 –, the 5.0 value yielded GO clusters with the highest biological accuracy. Finally, we classified the 15 largest clusters that contained 3 or more ontologies for all the clustered FEA results (Supplementary Tables S3, S4).

Our pipeline is summarized in Figure 1.

To evaluate how tools handle data with seemingly no biological context, we conducted FEA using multiple randomized lists of human protein-coding genes and analyzed the distributions of the number of enriched GO terms across input sizes. In general, when using the nominal p-value to determine significance, tools do yield enrichment results, and the number of enriched ontologies varies considerably across tools (Supplementary Figure S1).

Nevertheless, after adjusting for multiple comparisons, most results for all tested tools are not statistically significant (p-value <0.05) (Figure 2), reinforcing the importance of analyzing enrichment results using corrected p-values, as this dramatically reduces the number of false positives, especially when employing ORA (Hung et al., 2012; Wijesooriya et al., 2022). However, there are tools able to retrieve enriched GO terms despite the nature of the dataset and the statistical correction of p-values (Figure 2). Remarkably, ClueGO and goana consistently yielded the largest number of enriched categories for this dataset (Figure 2). In particular, ClueGO displayed a rather interesting behavior: there are substantially fewer enriched GO terms for the 500-lists in comparison to the other input sizes, despite it retrieving larger results when considering nominal p-values (Supplementary Figure S1). Such behavior could reflect how the p-value correction method is implemented and how detected but not enriched GO terms are handled. Ultimately, this suggests that both ClueGO and goana are susceptible to retrieving enriched GO terms that are either irrelevant or insufficiently descriptive.

These results reiterate that, while ORA-based methods can produce significant noise, this can be mitigated by appropriately correcting p-values. Likewise, it reinforces the need to test the coherence and filtering capabilities of biological information provided by ORA-based tools.

One of the main differences we observed when comparing results across enrichment tools was that the number of enriched GO terms varied significantly in the Hallmark, GOBP, and Contextual datasets (Figure 3). DAVID and topGO were the most conservative tools, retrieving considerably fewer enriched categories than the other tools throughout the three datasets. This could be a double-edged sword: while it identifies fewer enriched categories, potentially highlighting the ontologies closely related to the input, it may also overlook unexpected ontologies that could contribute important biological novelty.

In contrast, ClueGO, GOstats, PANTHER, ShinyGO, WebGestalt, clusterProfiler, g:Profiler, and goana yielded the largest results across the analyzed datasets. In this case, identifying too many enriched GO terms may be detrimental, as they are more prone to false positives and allow the user to choose from an excessively broad range of statistically significant terms, thereby biasing the interpretation of the results toward what is relevant to them. We argue that in such cases, a more stringent approach should be taken to analyze the enrichment results, such as using smaller p-value cutoffs and filtering by term annotation size and overlapping genes.

These results indicate that not all tools may be suitable for all potential research questions, as they retrieve highly heterogeneous sets of ontologies. In this sense, conservative tools might not yield biological novelty, as they are more likely to show expected results. In contrast, tools that retrieve thousands of GOs will not only eclipse relevant and expected results but also force researchers to browse overly extensive lists in search of useful findings.

To estimate the biological informativeness of the FEA results provided by each tool, we used the median and average of both the annotation sizes and depths of the enriched GO terms as measures of biological specificity (Tables 2, 3). We chose to limit our analysis to the top 20 and top 100 enriched ontologies, arranged by p-values in ascending order, to lessen the impact of greatly varying result sizes and to examine the tools’ ranking preferences. Additionally, since the final step of FEA inherently involves manual interpretation, it is reasonable to consider these portions of the results the most relevant. Both selected metrics were analyzed for the top 20 and top 100 ranked enriched ontologies for all FEA results with the Hallmark and Contextual datasets.

Remarkably, DAVID, Enrichr, clusterProfiler, and topGO displayed the smallest annotation sizes and the highest depths across all list sizes in both datasets (Tables 2, 3), indicating that these tools tend to obtain more specific enriched ontologies and provide more descriptive results. For example, their results with the 500-hallmark list include, among the top-ranking enriched GO terms, “nucleotide excision repair” (DAVID and Enrichr), “purine ribonucleotide catabolic process” (clusterProfiler), and “canonical glycolysis” (topGO), which are more informative compared to the broader ontologies “primary metabolic process” and “response to stimulus” (Supplementary Table S2).

Moreover, ClueGO’s output contained highly informative ontologies, especially for smaller list sizes (200, 100, and 50 for the Contextual dataset; 50 and 100 for the Hallmark dataset), but failed to do so with larger input lists. In this sense, as we tested the 500 lists, generic terms such as “metabolic process” and “biosynthetic process” became much more prevalent.

Conversely, g:Profiler and goana produced substantially less informative results, as they retrieved the largest and shallowest GO terms, which tend to be less precise in terms of biological significance. For instance, goana obtained the ontologies “regulation of biological process” and “metabolic process”, among multiple other generic descriptions, across all results with the Contextual and Hallmark lists, respectively. On the other hand, g:Profiler yielded more informative terms than goana, but still included rather broad GO terms, especially with the Contextual dataset, such as “response to stimulus” and “biological regulation”. However, we observed that both tools obtained smaller, deeper ontologies with the 200- and 100-lists in the Hallmark dataset, particularly in the top 20, compared to the top 100.

The remaining tools BiNGO, GOstats, PANTHER, ShinyGO, and WebGestalt provided intermediate results in terms of biological specificity with the Contextual lists. Nonetheless, a pattern similar to g:Profiler and goana with the Hallmark lists of sizes 200 and 100 can be observed. For all of these tools, these particular inputs appear to yield more biologically informative results and place them among the highest-ranked terms. Given that the 500 Hallmark list is relatively more heterogeneous, as it includes genes from three distinct biological hallmarks (DNA repair, Hypoxia, and Response to UV radiation), while the smaller lists are comprised solely of genes from the hypoxia hallmark, our results suggest that these tools tend to lose specificity when analyzing more heterogeneous gene lists. Furthermore, loss in specificity observed with the 50 Hallmark list – which also comprises only hypoxia-related genes – is likely related to the small size of the input. This reinforces our approach by showing that the nature and size of the input data will strongly affect the quality of the FEA results.

Furthermore, by analyzing the annotation size distributions for the top 20, top 100, and all enriched GO terms from each tool, we found that differences in distribution patterns are much more pronounced in the top 20 and top 100 than across all results. Specifically, DAVID, Enrichr, clusterProfiler, and topGO consistently identify smaller ontologies among their top results, suggesting these tools favor more specific GO terms at the top ranks. In contrast, others spread such precise ontologies throughout their results (Figures 4, 5). However, it is essential to note that the output size can influence the comparisons with the full extent of the results, especially in tools that yield larger results. The distribution of GO terms’ annotation size is skewed, with most ontologies having fewer genes annotated to them (Jelier et al., 2008). Consequently, as output size increases, the distribution of annotation sizes in FEA results typically shifts toward smaller values. Also, in the case of ShinyGO, which only outputs up to 1,000 enriched GO terms, distributions can be affected, as not all of the enriched ontologies are accounted for.

The findings above reflect a known issue in the field, widely discussed among users, regarding the relevance of the retrieved GO terms. Users are familiar with the challenge of finding relevant or more descriptive GO terms amidst the vast number of bioprocesses that tools might yield. It is unfeasible or impractical to manually analyze, for instance, more than 5000 GO terms (Figure 2) to identify relevant bioprocesses. Our results indicate that tools differ significantly in the GO terms they rank first and in their level of informativeness.

We also assessed the tools’ performance in terms of accuracy and FPR by using Wang’s SS to construct a confusion matrix for each FEA result with the GOBP dataset. DAVID, WebGestalt, and topGO were the tools that exhibited the highest accuracy values (Table 4). Since accuracy, in the context of this study and of how the confusion matrix was defined, fundamentally denotes how similar and dissimilar the enriched and not enriched GO terms are to the targets for each list, this indicates that the FEA results produced by the selected tools are highly consistent with the biological context of the input list.

Additionally, these tools displayed the lowest FPR values, indicating that the GO terms they enrich are highly related to the input data. ClueGO and ShinyGO were the least accurate software in this analysis (Table 4), primarily due to the generation of numerous false positives (i.e., enriched ontologies unrelated to the context of the input list) (Table 4). Moreover, both tools underperformed in terms of FPR due to the nature of their results, which contained few or none non-statistically significant GO terms, thereby biasing FPR towards 1. This result highlights that performance metrics are highly dependent on the method used to define them and the specificities of the tools. Therefore, despite being undoubtedly relevant performance measures, solely relying on such metrics to compare different tools can be problematic in the case of FEA.

To complement the confusion matrix analysis, we evaluated the identification and ranking of the target ontologies for each list size in the GOBP dataset. Regarding identification abilities, most tools were able to enrich 90% or more of the proposed target (Table 4), with the only exceptions being BiNGO and topGO, which failed to do so (66% and 69% identification rates, respectively). In terms of ranking capabilities, across the 100 and 50 lists, all software ranked the identified target GO terms near the top (Figure 6; Table 4). However, most tools lost this ability and increased the dispersion of the targets in the ranking. DAVID, Enrichr, clusterProfiler, and topGO were the only tools that could consistently position the targets within the top-ranking ontologies across all list sizes (Figure 6), suggesting that they tend to prioritize results closely related to the input.

We further applied MCL to cluster the enriched ontologies for each FEA result with the Hallmark and Contextual datasets. Headers were manually assigned to the 15 largest clusters to compare the biological information provided by each tool. Notably, with the Hallmark lists, most tools identified clusters directly related to the expected biological processes (Supplementary Table S3), namely, DNA repair, hypoxia, and response to UV radiation, since the inputs comprised their respective hallmark genes. The only exceptions were DAVID for the 50-lists and topGO for both 100 and 50-lists, which did not exhibit hypoxia clusters. Both tools did, in fact, enrich at least one hypoxia-related ontology; that enrichment, however, did not cluster with any other enriched GO terms.

As most tools yielded results that were at least partially expected in the Hallmark dataset, the most significant differences ultimately come down to the interpretability of their FEA results. Remarkably, DAVID and topGO, two of the tools that had previously produced the most informative results, also displayed the most coherent and precise clusters of GO terms. This observation aligns with the previously reported high accuracy values. In contrast, tools that generated larger and less specific outputs – for instance, BiNGO, ClueGO, GOstats, PANTHER, ShinyGO, WebGestalt, g:Profiler, and goana – enriched a greater number of terms that were weakly aligned with the biological context of the inputs. Clusters that pertained to “nucleotide metabolism”, “organ and system development”, and “regulation of cell fate and metabolism” were the most frequently produced groupings across all lists in this dataset. Although these are somewhat related to the expected biological processes, their GO terms are not as strongly associated as terms such as “DNA repair”, “response to oxygen levels”, or “response to radiation”. The presence of large numbers of these closely related yet unspecific ontologies might, depending on the tools’ ranking abilities, conceal the most insightful or expected descriptions, thus hindering the interpretability of the results.

Moreover, the behaviors observed in the clustering results for the Contextual dataset FEA analyses were consistent with those observed in the Hallmark dataset (Supplementary Table S4). DAVID still exhibited the most coherent and well-defined clusters for the 500-list, but did not display clustering results for the smaller lists, as it enriched none or only a few GO terms. Similarly, topGO also did not exhibit any clusters for any of the inputs in this dataset for the same reason as DAVID. Nonetheless, in general, the profiles of the biological groups yielded, and, thereby, the overall biological information contained in the enrichment results, were coherent across all tools.

As a consequence, this suggests that the main differences in the biological information contained within the FEA results across all tools lie mainly in their ability to prioritize what’s most relevant in the ranking and to enrich more specific GO terms at the top of the results.