gms | German Medical Science

In recent years, the improvements of high-throughput technologies for Omics experiments produce more data sets for analysis [1], [2], [3]. Consequently, the challenge is to interpret these data for a better understanding of diseases or a given phenotype [1], [4], [5], [6]. The most common outputs are gene expression profiles which are compared to examine the changes in gene expression between different test groups [5]. A standard analysis method on such data is the statistical enrichment with pathways, molecular functions, or biological processes [7]. There are various algorithms to compute enrichments which can be divided into over-representation analysis (ORA) [8], functional class sorting, and topology-based methods [9]. Unfortunately, no standard method or preferred category exists as such but ORA and functional class sorting methods are utilized in many publications [1], [2], [9].

Currently, a lot of different pathway databases exist like Kyoto Encyclopedia of Genes and Genomes (KEGG) [10], Reactome [11], WikiPathways [12], and PathBank [13]. Every pathway database has its focus and applications. The choice of pathway databases influences the results of enrichment [14]. That is the reason why recently the usage of multiple pathway databases or a merged pathway database became more prevalent [14], [15], [16], [17].

Here, a new tool of PharMeBINet’s web application is introduced allowing users to perform pathway enrichment analysis and investigate results in multiple ways. In PharMeBINet multiple pathway databases are connected and used for enrichment. For the enrichment analysis, two methods of the Fisher’s exact test are used, the standard Fisher’s exact test and a modified version of the Fisher’s exact test, which considers the ranking of the gene set. Four correction approaches for multiple testing are available to the user controlling the rate of false positives. The results are visualized as bar plots, heat map, and pie chart. The tool is available at https://pharmebi.net/#/analysis/enrichmentanalysis.

In the following, the development of the website is described as well as the databases utilized for the enrichment methods. The implemented enrichment and correction methods are explained. Afterward, a short explanation of how to use the website and analysis module is described as well as the use-case for comparison.

The web application is divided into frontend and backend. Vue.js 2.6.11 [18], an open-source JavaScript framework for user interfaces, generates the frontend. For the design of Vue, the library Vuetify 2.6.2 [19] is used intending to create intuitive and reliable user interfaces based on material design. Echarts 5.3.2 [20] generates different graphics and charts on the webpage. Node.js v14 LTS [21] is an open-source JavaScript framework which forms the backend. It communicates with the Neo4j database utilizing the library Neode 0.4.8 [22], a Neo4j object graph mapping (OGM) which takes care of creating, reading, updating, and deleting (CRUD) operations on data in the database. The communication between frontend and backend is implemented using Axios 0.27.2 [23], an open-source JavaScript library for HTTP requests.

The backend utilized the heterogeneous Neo4j database PharMeBINet [24]. Pathway information is derived from Comparative Toxicogenomics Database (CTD) (version 2022-04) [25], Pathway Commons (version 12) [26], Reactome (version 2022-03-31) [11], and WikiPathways (version 2022-04-10) [12].

The databases Pathway Commons and WikiPathways provide gene-pathway relationships. However, Pathway Commons is filtered for information with compatible licenses before integration. Reactome’s pathway information is mapped via Reactome identifier and name. Reactome pathways that do not map to Pathway Commons or WikiPathways are added as new nodes. CTD is mainly used to fill up additional gene-pathway relationships. Gene information is derived from Entrez Gene [27]. The resulting database contains 3,689 different pathway nodes extracted from four databases.

Fisher’s exact test is a common ORA method utilized for the pathway enrichment analysis [1]. A list of candidate genes is tested for over-representation in a list of genes connected to a given pathway. Fisher’s exact test calculates the statistical significance between the genes of interest and the genes of a pathway. The null hypothesis for each pathway is that a special set of genes is not over-represented in comparison to all genes of the pathway. For this, a contingency table (Table 1 [Tab. 1]) is utilized where N is the number of genes in the input list, n is the number of genes connected to the pathway, and M is the number of all genes.

The p-value is computed by the hypergeometric distribution:

For the right-tailed test, the p-value is computed with Pr(x≤k). The input list is significantly over-represented in a given pathway if the p-value falls below a threshold, typically set to p<0.05. The repeating calculation of Fisher’s exact test on different pathways increases the possibility for false-positive results [28], [29]. For this problem, multiple methods to control the type I error have been developed like Bonferroni correction (BC), Dunn-Šidák correction (DS-correction), Holm-Bonferroni method (HB-correction), and Benjamini-Hochberg procedure (FDR-correction) [28], [29]. Dunn-Šidák and Bonferroni correction adjust the threshold with a number of tests [29]. In contrast, for the Benjamini-Hochberg procedure and the Holm-Bonferroni method the p-values are ordered and assigned a rank. Consequently, the threshold is changed by the number of tests and depending on their rank [28], [29].

Usually Fisher’s exact test does not take the ranking of genes into account [15]. However, g:Profiler [15] developed an algorithm applying Fisher’s exact test on an ordered gene list. First, the gene list is ordered, for example, by fold-change. For each prefix sublist of the ordered gene list with length ≤N the algorithm computes the intersection x to the genes of a pathway as well as the p-value. Afterward, the minimal p-value and respective prefix are returned.

Additionally, the size of expected overlap by chance is calculated by considering the amount of genes from a given pathway compared to the background genes and multiplied by the number of genes in the input list. The fold-enrichment value is computed by dividing the overlap of genes between the pathway and the input list through the expected overlap by chance.

Not only genes can be utilized as input in the new PharMeBINet analysis module, but also diseases and chemicals (see Figure 1 [Fig. 1]). The tool can perform enrichment analyses for pathways, diseases, chemicals, biological processes, cellular components, molecular functions, and genes. First, the threshold for rejecting the null hypothesis can be set manually. Secondly, it can be selected if the threshold is used on the raw p-value or the corrected p-value. Thirdly, four different correction methods can be selected. It is possible to use multiple corrections at the same time. A pathway passes the threshold if at least one of the selected correction methods is below the threshold. Lastly, it is possible to select the background value from the database or to define a custom value. The background from the database is defined as the number of nodes from the input label connected to the enrichment label.

The right side of Figure 1 [Fig. 1] shows the possibilities for input data. Genes, diseases, and chemicals can either be put in using a text field or by uploading a tab-separated values (TSV) file with or without a header. For both methods, the first column represents an identifier field and the second column has some kind of weight. For genes, these weights could, for example, represent the fold-change.

Shin et al. [30] took 63 samples from 62 patients with cutaneous T-cell lymphomas (CTCL) to analyze the gene expression profile. CTCL represents a subset of the non-Hodgkin lymphomas, with T-cell lymphomas present in the skin [31], [32]. They analyzed samples from patients with cancer stages IA, IB, IIB, and III. The expression results are clustered using hierarchical clustering and self-organizing maps, from which three unique clusters are extracted. A unique selling point is that the first cluster contains patients from all stages of cancer. Additionally, the authors performed pathway enrichment on the first cluster. This is why upregulated genes of this cluster 1 (see Table 2 [Tab. 2]) are chosen for testing the pathway enrichment analysis presented here.

For the given gene expression data of upregulated genes of cluster 1, PharMeBINet finds 68 different pathways (Attachment 1 [Attach. 1]) of which the top fifteen are shown in Table 3 [Tab. 3]. The most significant pathway is “Generation of second messenger molecules” with an adjusted p-value of 2.858e-13 and a fold-enrichment of 57.25. Conditions of the previous study are applied to this data set in which only pathways are considered which contain at least 5 genes [30], reducing the result to 45 different pathways.

The enrichment with the modified Fisher’s exact test for ranked genes returns 104 different pathways shown in Attachment 2 [Attach. 2]. The most significant pathway is “Generation of second messenger molecules” too, but the adjusted p-value is better with 9.526e-14 and fold-enrichment of 61.57. The better values are because only the first 53 inputs were considered. Additionally, the top fifteen of the Fisher’s exact and the modified version are the same but some are ranked differently. Moreover, with the condition to have at least five genes in overlap with the input list only 57 pathways are left.

In addition to the tabular demonstration of the pathway enrichment, there is a graphical visualization. A bar plot (Figure 2 [Fig. 2] and Attachment 3 [Attach. 3]) shows the –log10(p-value) for the raw p-values and the different corrected (BC, FDR-, DS-, HB-correction) p-values for each pathway. It is possible to zoom in and out with the slide bar. The bar plot demonstrates that the BC is the strictest correction. In contrast, the FDR-correction is the least strict. However, the p-value is strongly reduced by each correction method. The FDR-correction demonstrates that with increasing rank the distance to raw p-values is decreased. Also, the black dashed line defines the threshold. The last ten pathways are only considered due to the FDR-correction.

Another visualization is a combination of a heat map and bar plot as demonstrated in Figure 3 [Fig. 3] and Attachment 4 [Attach. 4]. Each line represents a pathway and the heat map highlights the gene from the input list which is connected to this pathway. The bar shows the –log10(p-value) of the pathway. A slide bar makes it possible to zoom in and out of the heat map. The pathway “Immune System” includes the most genes from the input list but does not have the best p-value.

Further, a bar plot visualizes the input genes and their count of enriched pathways as seen in Figure 4 [Fig. 4] and Attachment 5 [Attach. 5]. It demonstrates that some genes do not appear in any of the enriched pathways. In contrast, the gene LCK appears in 37 and in the modified Fisher’s exact test in 52 different enriched pathways.

Another statistic is shown in Figure 5 [Fig. 5], demonstrating how many genes are in the input and how many have at least one connection to any pathway in PharMeBINet. The input lists consist of 62 different genes for this example and 5 have no connection to any pathway connection. The genes without any connection are POU2AF1, C1orf68, PSORS1C2, TRA@, and EVI2B.

In the following, the results are compared and discussed. First, the results of the PharMeBINet enrichment analysis for the use-case are compared with the published results as well as the resulting pathways validated for CTCL. The results are compared between the standard and the modified Fisher’s exact test and also between the different correction methods. Afterward, the results are compared with results from other available enrichment tools as well as their provided result visualizations.

PharMeBINet’s pathway enrichment analysis finds 68 relevant pathways for gene cluster 1. The 15 best hits mainly fit the “T cell receptor signaling pathway” (rank 5, 6, 10, and 13) and “T cell activation” (rank 11 and 12) comparable to the analysis of Shin et al. [30]. According to these pathways, “generation of second messenger molecules” (rank 1) is detected which is a specific part of the TCR-signaling in H. sapiens. Not only specific pathways are found but also the “immune system pathway” (rank 4) and “inflammatory response pathway” (rank 19), which represent both more general pathways and both contribute to the phenotype of CTCL. After TCR activation, naive CD4- T-cells differentiate into one of various T helper cells like Th17 cells [33]. In the Th17 cell differentiation pathway (rank 2), Th17 cells are required for immune responses to various extracellular bacteria and fungi. In CTCL, mycosis fungoides is the most common type. So, this suits genes that are highly expressed in this tissue. These cells produce Interleukin (IL)-17 among others, which describes a pathway Shin et al. found in their pathway-enrichment analysis. Not only similar results are found in terms of content, but also much more relevant pathways than in the analysis of Shin et al., which is due to the higher number and specificity of pathways. In PharMeBINet 3,689 pathways are included from three databases, whereas Shin et al. considered six databases with 408 pathways.

Some of the enriched pathways are opposites of each other like “PD-1 signaling” and “Cancer immunotherapy by PD-1 blockade’”. The PD-1 blockade pathway is the pharmacological treatment of cancer by immunotherapeutic inhibition (such as Cemiplimab) of reactions in the “PD-1 signaling” pathway.

Multiple enriched pathways are part of inflammatory processes of the immune system. As T-cells play an important role in the immune system response, the inclusion of pathways involving T-cells like “T-cell activation SARS-CoV-2” and pathways involving T-cell receptor CD4 like “Binding and entry of HIV virion” in the results are to be expected.[Fig. 5]

The three best pathways of both results are equal, however, the p-values, fold enrichments and the numbers of nodes expected by chance are different. This difference occurs due to the use of a prefix sublist in the modified test changing the value for N. However, the top 15 pathways of both enrichment tests are the same although their order, p-values, the number of nodes expected by chance, and the fold enrichments are different (see Table 3 [Tab. 3] and Attachment 2 [Attach. 2]). This is caused by the different prefixes of the input list by the modified Fisher’s exact test to get the best p-value and results in some genes of the input list not being considered even though they appear in the pathway. An example of this is the pathway “Innate Immune System” PC12 4895. In Fisher’s exact test it has an overlap of 15 genes but in the modified version the overlap is only 14. In the modified version only the first 47 genes are considered so the gene PTPRC from position 54 is not considered (see Attachment 1 [Attach. 1], Attachment 2 [Attach. 2], Attachment 4 [Attach. 4], and Figure 3 [Fig. 3]).

The pathway “Wnt signalling” is only part of the modified enrichment result. In the paper Shin et al. the association between the input gene and the pathway is described. This underlines that the modified Fisher’s exact test returns additional useful pathways.

The number of results for Fisher’s exact test and modified Fisher’s exact test would be nearly the same if the FDR-correction is not used. The number would be 40 for Fisher’s exact test and 43 for the modified test. This demonstrates that the other corrections reduce the number of results drastically for both enrichment methods.

For the enrichment analysis, four different correction methods are implemented: BC, DS, HB, and FDR. BC is the strictest correction of the four (see Figure 2 [Fig. 2] and Attachment 3 [Attach. 3]). This is caused by the fact that the family-wise error is below the given threshold [34]. In contrast, it increases the type II error filtering out significant pathways. The DS- and HB-correction are less strict than BC; however, for this example data set all three surpass the threshold for the same pathways (see Attachment 1 [Attach. 1], Attachment 2 [Attach. 2], and Attachment 3 [Attach. 3]). The loosest correction is the FDR-correction which results in 28 additional pathways for the Fisher’s exact test and 61 for the modified test.

“Measles virus infection” is the first pathway that is only included because of the FDR-correction. The paper Künzi et al. [35] shows that the use of the measles virus is a possibility to treat cutaneous T-cell lymphoma. This demonstrates that the pathway “Measles virus infection” is associated with cutaneous T-cell lymphoma and associated with the input list.

“T cell receptor and co-stimulatory signalling” is the first pathway of the modified Fisher’s exact test result where only the FDR-correction value is under the threshold. This pathway is similar to the pathway “T cell receptor signaling pathway”, which belongs to the top enriched pathways and is a good significant pathway for the input list.

Another pathway that belongs to both results but is only included because of the FDR-correction is “IL2”. IL-2 plays a role in cutaneous T-cell lymphoma [32]. This demonstrates that many of the enriched pathways are only accepted because the FDR-correction is still significant.

For comparing the results of PharMeBINet with other enrichment tools the same gene list is submitted to g:Profiler [15], ToppGene [36], DecoPath [16], Panther [37], and Enrichr [38]. If possible the threshold is set to 0.01 or manually filtered to have similar start conditions. Table 4 [Tab. 4] shows the number of enriched pathways returned for which enrichment method, database, and correction method.

Panther returns the lowest number of pathways for the input list. It is tested with Fisher’s exact test and binomial test with the FDR-correction. Both combinations return only three different pathways which are included in both PharMeBINet enrichments, as well (see Table 4 [Tab. 4]).

The next tool is g:Profiler which uses three different databases KEGG, Reactome, and WikiPathways with modified Fisher’s exact test and a custom correction method g:SCS threshold [15]. Each database returns its own pathways. The enrichment pathways from Reactome and WikiPathways are in the enrichment results of PharMeBINet. In contrast, four of the KEGG enrichment pathways are not in the PharMeBINet results (see Table 4 [Tab. 4]).

Enrichr uses Fisher’s exact test and the FDR-correction for multiple databases. Only the latest databases are utilized for the comparison: WikiPathways 2021 and KEGG 2021. The number of enriched pathways is higher than the number of g:Profiler. This may be because of the different database versions and/or because of the different correction methods. Not all WikiPathways 2021 results are included in the PharMeBINet results which may also be caused by the different database versions. The overlap between KEGG 2021 and PharMeBINet is low as seen in Table 4 [Tab. 4].

DecoPath uses Fisher’s exact test and FDR-correction in combination with the databases KEGG, Reactome, WikiPathways, PathBank, and their own DecoPath database. The returned enrichment pathways are grouped by their respective database as in the other tools. This causes some pathways to appear multiple times from different databases and the reason why in Table 4 [Tab. 4] the number of common pathways is split into common pathways from PharMeBINet and common pathways from other databases such as DecoPath. More than half of the pathways of DecoPath appear in the PharMeBINet results and the modified Fisher’s exact test has a greater overlap than the Fisher’s exact test.

The last tool, ToppGene uses Fisher’s exact test, FDR-correction, and multiple databases KEGG, MSigDB C2 BIOCARTA, BioCyc, Reactome, GenMAPP, Pathway Interaction Database, PantherDB, Pathway Ontology, and SMPDB. Also, like in DecoPath, all database results are in a single table increasing the number and causing the appearance of duplicated pathways. The overlap from ToppGene to PharMeBINet Fisher’s exact test is high and nearly all pathways are included in PharMeBINet results. The overlap to the modified Fisher’s exact test is even higher as the test returns more pathways.

In total, the results of PharMeBINet enrichment are on par with other tools. It does not yet include all pathways so a logical next step is to add more pathway databases such as KEGG.

Panther and ToppGene only provide table representations for the enrichment results. These tables are similar to the table of PharMeBINet.

Enrichr provides a tabular representation as well but also a log bar plot of the p-values/combined score rank/rank value for the top ten results. Additionally, it has a heat map for the top 30 enriched pathways with an additional bar plot in the graphic with log p-values/combined score rank/rank value. It is similar to the heat map and bar plot illustration of PharMeBINet. In contrast, in PharMeBINet all pathways are shown. However, Enrichr has additional plots on Appyter.

DecoPath also provides a tabular representation of the resulting pathways. Additionally, a pie chart shows the number of pathways which are concordant (multiple databases say it is significantly enriched), no mapping (pathways that are not mapped), and discordant (multiple databases say opposites). An additional table shows the pie chart’s information. This representation is unique. The last plot of DecoPath is a pathway hierarchy that highlights the significant pathways. This is a unique plot as well.

g:Profiler provides a Manhattan plot and a combination of table, log bar plot, and heat map. This tool has some plots in common with PharMeBINet but includes some unique plots, too.

However, none of them shows a graphical representation of how many of the genes are used or not. They only represent this in text form. Additionally, the last bar plot in PharMeBINet is unique, showing the number of enriched pathways the genes appear in.

Here, we present a new enrichment analysis module for the PharMeBINet website. It allows setting different parameters manually. The pathway enrichment for the use-case retrieved good results for Fisher’s exact test and modified Fisher’s exact test. Multiple figures visualize the results of the enrichment analysis. Not only the ranks of the pathways are shown but also the overlap of genes. Additionally, a diagram shows how many genes have no connection to a pathway. Furthermore, a bar plot demonstrates the number of enriched pathways each gene appears in. The enrichment already returned good results and got similar results in comparison to other online tools. However, the enrichment analysis could be improved if more pathway databases are considered. Further, it would be a good idea to add other enrichment algorithms than Fisher’s exact test as other methods are often used in research analyses as well.

• Cassandra Königs: 0000-0001-8991-9272

The authors declare that they have no competing interests.