Methods

From a list of 186 PubMed central articles using enrichment analysis that we examined previously [7], we selected 20 articles that fulfilled the following criteria; (1) the organism under study is human; (2) the tool used was DAVID (Database for Annotation, Visualization, and Integrated Discovery) [10] ; (3) the omics type was gene expression (array or RNA-seq); and (4) that the gene list was provided in the supplementary data. The basis for these selection criteria is that this combination represents the most common application of enrichment analysis according to our previous study, which could be reproduced without the need to reanalyse raw data. These articles are listed in the bibliography [11–30].

We screened reporting criteria of these articles, including whether the version of the application is given, whether p-values underwent correction for multiple testing, whether a background list was described and whether gene selection was described unambiguously.

The lists of genes were collected from the supplementary files and stored in a spreadsheet, paying special attention to preventing gene name errors from occurring [31]. The reproducibility analysis was performed using the respective methodologies stated in the methods sections of the 20 articles using DAVID v6.8, as at this time v6.7 was not available (conducted March-May 2022) [4].

Briefly, we used the DAVID v6.8 “Analysis Wizard” and pasted in the provided gene list into the text box, selecting this as “gene list,” together with the appropriate gene identifier type. Typically, this was the official gene symbol. The species was specified as Homo sapiens. Next, we clicked on the “Functional Annotation Tool” link, bringing up the “Annotation Summary Results” page. We followed the steps indicated in the respective papers to select appropriate gene sets investigated in the articles such as Gene Ontology (GO; biological process, molecular function), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, etc. Next, the “Functional Annotation Chart” was clicked on, opening a new window with a table of enrichment results, which were saved for future reference.

The pathway enrichment results presented in the paper were compared to those we obtained. It would be ideal to compare standardised effect sizes (fold enrichment scores) between original and reproduction results, however these were reported in the main text of only 3/20 articles.

Therefore we adopted a binary classification system, based on head-counting pathways that were presented as statistically significant in the original article. We defined a pathway as confirmed when it reached a similar level of statistical significance upon reproduction. We define “similar significance level” to be a value within 50% of the value presented in the original article. For example if an article shows a pathway with FDR=0.04, then an acceptable FDR in reproduction could be in the range of 0.02 to 0.06. Additionally, statements made in the abstract, results, and discussion sections based on specific pathway enrichment results were carefully examined and cross-checked with our results to understand whether those assertions survived reproduction or not.

Each article was critically evaluated and given a score of one (1; low) to three (3; high). Articles were rated 1 when <50% of pathways shown in the results were confirmed upon reproduction, and most conclusions derived from enrichment analysis were not supported by reproduction. Articles were rated 2 when ≥50% of pathways shown in the results were confirmed upon reproduction, and most conclusions derived from the pathway enrichment analysis were supported. Articles were rated 3 when ≥80% of terms shown in the original article were observed again in reproduction with a similar significance value, and all conclusions derived from the pathway analysis were supported in reproduction.

Results

As reproducibility is closely associated with the quality of documentation, we collected key methodological information about how the enrichment analysis was done and present in Table 1. From this, we can see that only seven articles reported the version of DAVID that was used. Version 6.7 was reportedly used in two articles while version 6.8 was specified in five. Correction of p-values for multiple testing was specified or inferred in 11 articles, while the results of the remaining nine articles did not take into consideration the importance of p-value correction for minimising false positives. Of the 20 articles, none mentioned explicitly that a custom background list was used for enrichment analysis. We also examined whether the gene selection criteria were unambiguously clear, which was true in 16 articles.

We attempted reproduction according to the methods described in the respective articles, with the limitation that we had to use DAVID v6.8 as v6.7 was no longer available. The 20 articles were classified based on the degree of reproducibility (Figure 1). Seven studies were classified as low reproducibility, while nine were classified as medium and just four were classified as high. Justifications for these classifications are provided in Table S1.

Figure 1. Examining reproducibility of enrichment analysis in a sample of 20 journal articles. Beside the PubMed Central identifier, the reported version of DAVID used in the original publication (blank if none given).

## png 
##   2

The ability to reproduce DAVID enrichment analysis using the supplementary gene lists was highly variable. The version of DAVID used is likely to be a big contributor to irreproducibility as we know that there are major differences in the underlying annotation set (knowledgebase) used by version 6.7 and 6.8. This is supported by the fact that none of the v6.7 studies received a high reproduicibility score. On the other hand these versions did not receive regular updates, so the knowledgebases were static for relatively long times; v6.7 (Jan. 2010 - Oct 2020) and v6.8 (Oct 2016 - Jun 2022). We can then deduce that knowledgebase version is likely not the cause of irreproducibility for PMC6561911, PMC6539328 and PMC6425008; studies that reportedly used v6.8 but received low reproducibility scores.

We note that DAVID version 6.8 was retired in June 2022, so these 20 studies are currently no longer reproducible with the original tools. The implications of this are outlined in the Discussion section below.

Method overview

We have developed a step-by-step protocol for how to conduct a highly reproducible enrichment analysis, which is attached as a supplementary file (Suppl Info 1) and published to protocols.io with a permissive license [36]. This protocol is prepared in accordance with the five pillars principles [35], and includes version control with git and GitHub, compute environment control using Docker, literate programming with R Markdown, persistent code and data sharing with Zenodo and documentation of all of these to empower others to reproduce the workflow. Central to this protocol is the provision of template Dockerfile and an R Markdown script, which are intended to be remixed and altered by end users to suit their project’s needs. The general steps involved are shown in Figure 2.

Figure 2. Flow chart of example workflow. An estimated amount of time to conduct each step is given.

## png 
##   2

The Dockerfile uses a base image from Bioconductor that corresponds to the latest available stable release version, on the Ubuntu operating system. The image contains instructions for installing a few useful utilities including nano and git for modifying scripts, and magic-wormhole for transferring data between computers. It installs R packages from CRAN and Bioconductor repositories to fulfill common tasks in transcriptome analysis such as fetching RNA-seq counts from DEE2 with getDEE2 [37], differential expression with DESeq2 [38], enrichment analysis with clusterProfiler, fgsea and mitch [39–41]. It also contains instructions to clone a GitHub repository, which contains the scripts that run the analytical workflow. Lastly, it copies a gene set library from Reactome [42] with a clear version history to guarantee provenance of the annotation data.

The R Markdown script called example.Rmd undertakes the transcriptome and enrichment analyses. R Markdown allows the unification of extended descriptions, code and results into a single document, which is becoming ubiquitous in bioinformatics and data science more generally [43]. R Markdown supports languages other than R, including Python, Shell, Julia and others. The ability to provide extended descriptions is helpful for providing contextualising information about the purpose of the work, detailed methods and interpretations derived from the analysis. The R Markdown script contains code chunks that perform specific tasks including fetching an example data set from DEE2, inspecting data quality, differential expression and two types of enrichment analysis: over-representation analysis (ORA) and functional class scoring (FCS) (Figure 3).

For each enrichment analysis conducted, a checklist is provided to comply with reporting guidelines [8,44]. An example checklist is shown in Table 2.

The software versions used in this example analysis are provided in Table 3.

Reproducing the example workflow is relatively quick and simple, requiring six commands and approximately 4-7 minutes starting from a base Linux install. The steps include installing Docker, pulling and running the container, running the R Markdown script, exiting the container and copy the newly created report from the container to the current working directory (Box 1).

sudo apt update && sudo apt install docker.io -y  # 30 sec
sudo docker run --platform linux/amd64 mziemann/enrichment_recipe
Rscript -e 'rmarkdown::render("example.Rmd")' # 150 sec
docker cp $(docker ps -aql):/enrichment_recipe/example.html . # instant
firefox example.html

Box 1. Linux shell commands to reproduce a transcriptome and enrichment analysis using a Docker image. Ubuntu 22.04 was used.

Demonstration

To demonstrate the utility of this workflow, we present the results of this example data analysis; it is a comparison of control and 5-azacitidine-treated AML3 cell transcriptomes with RNA-seq in triplicate [45]. We use this analysis to better understand the effect of DNA methyltransferase inhibitors on gene expression and elaborate upon previous findings [45], which mainly focused on the relationship between RNA expression and DNA methylation but did not show results of a transcriptome-only enrichment analysis. Furthermore, we provide a comparison of the two most common approaches for enrichment analysis; ORA and FCS, and give guidance for interpreting results of such tests.

For quality control the number of reads per sample was calculated and found to be in the range of 11-16M. Genes with fewer than 10 reads per sample on average were removed from the analysis, leaving 13,168 genes. A multi-dimensional scaling plot showed clustering of samples into their treatment groups, suggesting a robust effect of the chemical exposure (Figure 4A). Differential expression analysis revealed 3,598 differentially expressed genes with FDR<0.05, with 1,672 exhibiting higher expression and 1,926 exhibiting lower expression in 5-azacytidine treated cells (Figure 4B). Up and down-regulated gene sets were subjected to ORA separately, using the list of 13,168 genes as the background, and gene sets from REACTOME.

Figure 4. Differential expression and enrichment analysis with the example workflow. (A) An MDS plot of gene expression data for the largest two principal components. Blue and pink points represent control and azacitidine-treated samples respectively. (B) An MA plot. The x-axis represents the log₁₀ base mean expression, and the y-axis represents the log₂ fold change. Points in red indicate genes with statistically significant differential expression (FDR<0.05), while those in grey did not meet this threshold. (C) Bar plot of ORA-based enrichment test results. (D) Bar plot of FCS-based enrichment test results. (C,D) REACTOME gene sets with FDR>0.05 were discarded and the remaining gene sets were ranked by enrichment score (ES) The most extreme up- and down-regulated genes sets are shown (based on ES). Red bars indicate up-regulation and blue bars indicate down-regulation. (E) Euler diagrams of differentially regulated gene sets identified with ORA and FCS approaches.

With ORA, we observed 201 up-regulated and 99 down-regulated gene sets (FDR<0.05). We prioritised enrichment results by enrichment score, a surrogate measure for effect size, which tends to highlight smaller and more specific gene sets undergoing more drastic differential expression, rather than large gene sets that undergo small expression changes. This led to the identification of up-regulated pathways related to nuclear RNA import and glucokinase regulation; while down-regulated pathways were associated with interferon signaling, interleukin signaling and NADPH oxidase (host defence) (Figure 4C).

After FCS, we identified 239 up-regulated gene sets and 158 down-regulated gene sets (FDR<0.05). FCS up-regulated pathways were associated with prometaphase chromosome condensation, nuclear RNA transport and defective homologous recombination repair (Figure 4D). This finding is consistent with changes to the cell cycle observed in 5-azacytidine treated cells, where the proportion of cells in G2/M phase is reduced, and proxy measures of apoptosis appear slightly higher [45].

Down-regulated pathways related to endosomal/vacuolar (antigen presentation) pathway, interferon signaling and TLR/TNF signaling, which was not described in the original publication [45]. Interestingly, these pathways are reported to be upregulated in other cancer cell studies of DNA hypomethylating agents for similar durations [46–48].

Regarding up-regulated pathways, the similarity between ORA and FCS was relatively high, exhibiting a Jaccard index of 0.58, while this agreement was lower for down-regulated pathways, with a Jaccard index of 0.46 (Figure 4E).