Comparative study on gene set and pathway topology-based enrichment methods

Analysis of gene expression experiments typically yields long lists of differentially expressed genes (DEGs) [1]. To ease interpretation of such results it has been proposed to analyse expression data in the context of biological pathways or other biologically meaningful gene sets instead of individual genes [2]. There are many tests, which aim to detect pathways significantly enriched between two experimental conditions [3–6]. These tests were implemented into a plethora of methods that differ in many aspects ranging from null hypothesis formulation, through statistical framework up to pathway data encoding, database support and software availability.

We were interested to evaluate different methods for analysing gene expression data in the context of signalling pathways, with a special focus on comparing methods with different pathway representation strategies. There are different approaches how to handle a pathway in enrichment analysis. The traditional approach considers a pathway as a simple gene list omitting any knowledge of the gene and protein interactions. Methods belonging to this category are commonly called gene set (GS) analysis methods. Another way of integrating a pathway into enrichment analysis takes into account its graphical structure and/or topological measures. One of the first pathway topology-based (PT-based) methods was impact analysis proposed by Draghici et al. [7]. Since then this has approach become very popular, resulting in a number of PT-based algorithms being published [8–11].

Enrichment methods can be divided into two groups based on the null hypothesis: (i) competitive methods which are naturally linked with gene sampling for p-value calculation and (ii) the group of self-contained methods which is associated with subject sampling [12]. The basic difference between these two groups is that competitive methods compare genes in a pathway to its complement usually represented by the rest of the genes measured in the experiment, while self-contained methods consider only genes within a pathway and test their association with the phenotype. Both approaches have their limitations. On the one hand, competitive methods coupled with a gene sampling model for p-value calculation assume independence of genes, which is known not to be true. On the other hand, self-contained methods were criticized for being too powerful and thus yielding too many significant gene sets. Also the number of replicates in the experiments is often too low for subject sampling [12].

Another classification of enrichment methods separates them into over-representation analysis (ORA) and functional class scoring (FCS) approaches [13]. ORA is commonly referring to 2×2 table methods such as Fisher’s exact test, hypergeometric test and chi-squared test. It represents exclusively the competitive approach. The ORA methods require a strict cut-off in the list of DEGs and therefore the results are strongly dependent on the chosen threshold. FCS comprises methods, which first score genes and then transform gene level scores into pathway level scores. This group includes both competitive and self-contained methods, depending on the pathway-level transformation and significance assessment of pathway level score. Most of the PT-based methods can be classified into either ORA or FCS, however in some cases it might be difficult to draw a strict line.

From the user’s point of view there are several features of the enrichment methods to be considered besides classification of the methods according to methodological and statistical framework. For instance, the implementation and availability of the software can determine how popular a method becomes. Furthermore, several pathway databases such as Kyoto Encyclopedia of Genes and Genomes (KEGG) [14], Reactome [15], Pathway Interaction Database [16] or BioCarta [17] make data available for pathway analysis. While it is rather straightforward to supply simple gene set of a pathway for GS methods, providing a pathway with topology information and integrate it into a PT-based method can be more difficult. Therefore, the utility of PT-based methods also depends on their support of different pathway databases.

Several studies comparing different GS analysis methods were already published [18–20], as well as review papers offering surveys of PT-based methods [13, 21]. However, to our best knowledge there is no paper comparing the GS against the PT-based approach systematically.

To perform a fair comparison of the different methods we faced two major challenges: First, it was necessary to provide the same pathway data inputs for all methods, which was in particular challenging for PT-based methods. The integration of topological pathway data is described in detail in the third section of the Methods. Second, to simulate expression data for evaluation of PT-based methods was a comprehensive problem, which resulted in a complex simulations scheme. Hence, we described Simulations as a separate chapter.

In this chapter we present in detail 3 GS and 4 PT-based enrichment methods and their use within this work. Further, we describe how pathway data were processed and integrated into the methods. The last section of the Methods introduces the analysis of benchmark data.

Pathway data

Real pathway data available from public sources were used in this work. Within R various approaches to integrate pathways into enrichment analyses are available [25]. The three R-packages of PT-based methods comprise default pathway data of KEGG and other databases (see Table 1) in a preprocessed and suitable format. However, to ensure comparability of the evaluated algorithms we supplied all methods with the same KEGG pathway data.

Parsing the KEGG database

A BioPAX level 3 export of KEGG pathway database was downloaded on March 2013 and parsed into R using the rBiopaxParser R-package [26] (Fig. 1a). Non-metabolic pathways were selected and transformed into interaction graphs in which the nodes represent genes and the directed edges represent activation or inhibition processes between the genes.

After mapping HUGO Gene Nomenclature Committee (HGNC) gene symbol IDs on the graph nodes several editing steps had to be taken (Fig. 1b). First, in case that several nodes were annotated with the same gene symbol, these nodes were merged into a node, which shared all incoming and outgoing edges of the original nodes. Next, protein complexes or gene families often occupied a single node resulting in multiple symbol IDs embedded in this node. Such a node was split into multiple nodes and each one was assigned with a single gene symbol. Finally, nodes representing small molecules and other non-gene components were removed in a fashion that the parents and children of such a node stayed connected.

Transformation into suitable inputs

Depending on the method, suitable pathway input has to be provided implying graphs transformations [27]. For GS methods pathway graphs were converted into simple lists of genes. However, for the PT-based methods specific pathway topology inputs were required (Fig. 1c).

To create the SPIA pathway input, a graph of each pathway was transformed into a list of 2 adjacency matrices for activation and inhibition processes. Accordingly, a vector of weights β was set to β = {1, −1} to reflect activation and inhibition. For CePa a pathway catalogue was constructed comprising a list of pathways with the interaction IDs, and a table with columns for the interaction IDs and for the corresponding input and output interaction components. The pathway catalogue also included a mapping table. However, in our pathway data a single gene ID was always mapped onto one node. Input for PathNet consisted of an adjacency matrix of a pooled pathway, which was created by merging all database pathways, an interaction table with pathway IDs and a mapping table. The mapping table was required as PathNet algorithm transforms gene IDs into numeric objects and artificial integer IDs had to be generated.

All pathway data inputs were regenerated for the second simulation study (see Simulations chapter). In this study pathway graph nodes were relabelled with new synthetic IDs to construct non-overlapping pathways with unique components, while the topology of the pathway remained intact.

We generated synthetic expression data for two simulation studies each comprising 5 simulation types. Simulation types differed in the topology designs for pathway deregulation and in a choice of the parameter whose ranges were investigated (see Fig. 2).

The two distinct simulation studies differed in their gene ID annotation. While in the first study real gene HGNC symbols were used, in the second study these symbols were replaced by unique synthetic IDs for all nodes of all pathways in order to prevent overlapping of the pathways.

In both studies expression data were drawn from a multivariate normal distribution, always consisting of 10 control and 10 treatment samples. As treatment samples we refer to the group in which we introduced changes in the p-dimensional mean vector. By introducing changes in the mean vector we controlled expression levels of the genes reflecting fold-changes between the two groups. The genes that had an increased expression above 0 are called affected. In the (p × p)-dimensional covariance matrix employed for the multivariate normal distribution, we introduced a correlation of 0.8 between genes in the same pathway and 0.05 otherwise (see Additional file 1B). Variance in the matrix was set to 2. The same covariance matrix was used for the generation of both treatment and control samples.

Topology designs for pathway deregulation

In order to deregulate a pathway within a stimulation we needed to affect some of its genes. We examined 3 approaches how to reflect pathway topology to allocate affected gene in a deregulated pathway: community, betweenness and neighbourhood approach (Fig. 3).

Communities in a given pathway were detected by an algorithm for greedy optimization of modularity implemented in the igraph R-package (fastgreedy.community function) [31]. It aims to find modules with dense connection between the module nodes and spare connections between nodes of different modules. For each pathway we searched for a community which represented 45-55 % detection call (DC). However, for some pathways several communities had to be joined or a too big community had to be cut to get the appropriate DC (Fig. 3b).

To select affected genes in the betweenness deregulation design we considered the top highest scored betweenness nodes (Fig. 3c). Betweenness of a node is defined as the number of all shortest paths in a directed graph going through a given node. The required number of the top scored nodes for the given DC was assigned as affected.

In the neighbourhood deregulation approach one node of a pathway is chosen and then all nodes within certain distance create the neighbourhood (Fig. 3d). In the search for the neighbourhoods representing different ranges of DC we calculated the neighbourhoods of several orders for each pathway nodes and chose the one best fitting the required DC level (+/− 5 %).

In order to place long lists of differential genes into a context of biological processes and pathways the enrichment analysis strategy is widely used. In this work we compared standard gene set (GS) enrichment tests with methods, which attempt to incorporate pathway topology information into analysis.

Comparing GS methods and pathways topology-based (PT-based) enrichment approaches is challenging in two regards. First, the pathway input data have to be corresponding for all methods to ensure a fair comparison. Therefore, we parsed a BioPAX export of the KEGG database and transformed the pathways into gene sets as well as various interaction graph representations for the different methods. Second, it is difficult to simulate expression data for PT-based methods. The signal in a pathway is mediated via multiple biological mechanisms. Capturing a pathway deregulation in a graph and reflecting this alteration in the synthetic expression data is an intriguing task. To our knowledge there is no coherent work published on this topic. We made one of the first attempts to generate expression data mirroring the topology of the deregulated pathways. However, there are several limitations in our simulation scheme. For instance, multivariate approaches usually require more complex covariance matrices for a fair and unbiased evaluation, therefore, we have not included multivariate self-contained methods.

We decided to use 3 approaches to allocate affected genes to a pathway: gene community, betweenness and neighbourhood. However, in a certain setting some methods might be favoured due to their inherent algorithms, e.g. both CePa methods by betweenness and PathNet by neighbourhood designs. Further, it could be questioned how well these network measures reflect the real perturbation signals in a pathway.

In the first simulation study we could not identify any PT-based method as outstanding neither for any parameter configuration, nor for the overall accuracy measure. The pathway input for this study was represented by original KEGG pathways, which exhibit considerable gene overlap among each other. When a pathway was called deregulated it received a certain number of genes assigned as affected according to the detection call level. However, a single gene was on average present in 2.4 pathways. This resulted in a number of pathways, which were not considered as deregulated but also contained affected genes. These pathways then showed up as false positives and consequently led to a very low specificity and accuracy of all methods in study 1. However, one could assume that these ‘accidental’ affected genes were allocated randomly in the non-deregulated pathways, whereas in the deregulated pathway they were placed in a topological context. Therefore, it could be expected that the PT-based methods would overcome the problem of these false positives, however, this seems not to be the case.

It was already demonstrated that ranking methods are able to detect modest signals of deregulation better than the ORA approach [20]. We further showed that GS ranking methods were more sensitive than any of the PT-based methods when gene expression changes were weak (mean = +/− 1). The CePa GSA method exhibited a clearly different behaviour than the others, which was expected from the single self-contained method included in this work. Prominently different performance of this self-contained method could be further enhance by the fact that CePa computes for each pathway five p-values corresponding to different centrality measures and the authors of the method recommend to use the lowest value. Therefore, its overall sensitivity in the first study was higher than all the other competitive methods [12]. On the other hand, its specificity was utterly lost when a lot of pathways were deregulated (study 1, N = 70). Similarly, in the diseased benchmark datasets CePa GSA was the best in identifying target pathways but it also identified more than half of all KEGG pathways as significant. Depending on a disease, multiple pathways could be altered. However, it could be questioned whether such a number of pathways is realistic or it reflects lack of specificity of this method.

In the second simulation study with non-overlapping pathways the overall specificity markedly increased, mainly for the PT-based methods. This was also reflected in the best accuracy for these four PT-based methods.

Interestingly, the overall specificity of KS was the lowest among all methods and it was also always decreasing when more genes were affected (see Additional file 4, DC = 70 %, N = 70, size = big). That could stem from its null hypothesis, as it actually tests whether a gene set or its complement is significant [32]. In a setting when it was the complement it resulted in the false positive pathways in our simulations.

Because of the limitations of the simulations, which can never completely capture complex biology, we compared methods also on the real expression data. It is not possible to judge overall accuracy on these 36 benchmark datasets but PathNet and WRS showed reasonable ability to identify the target pathways. This suggests that simple GS approach might be sufficient for detecting enriched pathways.

Enrichment analysis is a useful and popular approach to ease interpretation of high throughput data. We comparatively investigated 3 gene set (GS) and 4 pathway topology-based (PT-based) methods in extensive simulation studies and on real benchmark data.

Our results are consistent with some already published findings: As expected, the self-contained method behaved differently than the competitive ones. The ranking methods were more sensitive than over-representation analysis methods when expression changes were moderate. Further, we demonstrated that on original KEGG pathways none of the PT-based methods was clearly outperforming the GS approach, neither in simulations nor in benchmark testing. This changed in the simulation study with non-overlapping pathways, where the PT-based approach exceeded simple GS tests. This suggests that further conceptual work has to be done to deal with the problem of pathway overlaps.

From the gene set methods Wilcoxon rank sum test performed comparably as topological methods in most settings, and it is much easier to use and re-implement with different input data than PT-based approaches.

Including topology to test for enrichment of a pathway is an interesting approach. It introduces additional information into the analysis and allows attractive visualizations of the pathway graphs. However, it needs to be further explored whether the additional biological insights justify the increased complexity of the PT-based enrichment analyses.

The data sets supporting the results of this article are included within the article and its additional files.