Guidelines for reliable and reproducible functional enrichment analysis

An example gene set analysis

To demonstrate the effect of misusing gene set analysis, a publicly available RNA-seq dataset (accession SRP128998) was downloaded from DEE2 on 13th August 2021 (Ziemann et al, 2019). This data consists of immortalised human hepatocytes cultured in standard or high glucose, first described by Felisbino et al (2021). Transcript level counts were aggregated to genes using the getDEE2 R package v1.2.0. Next, genes with an average of 10 reads per sample were discarded. Differential expression statistical analysis was conducted with DESeq2 v1.32.0 (Love et al, 2014) to identify genes altered by high glucose exposure. For gene set analysis, Reactome gene sets (Jassal et al, 2020) were downloaded (URL: https://reactome.org/download/current/ReactomePathways.gmt.zip , accessed 13th August 2021). FCS was performed using the mitch R package v1.4.0 with default settings, which uses a rank-ANOVA statistical test (Kaspi and Ziemann 2020). Genes with FDR<0.05 were used for ORA analysis using the clusterProfiler R package (v4.0.2) enricher function that implements a hypergeometric test (Yu et al, 2012). No fold-change threshold was used to select genes for ORA. For ORA, two types of background gene sets were used (i) all detected genes, or (ii) all genes in the annotation set. For genes and gene sets, a false discovery rate adjusted p-value (FDR) of 0.05 was considered significant. Analyses were conducted in R version 4.1.0.

Survey of enrichment analysis

We collated 2941 articles in PubMed Central published in 2019 that have keywords “enrichment analysis”, “pathway analysis” or “ontology analysis”. We sampled 200 of these articles and collected the following information from the article, searching the methods sections and other parts of the article including the supplement.

• Journal name

• Type of omics data

• Gene set library used, and whether a version was reported

• Statistical test used

• Whether p-values were corrected for multiple comparisons

• Software package used, and whether a version was reported

• Whether a background gene set was used

• Code availability

• Whether gene profile was provided in the supplement

• Whether the analysis had any major flaws that might invalidate the results. This includes:

1. background gene set not stated or inappropriate background set used,
2. lack of FDR correction,
3. no enrichment data shown,
4. inference without performing any statistical test, and
5. misinterpreting p-values by stating results were significant when FDR values indicate they aren’t.

We excluded articles describing novel enrichment analysis techniques/tools, review articles and conference abstracts. Some articles presented the results of >1 enrichment analysis, so additional rows were added to accommodate them. These data were entered into a Google Spreadsheet by a team of five researchers. These articles were cross checked by another team member and any discrepancies were resolved.

For analyses using GSEA, we scanned the articles to identify whether key methodological steps were described, including (i) the gene weighting parameter, (ii) test type, ie: permuting sample labels or genes, and (iii) method used for ranking genes.

For ORA of differential gene expression analyses, we inspected the articles to see evidence that researchers were conducting enrichment of up and down regulated genes separately or combined.

Cleaned data were then loaded into R v4.1 for downstream analysis.

For assessment of enrichment analysis quality with journal metrics and citations, we required a larger set of studies so we sampled and analysed a further 1300 articles from PMC. Results from this sample were not double-checked and so may contain a small number of inaccuracies. We rated each analysis with a simple approach that deducted points for missing methodological details and awarding points for including extra information (Table 1).

Scimago Journal Rank (SJR) data were downloaded from the Scimago website (https://www.scimagojr.com/journalrank.php) and used to score journals by their citation metrics. Using NCBI’s Eutils API, we collected the number of citations each article accrued since publication. We used Pearson correlation tests to assess the association with the analysis scores we generated. All data analyses were conducted in R v4.1 and the analysis scripts are available at GitHub (https://github.com/markziemann/SurveyEnrichmentMethods).

An example of gene set analysis misuse

In order to demonstrate the effect of gene set test misuse, we used an example RNA-seq dataset examining the effect of high glucose exposure on hepatocytes. Out of 39,297 genes in the annotation set, 15,635 were above the detection threshold (>10 reads per sample on average). Statistical analysis revealed 3,472 differentially expressed genes with 1,560 up-regulated and 1,912 down-regulated (FDR<0.05) due to glucose exposure.

FCS revealed 95 up and 316 downregulated pathways (FDR<0.05). ORA with a background gene list consisting of detected genes revealed 55 upregulated and 238 downregulated gene sets (FDR<0.05). There was a strong overlap between FCS and ORA results with a Jaccard statistic of 0.66 (Figure 1A).

We then performed ORA using the whole genome gene list (indicated as ORA*), which resulted in 147 up and 530 down regulated gene sets (FDR<0.05), The overlap of ORA with ORA* was relatively smaller, with a Jaccard statistic of 0.41 (Figure 1B). Interestingly, 26 gene sets were simultaneously up and downregulated with this approach.

Then we performed ORA analysis after combining up and downregulated genes as sometimes observed in the literature (indicated as ORA comb) with the standard background and again with a whole genome background (indicated as ORA* comb). ORA comb revealed only 22 differentially regulated gene sets with a Jaccard similarity to ORA of only 0.04 (Figure 1C). ORA* comb yielded 64 differentially regulated gene sets with a Jaccard similarity to ORA of only 0.08. These results demonstrate that commonly observed practices of combining the up and downregulated gene lists and using the whole genome background gene list severely distort ORA results.

Figure 1. Example of enrichment analysis misuse. (A) Comparison of FCS and ORA methods. (B) ORA with standard and whole genome background. (C) Effect of combining up and downregulated genes on ORA results.

A screen of functional enrichment analyses in open access journals

A search of PubMed Central showed 2941 articles published in 2019 with the keywords “enrichment analysis”, “pathway analysis” or “ontology analysis”. From these, 200 articles were screened for methodological errors. We excluded 14 articles from the screen because they did not present any enrichment analysis. Those excluded articles included articles describing novel enrichment analysis techniques or tools, review articles or conference abstracts. As some articles included more than one enrichment analysis, the dataset included 235 analyses from 186 articles; this data is available in Table S1. A flow diagram of the survey is provided in Figure 2.

Figure 2. A summary of the survey of functional enrichment analyses.

There were articles from 96 journals in the sample, with PLoS One, Scientific Reports, and PeerJ being the biggest contributors (Figure S1A). There were 18 different omics types, with gene expression array and RNA-seq being the most popular (Figure S1B). There were 31 different species under study, but human was the most common with 157 analyses (Figure S1C).

We recorded the use of 26 different gene set libraries, with GO and KEGG being the most frequently used (Figure 3A). There were 14 analyses where the gene set libraries used were not defined in the article. Only 18 analyses reported the version of the gene set library used (Figure 3B). There were 12 different statistical tests used, and the most common reported tests were Fisher’s Exact, GSEA and hypergeometric tests; but the statistical test used was not reported for the majority of analyses (Figure 3C). Out of the 225 analyses that performed a statistical test, only 119 (53%) described correction of p-values for multiple testing (Figure 3D).

There were 50 different tools used to perform enrichment analysis, with DAVID and GSEA being the most common; while 15 analyses (6.4%) did not state what tool was used (Figure 3E). The version of the software used was provided in only 68 of 235 analyses (29%) (Figure 3F).

For analyses using ORA methods, we studied what background gene set was used (Figure 3G). This revealed that in most cases, the background list was not stated or it was clear from the article methods section that no background list was used. In a few cases, background list was mentioned but was inappropriate, for example using a whole genome background for an assay like RNA-seq. In only 8/197 of cases (4.1%), the appropriate background list was described in the article.

Of the 47 analyses which used a scripted computer language, only 3 provided links to the code used for enrichment analysis (6.4%) (Figure 3H). For 93 of 235 analyses (40%), the corresponding gene lists/profiles were provided either in the supplement or in the article itself (Figure 3I).

Next, we quantified the frequency of methodological issues and reporting that would undermine the conclusions (Figure 3J). Lack of appropriate background was the most common issue (179 cases), followed by FDR (94 cases), then lack of data shown (13), inference without test (11 cases), and misinterpreted FDR values (2 cases). Only 35 analyses (15%) did not exhibit any of these major methodological issues.

During the course of this survey, we also noticed some studies that grouped up- and down-regulated gene lists together prior to ORA, a practice we were not expecting. To assess this in a systematic way we assessed differential gene expression studies to determine how common it was to combine up- and down-regulated gene lists prior to ORA (Figure 3K). From 109 analyses, just 28 studies performed separate analysis of up- and down-regulated gene sets.

We also looked at studies performing GSEA, and whether three important analytical choices were described in the methods. These are (i) the gene weighting parameter, (ii) test type, ie: permuting sample labels or genes, and (iii) method used for ranking genes. These parameters were not stated in more than half of the analyses (Figure 3L).

Taken together, these data suggest methodological and reporting deficiencies are widespread in published functional enrichment analyses.

Figure S1. General make up of the sample. (A) Journals. (B) Types of omics analysis applied. (C) Organisms under study.

Figure 3. Findings from the survey of published enrichment analyses. (A) Representation of different gene set libraries. (B) Proportion of analyses reporting gene set library version information. (C) Representation of different statistical tests. (D) Proportion of analysis conducting correction of p-values for multiple comparisons. (E) Representation of different software tools for enrichment analysis. (F) Proportion of analyses that reported version information for software used. (G) Background gene set usage and reporting. (H) Proportion of scripted analyses that provided software code. (I) Proportion of analyses that provided gene profiles. (J) Proportion of analyses with different methodological problems. (K) Whether up- and down-regulated genes were analysed separately or combined. Limited to ORA of differential expression only. (L) Reporting of GSEA parameters.

Analysis scores

In order to make an accurate assessment of association between analysis quality and bibliometrics, we required a larger sample of articles. We therefore analysed a further 1300 articles, bringing the total number of analyses described up to 1630; this dataset is available in Table S2.

We then scored each analysis based on the presence or absence of methodological issues and included details. The median score was -4, with a mean of -3.5 and standard deviation of 1.4 (Figure 4A). Next, we assessed whether these analysis scores associate with journal citation metrics (Figure 4B). There was no statistically significant association of analysis score and SJR (Pearson r=0.03, p=0.23).

Next we wanted to know which journals had the highest and lowest scores. Only journals with five or more analyses were included. The best scoring journals were Transl Psychiatry, Metabolites and J Exp Clin Cancer Res (Figure 4C), while the poorest were RNA Biol, Mol Med Rep, Cancer Med and World J Gastroenterol (Figure 4D), although we note that there was a wide variation between articles of the same journal.

Then we assessed for an association between mean analysis score and the SJR, for journals with 5 or more analyses (Figure 4E). Again, there was no association between analysis score and SJR (Pearson r=-0.008, p=0.95). Next we assessed whether there was any association between analysis scores and the number of citations received by the articles. After log transforming the citation data, there was no association between citations and analysis scores (Pearson r=0.002, p=0.95) (Figure 4F). Taken together, these findings suggest that these methodological issues are not limited to lower ranking or poorly cited articles.

Figure 4. Comparison of Analysis scores to bibliometrics. (A) Distribution of analysis scores. (B) Association of analysis scores with Scimago Journal Rank (SJR). (C) Journals with highest analysis scores. (D) Journals with lowest analysis scores. (E) Association of mean analysis score with SJR for journals with 5 or more analyses. (F) Association of analysis scores with accrued citations.

Essential minimum standards for functional enrichment analysis

The following guidelines are selected to avoid the most severe methodological issues while requiring minimal effort to address.

1. Report the origin of the gene sets and the version used or date downloaded. These databases are regularly upgraded so this is important for reproducibility.

2. Report the tool used and its software version. As these tools are regularly updated, this will aid reproducibility.

3. If making claims about the regulation of a gene set, then results of a statistical test must be shown (including measures of significance and effect size). The number of genes belonging to a pathway on its own is insufficient to infer enrichment. “Data not shown” is not acceptable when it is possible to attach supplementary files.

4. Report which statistical test was used. This will help long term reproducibility, for example if in future the tool is no longer available. This is especially relevant for tools that report the results of different statistical tests.

5. Always report FDR or q-values (p-values that have been corrected for multiple tests). This will reduce the chance of false positives when performing many tests simultaneously (Hung et al, 2012). Report what method was used for p-value adjustment. Bonferonni, FDR and FWER are examples.

6. If ORA is used, it is vital to use a background list consisting of all detected genes, and report this in the methods section. Avoid tools that don’t accept a background list.

7. Report selection criteria and parameters. If performing ORA, then the inclusion thresholds need to be mentioned. If using GSEA or another FCS tool, parameters around ranking, gene weighting and type of test need to be disclosed (eg: permuting sample labels or gene labels). Any parameters that are vary from the defaults should be reported.

8. If using ORA for gene expression analysis, be sure to conduct ORA tests separately for up- and down-regulated genes before analysis. Combining up- and down-regulated genes into a single list will limit sensitivity.

9. Functional enrichment analysis results shouldn’t be considered proof of biological plausibility nor validity of omics data analysis (Timmons et al, 2015). Such tools are best used to generate new hypotheses; informing subsequent biochemical or signaling experimentation.

Additional recommendations for enhanced reproducibility

The following guidelines are directed at going over and above the minimum standards towards a gold standard of reliability and reproducibility.

10. Preference FCS and pathway topology tools over ORA tools. ORA tools have a lower sensitivity. In practice, they are rarely used properly.

11. Include the gene profile data (including any background lists) in the supplementary data in TSV or CSV formats. Excel spreadsheets are not recommended for genomics work (Ziemann et al, 2016)

12. Scripted analysis workflows are preferred over analysis involving graphical user interfaces because they can attain a higher degree of computational reproducibility, so consider investing in upgrading your analysis. Similarly, free and open source software tools are preferred over proprietary software.

13. If using a scripted analysis, provide code by depositing it to a repository like GitHub or Zenodo. Link the data to the code, so anyone can download and reproduce your analysis (Peng 2011).