Pathway-GPS and SIGORA: identifying relevant pathways based on the over-representation of their gene-pair signatures

Algorithm

As illustrated in Fig. 2, our approach to the problem consists of two phases: In an offline phase, we compile a set of weighted markers/Signatures for each pathway in a repository, which we call Pathway Gene-Pair Signatures (Pathway-GPS). Subsequently, in an online phase, the method identifies the statistical over-representation of such Signatures in a user-specified gene list using an adapted version of the hypergeometric test.

Given a pathway repository (e.g., KEGG), for each gene-pair in a pathway, SIGORA investigates the co-appearance of the two genes in other pathways of the repository. A gene-pair that uniquely occurs in a single pathway is considered a Signature of that pathway and is assigned a weight. The weight of a Signature (from [0, 1]) quantifies the average commitment of the components of the GPS towards the common pathway, i.e., the weight scores the reliability of the Signature as evidence for the associated pathway. For hierarchically organized repositories (like REACTOME (Matthews et al., 2009)), this process is repeated iteratively after the removal of pathways on the top level of the repository, i.e., in each iteration, new weighted Signatures are identified for the pathways on the lower, more specific levels of the hierarchy. Once this offline stage is completed, the resulting sets of weighted gene-pairs that represent each pathway are non-overlapping and can be re-used for pathway analysis of any user-specified gene lists.

When presented with a gene list of interest (e.g., genes that are differentially expressed), SIGORA determines which of the pairs from its (pre-compiled) Signature repository can be reconstructed from the genes in the list. A Signature is considered “present” only if both of its constituent genes are found in the user-specified query list. This inherently leads to the selection of the more relevant roles of a gene in the experimental context, as SIGORA relies on the status of the other genes in the pathway for the reconstruction of the Signatures. For each pathway, the weights of present Signatures are summed up and hypergeometric probabilities are used to assess the statistical significance of the observed Signature sets.

Pathway Gene-Pair Signatures (Pathway-GPS)

A pathway database/repository contains (at least) two types of entities: pathways and genes. This can be represented by a bipartite graph (or bipartite network) B = (V g, V p, E) with two distinct sets of nodes (Vg: gene nodes and Vp: pathway nodes) where the edges in E connect the genes to the pathways and signify the annotation of a gene in a particular pathway (Fig. 3A). In this graph, the degrees (number of incident edges) of the pathway nodes correspond to pathway sizes (i.e., the number of genes annotated in the pathway) and the degree of the gene nodes corresponds to the number of different pathways a gene is annotated in. In particular, genes with degree one are exclusively annotated within a single pathway (‘Pathway Unique Genes (PUGs)’).

Using the igraph (Csardi & Nepusz, 2006) package of the R-Statistical Framework (R Development Core Team, 2008), B can be manipulated and transformed. A weighted one-mode projection (Newman, 2001) of B into the gene dimension, yields a new graph Proj = (V g, E′, W) with only one type of nodes (only gene nodes Vg), where two genes are connected by an edge in E′ if (and only if) they are co-annotated in one or more pathways, and each edge is associated with a weight (in W) which signifies the number of pathways that are shared between the two incident genes (i.e., the two genes connected by the edge). Two genes connected by an edge of weight one in Proj are, as a combination, unique to a single pathway. We call such pairs Pathway Gene-Pair Signatures (Pathway-GPS) of that pathway. The process of identifying PUGs and Pathway-GPS for all pathways in a given repository is termed the ‘Signature Transformation’ of the repository.

Hypothetically, one could go beyond gene-pairs and also consider n-tuples (with n > 2) of genes that co-occur in a single pathway P as signatures of P. The possible benefits of such extensions, however, do not currently seem to justify the associated complications to the computational and methodological framework. This point is further discussed in the Text S1.

Assignment of weights to GPS

As, by definition, each Pathway-GPS is uniquely associated with a single pathway, identifying such Signatures in a gene list of interest (e.g., observing that both constituent genes of a Pathway-GPS are in the list of differentially expressed genes) can serve as an indicator of the activation/perturbation of the associated pathway. Each pathway can (and usually does) have multiple possible Pathway-GPS, and (as discussed below) the method does not rely on the observation of an individual GPS but rather on the statistical over-representation of multiple GPS in comparison to the expected proportion.

Yet before the over-representation of GPS can be used for the identification of relevant pathways, we need to emphasize that different GPS vary in their reliability as indicators of a particular pathway. This is due to the fact that, while the two genes comprising a Pathway-GPS can only be co-annotated in a single pathway, each of the two genes – considered individually – can be a member of several distinct pathways.

Consider GPS of the form (g1, g2) for a pathway P, where g1 is annotated in i pathways and g2 is in j pathways. Intuitively, a GPS that consists of two PUGs (i.e., a case where i = j = 1) is a more appealing signature than a GPS that consists of two ‘multifunctional’ genes (say, i = 4 and j = 3), where the simultaneous observation of the two genes might be due to other factors (e.g., simultaneous activation of two different pathways).

To address this issue, a weight is assigned to each GPS to quantify its reliability as an indicator of its associated pathway. Let (g1, g2) be a GPS associated with a pathway P and let i and j be the number of individual pathways annotations of g1 and g2, respectively. The weight of the GPS is ${\displaystyle w_{i,j}=\frac{1}{2}\left(\frac{1}{i}+\frac{1}{j}\right)=\frac{i+j}{2*i*j}.}$ A motivation for this weighting scheme and a discussion of alternative weighting strategies can be found in the Text S1. In the implementation, the user also has access to several alternative weighting schemes and user defined schemes are also supported. Figure 3 summarizes the main ideas behind the compilation of weighted Signatures for a given pathway repository.

Identifying statistically over-represented Pathway-GPS

Analogous to traditional (individual gene) ORA (IG-ORA) methods, the distribution function of hypergeometric probabilities is used to calculate p-values indicating the statistical enrichment of Pathway-GPS in a user-specified gene list, which is given by: ${\displaystyle p(k,n,m,N)=p(x\ge k)=1-\sum _{x=0}^{k-1}\frac{\left(\begin{array}{c}\hfill m\hfill \\ \hfill x\hfill \end{array}\right)\left(\begin{array}{c}\hfill N-m\hfill \\ \hfill n-x\hfill \end{array}\right)}{\left(\begin{array}{c}\hfill N\hfill \\ \hfill n\hfill \end{array}\right)}.}$ In individual gene ORA, k denotes the number of query genes in the tested pathway, n the number of the genes in the pathway, N the number of all (assayed or annotated) genes, and m the length of the query list (i.e., the number of the genes with interesting status). In contrast to these traditional approaches, however, the parameters of the hypergeometric function in SIGORA are calculated as sums of GPS weights rather than frequency statistics of individual gene annotations. The Signature ORA parameters for p-value calculations are summarized in Table 2.

Note that – strictly speaking – hypergeometric probabilities are only defined on natural numbers while the sums of GPS weights are positive floating point values, which is why the table refers to rounded values of the sums (to the closest integer). Theoretically, this rounding could lead to a ‘blurring of the weights for pathways with few GPS’, however, in practice this issue hardly materializes, as such pathways are very rare (e.g., 219 out of 226 KEGG human pathways in our repository are associated with 10 or more possible GPS). In tests with simulated and real biological data sets, different choices of the rounding strategy (floor, ceiling or nearest integer) did not substantially affect the significance or the rank order of the identified pathways. As an aside, in the widely popular statistical framework R, the phyper function (which computes the distribution function of the hypergeometric distribution) does accept non-integer (floating point) parameters and handles such input by applying the following rounding strategy: the number of successes is rounded down, all remaining parameters are rounded to the closest integer value.

In our implementation, the user has the option to restrict the GPS-sets for the universe (N), and the success states (n) by providing a list of assayed genes (background). Furthermore, all PUGs are by default considered to represent a GPS of weight 1 (as a combination of the PUG with itself), but the user has the option to restrict the analysis to pairs of genuinely distinct genes.

Multiple testing correction

Undertaking pathway analysis generally involves a large numbers of significance tests. As testing a multitude of hypotheses will inevitably lead to some ‘significant’ results, adjustment of p-values for multiple testing is a crucial feature of any pathway analysis tool. Bonferroni’s method is used by default in SIGORA for multiple testing correction (MTC). It is, however, relatively easy to change the MTC procedure, if a user prefers to explore other adjustment methods (see the implementation section).

Selection of cut-off threshold for statistical significance of pathways

The choice of a reasonable cut-off threshold for statistical significance of the hypergeometric test results is an open methodological question in IG-ORA. In practice, values smaller than 0.1 or 0.05 after correction for multiple testing are commonly considered significant. Shifting the perspective from individual genes to the weighted Gene-Pair Signatures brings an additional challenge: as the size of the universe for the hypergeometric test dramatically increases (we move from a few thousand genes to up to a few hundred thousand weighted gene-pairs), the calculated p-values become by several orders of magnitude smaller than those observed in a typical IG-ORA analysis.

Based on our experience with simulated and biological datasets, we recommend a significance threshold of 0.001 after MTC (by Bonferroni). In the implementation, the default output of the analysis is the ranked list of pathways that achieve a corrected p-value up to this value. The user can also export the entire results table (including the p-values, corrected p-values and the parameters of the hypergeometric test) and the evidence (lists of present PUGs, list of the genes involved in present GPS or list of all present GPS along with their weights).

Dealing with redundancies of semantic origin

Thus far, we have described the motivation for Signature Transformation from a biological perspective. However, in some cases, there are additional semantic reasons for sharing of components among pathways. In particular, some repositories (e.g., REACTOME (Matthews et al., 2009), INOH (Yamamoto et al., 2011) and Gene Ontology (Ashburner et al., 2000)) are organized in a hierarchical structure, where all genes associated with one pathway (child) are also included in a more general pathway (parent). This poses a general challenge for pathway analysis tools that ideally should identify the most relevant level of the hierarchy (Alexa, Rahnenführer & Lengauer, 2006; Grossmann et al., 2007; Jupiter, Sahutoglu & VanBuren, 2009).

The hierarchical nature of such repositories poses a special challenge for our method. As any gene-pairs from a child category also co-occur in the parent pathway, they would thus be excluded from being identified as a possible Signature. This would have the undesirable effect that all child pathways on the lower levels of the hierarchy would be left without any Signatures at all in the offline stage and hence be undetectable by SIGORA in the online stage.

To address this issue, we deploy the following iterative top-down strategy in the offline stage (Signature Transformation):

1. Set level = 1. Compile the repository Signatures as described for the non-hierarchical case. Assign the compiled signatures to level 1.

2. Remove all pathways from the top level of the hierarchy, increase the level and recompile Signatures for the remaining pathways. Assign the GPS to the current level. Iterate this step until no further hierarchical levels can be removed.

Figure S1 illustrates this procedure on a simplified hierarchical repository. In the online stage (the identification of significantly over-represented Signatures in a user-specified gene list) the user can specify how many levels of the hierarchy should be considered in the analysis. Any GPS (and any pathways) that are deeper down the hierarchy (i.e., are at a higher threshold level) are left out of the universe (and the analysis). If the user (for instance) asks for analysis up to the second level, then only GPS from levels one and two are considered.

The effect of this simple modification (i.e., the iterative strategy for Signature Transformation) is similar to a combination of the Elim (Alexa, Rahnenführer & Lengauer, 2006), and TreeHugger (Jupiter, Sahutoglu & VanBuren, 2009) algorithms in dealing with GO’s structure. Elim essentially excludes genes in more specific categories from consideration in more general categories, whereas TreeHugger weakens the contribution of such genes to the higher levels.

Complexity and computational cost of the method

The complexity of the transformation step (in an implementation based on the outlined bipartite network interpretation of gene-pathway membership) is dominated by the complexity of the bipartite projection, which is given as: ${\displaystyle O(\Vert V\Vert *d^2+\Vert E\Vert )}$ With ‖V‖: number of nodes in the network (number of genes + number of pathways), d the average degree of nodes, ‖E‖ the number of edges in the bipartite network.

For most inputs, the online phase (i.e., analysis of user specific input lists using the precompiled Signatures) completes within 15 s on a standard laptop (1.8 GHZ, 2 GB of RAM).

Implementation

SIGORA is implemented as an R package. The package and a detailed manual are available for download from http://sigora.googlecode.com/svn/. The following highlights the most important steps in a typical work-follow in R (requires R version ≥2.10):

## install and load the downloaded package 
>install.packages(‘sigora_0.9.8.tar.gz’, type=‘source’, repos = NULL) 
>library(‘sigora’) 
## (please note that all of the following commands require that the package is already loaded) 
## import the query list from a file (assuming the list is given as Ensembl gene IDs): 
>myquerylist <-ens_converter(scan(‘myfile.txt’,what=‘character’)) 
## Alternatively, if the file consists of Entrez gene IDs 
>myquerylist <-entrez_converter(scan(‘myfile.txt’)) 
## perform signature over-representation analysis, using KEGG GPS 
>sigs(myquerylist,‘k’,markers=1,level=2) 
## multiple testing correction is done by Bonferroni and FDRs are also provided. In order to add Hommel’s method: 
>cbind(summary_results,p.adjust(summary_results[,5],‘hommel’)) 
## export the results into a file 
>export_results(filename=‘my_results.csv’, genes=T) 
## help (on Windows systems, help is shown in the web-browser) 
>help(sigs) 
## the following demo is also available 
>demo(sigora)

Evaluation methods

We evaluated the performance of SIGORA by comparison to several other analysis tools on simulated and published biological datasets. Three of the methods compared are based on individual gene over-representation (DAVID (Huang, Sherman & Lempicki, 2009; Jiao et al., 2012), gProfileR (Reimand et al., 2007), InnateDB (Lynn et al., 2008)) and the remaining three (GSEA_Preranked, GSEA (Subramanian et al., 2005), GSEA_AF (Ma, Sartor & Jagadish, 2011)) are GSEA based. Acknowledging the inherent challenges associated with comparison of p-values across different statistical frameworks, for each tool we follow the recommendations of the authors of that tool regarding the choice of the most appropriate significance threshold and multiple testing correction (MTC) method (Text S1).

Among the methods listed above, we use DAVID, gProfileR, GSEA_PRERANKED for the simulation study, as these tools (like SIGORA) can be run on any pre-selected gene list and in ‘batch mode’. GSEA and Appearance frequency modulated GSEA (AF) are limited to particular experimental designs and have specific input data format requirements, and are used here only in the analysis of three biological dataset for which data in the required input format was available. InnateDB is used as a reference point in evaluation of the biological datasets, because SIGORA’s GPS are based on pathway annotation data as present in InnateDB. The rationale for selecting GSEA is its popularity; while AF was chosen because it attempts to address similar issues as SIGORA.

A short summary of the relevant characteristics of each of these methods is given below.

Simulation experiment

Biological datasets

We further compare each of the different pathway analysis tools by examining their results when applied to three different gene expression datasets. These datasets incorporate the results of microarray studies investigating the host response to a parasite infection (Experimental Cerebral Malaria), a viral infection (Dengue Fever), and a bacterial infection (Tuberculosis) (Lovegrove et al., 2007; Thuong et al., 2008; Hoang et al., 2010). We map the lists of differentially expressed genes in each experiment (as provided by the authors of the respective studies) to unique Ensembl/Entrez IDs and use the resulting sets as input lists for four over-representation based methods (InnateDB, DAVID, gProfiler and SIGORA), ensuring that all methods are run on identical input lists. Additionally, two GSEA-based methods (GSEA and AF) are also applied to the corresponding expression datasets obtained from the gene expression omnibus (GEO): GSE25001, GSE11199, GSE11199.

Results on simulated gene lists

To evaluate the performance of SIGORA, we compared it to three other popular pathway analysis methods (DAVID, gProfiler and GSEA_Preranked) applied to simulated data, where we know a priori which are the significant pathways. The simulated input data was created by randomly choosing five KEGG pathways and selecting a fraction (alpha, 50% or 15%) of the genes in each of these pathways as being “differentially expressed”. Each of the methods was then applied to the selected gene list to determine the statistically significant pathways (using the respective recommended significance threshold and MTC approach, see Text S1).

If the sharing of genes between different pathways was not a factor, we would expect that each method should identify only the five preselected pathways as significant. As can be seen in Table 3, in our experiments (using 1,000 simulated datasets with alpha = 50% and 1,000 datasets with alpha = 15%), this was not the case: gProfileR, for example, identified more than 60 pathways on average as being significant at alpha = 50%, despite only 5 pathways being simulated as significant in the input data. SIGORA performed best by this measure and identified on average 8 pathways as significant across the different datasets.

The Recall and Precision metrics in the third and fourth columns of Table 3 capture the relationship between originally preselected ‘target’ pathways and the identified pathways. More specifically: Recall describes the fraction of the target pathways that were identified as significant and Precision signifies the fraction of statistically significant pathways that were among the originally selected pathways. Neither of these two metrics by itself is decisive, and there is a certain trade-off between the two metrics. This trade-off is captured by the F1-score, the harmonic mean of Precision and Recall. As can be seen in Table 3, SIGORA had the best F1 scores when compared to the other methods.

Finally, note that although this analysis is based on recommended significance thresholds in each method, Precision, Recall and F1 score are all dependent on the – ultimately arbitrary – choices of significance thresholds. The last column in Table 3 describes the results according to a less threshold-dependent measure. Ideally, the preselected pathways would occupy the first five positions in the list of the identified pathways, resulting in an average rank of 3. As can be seen in the last column of Table 4, at both choices for alpha (15% or 50% of the genes in each pathway selected) the average rank of the originally preselected pathways (‘target pathways’) in SIGORA’s results (3.6 and 3.8 respectively) is very close to this ideal value, whereas the other methods tend to identify several additional pathways as more significant than the target pathways, resulting in higher average ranks of target pathways in these methods.

Results on published datasets

In addition to the simulated data, we also compared SIGORA to five different methods (InnateDB, DAVID, gProfileR, GSEA and AF) applied to real biological data, in this case three different gene expression datasets. Full details of the input genes, p-values and highlighted pathways obtained by each method can be found in the Tables S1–S3.

Tuberculosis

SIGORA was compared to five different pathway analysis methods (InnateDB, DAVID, gProfileR, GSEA and AF) applied to a gene expression dataset (GSE11199) which measured the host transcriptional response in human macrophages infected with Mycobacterium tuberculosis (Thuong et al., 2008). 1,250 transcripts (corresponding to 1,100 distinct Ensembl genes) were identified by Thuong et al. as being induced in response to this infection. Figure 4 shows the pathways that were identified as statistically significant by each of the six methods. The first thing that one notes is that the GSEA-based methods tended to predict large numbers of pathways as statistically significant (the AF method predicted nearly a third of the KEGG database as significant in this example). This is somewhat by design, as GSEA methods attempt to identify subtle but coordinated changes in gene expression. This may be very helpful in investigating cases where there are only subtle differences between conditions but in a dataset like this one, it leaves the biologist bewildered as to which pathways should be followed-up on experimentally. In the other extreme is DAVID, which predicted only 4 pathways as statistically significant. SIGORA, on the other hand, identified 12 pathways as statistically significant; 10 of which were also identified by at least two other methods. Comparing the pathways identified by SIGORA as significant to the significant pathways identified by the other methods, one can see (Fig. 5) that many of the pathways identified by other methods as significant but not by SIGORA share many genes with the SIGORA pathways. Notably, after removing the multifunctional genes that are involved in the pathways identified by SIGORA from the input list, the individual gene over-representation based methods (DAVID, InnateDB and gProfileR) did not return any significant pathways at all. This reinforces our observation from the simulated data that SIGORA will identify truly significant pathways but avoid identifying pathways that are significant because they share genes with other more relevant pathways.

Interestingly, SIGORA may also be able to identify some important pathways that are not significant using other methods. One pathway was identified as significant in this dataset only by SIGORA; FcγR-mediated phagocytosis. Fcγ receptors regulate immune activation and susceptibility during Mycobacterium tuberculosis infection (Maglione et al., 2008; Maertzdorf et al., 2011) and it has been implied that “entry through Fcγ receptors may specify a distinct intracellular trafficking pathway for virulent M. tuberculosis” (Ernst, 1998). Finally, DAVID’s top functional cluster for this dataset (Enrichment score 0.99) contained 16 pathways, 12 of which are different cancer subtypes (Table S1).

Table S1 lists all Pathways identified as statistically significant by each of the considered methods and their respective ranks (by p-value) in analysis of this dataset.

Experimental cerebral malaria

Example 2 is a mouse cerebral malaria (ECM) dataset, comparing the whole-brain transcriptional responses of genetically susceptible (C57BL/6) and resistant (BALB/c) inbred mouse strains 6 days after infection with Plasmodium berghei ANKA (NCBI GEO: GSE7814) (Lovegrove et al., 2007). We fist performed a differential expression analysis using Geo2R to obtain a list of up-regulated genes at FDR <0.01 (637 Ensembl genes, Table S2). We use this list as input for SIGORA, InnateDB, gProfileR and DAVID, and apply GSEA and AF to the corresponding expression matrix.

Similar to the previous example, the number of KEGG pathways that different methods identified as statistically significantly enriched in this dataset varied widely: AF identified 59 pathways (out of 185 in its repository), while DAVID highlighted just two pathways. SIGORA identified 14 pathways as significant, 13 of which were also reported by at least two other methods (Fig. 6, Table S2). The remaining pathway is PPAR signaling pathway (discussed further below).

Aside from this big-picture view, the strong changes in the rank orders of the following individual pathways are also worth mentioning:

Dengue fever

As a third evaluation set, we re-examined a list of 483 up-regulated genes in the whole blood transcriptome of patients infected with dengue virus (NCBI GEO: GSE25001) (Hoang et al., 2010). More specifically, we compared expression profiles of hospitalized patients with uncomplicated Dengue during acute phase (≤72 h of illness history) to follow up samples of such subjects two weeks after discharge (n = 72).

Here, for the most part, the SIGORA results contain well-defined immunity related pathways. As before, some additional, potentially spurious pathways that are identified by other methods are not significant in the SIGORA analysis. Examples include “Staphylococcus aureus infection” (the second ranked pathway in InnateDB results) and the “Prion Disease” pathway (identified by both InnateDB and AF). These two pathways share components with the complement pathway (4 and 6 up-regulated genes respectively, Fig. 7), which has been shown to have a role in neutralising Dengue (Shresta, 2012).

Again, this example contains a possibly relevant pathway that was significant in SIGORA’s results, but overlooked by the other methods: “Lysosome”, at position eight of SIGORA’s results (Table 4, Table S3). Recent experimental evidence suggests that manipulation of the host’s autophagolysosomes by the dengue virus is an important part of the virus’s life cycle (Khakpoor et al., 2009; Heaton & Randall, 2010).

DAVID does not return any functional clusters for this dataset.

Discussion and related work

The existence of shared components between pathways poses a challenge for pathway analysis methods: which of the statistically significant pathways associated with such components are the most biologically significant? In 2005, Khatri and Drăghici surveyed the state of the art analysis methods and tools of the time and outlined several limitations as the challenges for the next generation of analysis tools. One such challenge “is related to genes that are involved in several biological processes. For such genes, all current tools weight all the biological processes equally. At the moment, it is not possible to single out the more relevant one by using the context of other genes [of interest] in the current experiment” (Khatri & Drăghici, 2005).

The years since have seen many interesting developments in the functional analysis of high throughput biological data: Several methods (Alexa, Rahnenführer & Lengauer, 2006; Grossmann et al., 2007; Jupiter, Sahutoglu & VanBuren, 2009) have been described to highlight processes at the appropriate level of specificity and to reduce the redundancy of the results in the context of Gene Ontology (GO) analysis, where the overlap between categories are often due to the hierarchical organization of the ontology. These methods, however, fail to account for overlap between pathways that don’t involve the full inclusion of all members of a pathway in another pathway. To deal with such cases, Antonov et al. (2008) proposed the creation of new functional categories as complex Boolean combinations of available GO terms. Unfortunately, such combinations are often not easy to interpret and a comprehensive search over all possible combinations is computationally infeasible.

Outside of GO, a few methods have been proposed that indirectly tackle the issue by either discriminative treatment of individual genes or alternative representation of the pathway repository in specific scenarios. An example of the later approach (alternative representation of the pathway repository) is Bayesian Pathway Analysis, BPA (Isci et al., 2011). BPA transforms each pathway in a pathway repository into a separate Bayesian Network (BN) and scores the fit of each model with the experimental (expression) data. BPA is expected to leverage the expression status of other genes in the experimental context, as BNs, in contrast to simple lists of genes, are deemed capable of accommodating local interactions between genes. Although highly sophisticated, BPA is computationally intensive and by design limited to the interpretation of expression datasets.

An early example of the former approach (non-egalitarian treatment of individual genes) is impact-analysis (Draghici et al., 2007). Impact analysis integrates the magnitude of each gene’s expression change along with the type (e.g., receptor, transcription factor) and position of each gene within the given pathways and their interactions into the statistical framework, however, the authors do not explicitly address the problems related to component sharing among pathways.

More recently, it has been proposed to add an appearance frequency based parameter to the statistical framework of GSEA. This additional parameter is intended to weaken the contribution of genes with multiple pathway memberships to the statistical significance of all of their associated pathways (Ma, Sartor & Jagadish, 2011). While this is a significant step in the right direction, the addition of such a parameter does not exploit the status of other genes in the experiment for the selection of the most relevant function of a gene in the experimental context. As exemplified by the gene, BRCA1 in Khatri & Drăghici (2005), even key players of one process (maintaining genomic stability) can have several less prominent roles in other unrelated pathways (response to nutrient and brain development). Nor can appearance frequency distinguish between the causes of appearance of a gene in several pathways, which, aside from functional pleiotropy of genes, can be partly due to the hierarchical organization of some pathway repositories like REACTOME and GO.

Here we introduced the concept of Pathway-GPS, as genes or gene-pairs that (as a combination) are specific to a single pathway, and we described Signature Over-representation Analysis (SIGORA) as a novel (and comparably efficient) approach to pathway analysis. SIGORA uses Pathway-GPS to bridge the gap between the context sensitive and collaborative nature of biological processes on one side and the universal and discrete statistical framework of over-representation analysis on the other side. Although each Signature’s weight is fixed in advance, the net contribution of an individual gene G to the measured success (parameter k in Table 1) of each of its associated pathways is not fixed and depends explicitly on the status of its partners (genes that together with G form a GPS).

In contrast to the GSEA-based solutions, SIGORA inherits the versatility of the ORA statistical framework and is applicable to lists of genes of interest obtained in any type of high throughput experimental set-up (e.g., copy number variations from cellular profiling, lists of epigenetically silenced genes from promoter methylation analysis, differential gene expression data from NGS and microarrays or SNPs from GWAS experiments) without the need for adaptation of the computational method. This is especially notable in situations where ranking of the entire dataset by a single biological parameter (as required by GSEA) is not feasible (see Huang, Sherman & Lempicki, 2009, for a few examples).