 Research article
 Open access
 Published:
Robust extraction of functional signals from gene set analysis using a generalized threshold free scoring function
BMC Bioinformatics volume 10, Article number: 307 (2009)
Abstract
Background
A central task in contemporary biosciences is the identification of biological processes showing response in genomewide differential gene expression experiments. Two types of analysis are common. Either, one generates an ordered list based on the differential expression values of the probed genes and examines the tail areas of the list for overrepresentation of various functional classes. Alternatively, one monitors the average differential expression level of genes belonging to a given functional class. So far these two types of method have not been combined.
Results
We introduce a scoring function, Gene Set Zscore (GSZ), for the analysis of functional class overrepresentation that combines two previous analysis methods. GSZ encompasses popular functions such as correlation, hypergeometric test, MaxMean and Random Sets as limiting cases. GSZ is stable against changes in class size as well as across different positions of the analysed gene list in tests with randomized data. GSZ shows the best overall performance in a detailed comparison to popular functions using artificial data. Likewise, GSZ stands out in a crossvalidation of methods using split real data. A comparison of empirical pvalues further shows a strong difference in favour of GSZ, which clearly reports better pvalues for top classes than the other methods. Furthermore, GSZ detects relevant biological themes that are missed by the other methods. These observations also hold when comparing GSZ with popular program packages.
Conclusion
GSZ and improved versions of earlier methods are a useful contribution to the analysis of differential gene expression. The methods and supplementary material are available from the website http://ekhidna.biocenter.helsinki.fi/users/petri/public/GSZ/GSZscore.html.
Background
The analysis of differential gene expression between two sample types, such as pathological and healthy tissues, is one of the cornerstones of the modern biomedical science. Here typically the up or downregulation of each gene in the pathological samples is measured. The obtained expression data can be considered as Ordered Gene List, OGL, by sorting it according to gene regulation. The upper end of the OGL represents the strongest upregulation and the lower end the strongest downregulation in the pathological sample. The middle area of the list represents genes with insignificant regulation. Similar gene lists can also be generated with various other data sources, like sequence similarity searches, high throughput screening of gene knockouts, or protein expression arrays.
Originally, the analysis of differential expression focused on the few genes with the most extreme up or down regulation between the samples. This is sensitive to potential measurement errors, and prone to produce false positive findings, such as genes reacting to any sample handling. Furthermore, it cannot detect biological processes showing consistent but weak regulation between the two datasets. A central improvement was to carry out the analysis using predefined gene classes or gene sets, such as functional classifications [1–3]. The predefined gene classes usually show more stable behaviour as a group than at the single gene level, and classes also simplify the biological analysis [4, 5]. Functional classes are typically used to study a subset of the genes, like the ones whose expression exceeds a predefined threshold. This subset is usually analyzed for overrepresented functional classes, using Fisher's exact test. Currently several web tools, such as DAVID [6], implement overrepresentation analysis.
Although the utilization of the gene annotation classes moves the focus away from the actual single gene level to the more robust biological process level, it is strongly dependent on the definition of the used threshold. This has motivated many authors to propose threshold free gene set analysis methods (also called as Gene Set Enrichment Analysis [5] and Threshold free ontological analysis [7]). These methods monitor differential gene expression at the class level. They can be roughly divided into two categories: Signal summary methods take the whole gene list into account and calculate a mean or median based score for all the member genes of the class [8, 9]. Ranked list based methods [5, 10, 11] analyze overrepresentation of the gene class by going through the whole OGL and test every possible position for the threshold. Usually the threshold with the strongest score is selected. The KolmogorovSmirnov score (KS) [5] and the hypergeometric pvalue [11, 12] represent two popular test scores for this analysis. The latter, hence called here iGA (iterative Group Analysis, as in [11]) represents a simple extension of the standard threshold based analysis. We refer to these functions, used to monitor differential gene expression at class level, as class level scoring functions.
Both types of methods have different drawbacks. Signal summary scoring functions assume that the whole gene set will show homogeneous behaviour. This is violated, for instance, when the gene set constitutes a pathway with genes repressing and inducing the amount of endresult, or when the gene set includes a significant amount of genes misclassified as members of the functional class. On the other hand, the ranked list based methods omit the score associated with each gene, and consider the step between consecutive genes always constant. Yet, in the real datasets one can see very drastic regulation at the tails of the gene list with clear differences between genes and small differences in the middle region with more or less random gene ordering. Indeed, one of the ranked list statistics has been criticized for being too sensitive to the signals in the central area of the gene list [13].
Here, we propose an improvement that avoids these problems by merging the two method groups. The main principle is to apply a statistic that is based on under or overrepresentation and to analyze the whole list for the strongest signal. The difference to standard ranked list methods is that the actual scores for the differential expression, assigned to each gene, are taken into account. A similar idea has already been presented [13], but the unique feature here is the derivation of the mean and the standard deviation of the score, allowing direct normalization as a Zscore (Gene Set Zscore, GSZscore). As a result, the regularized version of the GSZscore shows very stable behaviour, especially when compared to the similar earlier method [13]. GSZscore can be linked with two standard test scores, the hypergeometric test and correlation between the gene set labels and the expression scores, and with two recent scoring functions, Random Sets by Newton et al. [14] and maxmean used by Gene Set Analysis [15]. Thus GSZscore represents a novel unification between these different relevant scoring functions. It allows monitoring over and underrepresentation of the gene class within the subset of the gene list and also consistent shift of the gene class as a whole. The regularized version of the obtained GSZscore is stable across different threshold positions in the gene list and across the different GO class sizes. To our knowledge, this is the first publication to analyze these stabilities. Furthermore, we present extreme value distribution (EVD) as an asymptotic distribution for GSZscore pvalue calculation. The presented method and test datasets are available from our web site http://ekhidna.biocenter.helsinki.fi/users/petri/public/GSZ/GSZscore.html.
We comprehensively evaluate the performance of GSZscore and various other frequently used class level scoring functions. We focus our testing on class level scoring functions, by keeping everything else in the analysis pipeline identical. This ensures that the observed differences come from the evaluated functions. The same cannot be guaranteed when comparing software packages. Similar modular analysis of gene set analysis methods was recently recommended elsewhere [16]. We first test class scoring functions with artificial datasets. Differently from earlier similar work [7, 10], we tested scoring functions with classes showing upregulation and also with classes showing simultaneous up and downregulation. GSZscore shows the best overall performance in these tests. In addition, we show optimal regions from the artificial signals for different scoring functions. Here, GSZscore has a good performance with signal types that were considered important in earlier publications [5, 7].
Next, we do repetitive testing of class scoring functions on real data set. As the correct positive and negative classes are not known in real data, we used the evaluated functions to predict them. We used Cross validation style evaluation, where data was split to halves, with the first half used to predict the positive and negative classes and the other half used to test the competing methods. Test situation was replicated by (i) using each function itself to select the positive classes from the other half and (ii) by combining the results of all the methods from the other half and using this same set for all methods. Both tests were repeated with different number of positive classes. In addition, we further confirm our results by monitoring the correlation between the test and control results. GSZ constantly shows the best performance over different test situations, whereas the performance of other functions varies between the test situations.
We further evaluate the analyzed class scoring functions by generating large set of randomizations. Two very different datasets were used here. We monitored the magnitude of different empirical pvalues, in order to see the risk of potential false positive findings for each function. Results showed the best performance for GSZ with small difference on one dataset and with very clear difference on the other dataset. We point out that all the observed differences in numerical data evaluations resulted from the class scoring functions, establishing GSZscore as outstanding in this group. In addition, real data analysis shows weak performance by some of the popular scoring functions. This is, to our knowledge, one of the most detailed numerical evaluations of class scoring functions on real datasets.
We also carried out a detailed biological evaluation of scoring functions. This was done by ordering GO classes with empirical pvalues and monitoring top GO classes. Next, we defined biologically relevant classes for both datasets and analyzed how different methods were able to find them. GSZscore clearly shows outstanding performance in finding the important biological themes on both datasets. Again some of the competing scoring functions evidently show weak performance on one or both of the datasets. Results were further analyzed by monitoring GO classes with largest pairwise difference between the GSZscore and competing scoring function. Classes where GSZscore outperformed others showed stronger differences in pvalues and were clearly linked to the research setup.
We also performed a comparison with some of the most widely used software packages in this research topic. This analysis is more difficult to evaluate due to the other potential sources of differences between software packages. However, this ensures that we evaluate the state of the art implementations of class scoring functions. We concentrate here on the evaluation of the reported list of GO classes, and mostly omit the reported pvalues. GSZ again showed best performance on two tested datasets, when evaluated with biologically relevant classes.
Results
A typical pipeline for gene set analysis starts by looking at the differences in the gene's behaviour in the two datasets, like in treatment and control. This can be done, for example, by calculating a score based on the difference between the mean expression values, like ttest score, with standard ttest or its variations [17, 18]. We call these primary tasks differential expression tests and their scores as differential expression test scores for clarity. The obtained gene list is then sorted using the selected score. In the next step, gene classes, one at the time, are taken to the analysis and the gene list is analyzed using either rank based methods or signal summary methods. We call these tests class scoring functions and their results class expression scores. Our aim is to improve the existing class scoring functions. For more information on this topic see earlier publications [5, 10, 11, 13].
Novel Gene Set Zscore
We start with a ranked list based analysis, where a threshold is placed repetitively in between every consecutive pair of genes in the ordered list. The obtained subset is then used to calculate a selected statistical score. This test score can be a difference between two empirical cumulative distributions like in KolmogorovSmirnov statistics [5] or a log(pvalue) calculated from Hypergeometric distribution (iGA, [11]) among others.
We propose to improve these tests by taking the differential expression scores into account. Let M denote the total number of genes, and let X_{ i }, i = 1,..., M, denote the differential expression test score for the i th gene. Furthermore, let S_{ N }be a subset of N genes from the upper end of the OGL (i.e. N genes with the highest differential expression test score values X_{ i }). We use a simple function:
This denotes a difference between the sums of scores for members (positive gene group, pos) and nonmembers of a gene class (negative gene group, neg), among the N genes with the highest differential expression test scores. This will then be calculated separately with each threshold position N, and analogously for the lower part of the OGL. To simplify the notation, the subset S_{ N }will be kept implicit in our notation by dropping the subscript N in the sequel. A similar but more complex idea has been proposed in the article that motivated this research [13] (see also Methods in supplementary text S1 [additional file 1]). However, we decided to use the simpler equation as:

Zscore normalization is easier to perform for the plain difference.

The obtained Gene Set Zscore (GSZscore) includes other popular methods as limiting cases (discussed below).

The resulting score is only affected by differential expression test score values within the selected subset, whereas in the equation in [13], the whole distribution of class member gene scores affects the weighting of the class members. Our idea should therefore, in principle, be more efficient in observing simultaneous up and downregulation.

As a more theoretical benefit, the obtained GSZscore treats both class members and class nonmembers similarly. So the results stay identical if we flip the class member and class nonmember classifications.
The selection of the starting equation 1 is still somewhat an open issue and our main emphasis is rather on the Zscore normalization for gene set presented later. The similarities with other class scoring functions are discussed more in detail in the supplementary text S2 [see additional file 2], where we show that GSZscore links together Random Sets scoring function [14], maxmean for Gene Set Analysis [15], hypergeometric Zscore and correlation between the class labels and the gene expression scores.
The raw difference is an unusable statistic due to the biases caused by different sizes of classes, different sizes of subsets, and different variances of the scores in different subsets. However, these can be corrected by calculating estimates for the expected value (mean, E(Diff)) and variance (D^{2}(Diff)) for the eq. 1 under the null hypothesis, that the class members and nonmembers are distributed randomly across the list. This allows us to use a Zscore normalization with equation:
where E(Diff) = 2E(X)E(N)  ME(X) and
where N is the number of positive genes in the subset, M is the size of the subset, E(N) is the mean and D^{2}(N) is the variance of the hypergeometric distribution of the number of positive genes N for the analyzed subset. D^{2}(X) and E(X) are the variance and the mean of the differential expression test scores for the subset. The derivation of E(Diff) and D^{2}(Diff) is shown in Methods. As a further improvement to our GSZscore we consider a regularized version:
This is otherwise similar to eq. (2) but we add a prior variance k to the variance estimate inside the square root. This was done due to the too unstable behaviour with small subsets (see supplementary text S1, [additional file 1]). The regularization represented here should not be mixed to normalization that is performed later with the randomizations. Note that prior variance is regularly used with ttest in the gene expression data analysis [17, 18]. Supplementary text S1 [additional file 1] shows that the regularized version has a more stable behaviour. We have used two rules for defining the prior variance k: We take the median of the variance estimates, obtained with the analyzed class across the whole gene list, and multiply it with a weight, w2 (0 ≤ w 2 ≤ 1). We also take the median of the variance estimates, obtained with the class of size 10 across the whole gene list, and multiply it with a weight, w1 (0 ≤ w 1 ≤ 1). Next the weighted variance medians are summed to obtain a value for k. Class size of 10 is used here as a reference class size, penalizing especially classes smaller than 10. This class size selection is 'ad hoc', but it has performed well in our analysis. We tested 12 different combinations for these two weights, with artificial data, and selected the best performing parameter values for real dataset analysis. The calculated GSZscores form a profile of score values over the OGL obtained from each cutoff position, including also the values corresponding to the whole gene list as a subset. From the obtained scores we select the maximum value similarly to earlier works [5, 11, 13]. Other options are stated in discussion.
Stability of various class scoring functions with row and column randomizations
First we tested the stability of various class scoring functions with randomized real datasets. The tested scoring functions in the results are: GSZscore, standard and modified KolmogorovSmirnov test (KS, modKS [5, 13]), threshold free implementation of hypergeometric test (iGA, [11]) and ttest (with added prior variance) calculated between the expression values of class members and class nonmembers. Ttest reflects here the behaviour of signal summary methods [8, 9]. We use two types of randomizations: in row randomization the gene classifications associated with dataset rows are randomized, and in column classifications the sample labels associated with data columns are randomized. The aim was to see what type of biases the methods show as the GO class size varies, or as the threshold moves through the gene list. These results are shown in the supplementary text S1 [see additional file 1]. Here, we highlight a few details: The use of the prior variance k significantly stabilizes the GSZscore, with a group of tested parameter values for the k obtaining equally good results (see supplementary figure S1, [additional file 3]). Regularized GSZscore has a stable signal distribution over all the threshold positions with a randomized dataset, whereas two other evaluated methods, KS and modKS show significantly less stable behaviour over different threshold positions (see supplementary figure S2, [additional file 4]). Furthermore, a larger randomization shows very good stability for GSZscore in row randomizations as the size of the class is varied. A more detailed analysis is in supplementary text S1 [see additional file 1].
In addition, we show in supplementary text S1 [see additional file 1] the dramatic biases caused by gene level correlations, discussed also earlier [13, 15, 19]. These can be seen when analyzing the column randomizations. We show the bias only for the most stable methods from earlier row randomizations (see supplementary text S1, [see additional file 1]), but it can be expected to affect all the methods. Here, mean and STD estimates from the column randomizations were tested for normalization of the results. The stability of the estimates was evaluated by dividing the randomized dataset into testing and learning data, where learning data was used to obtain the estimates and the testing data was normalized using the estimates and visualized to see the stability of the normalization with randomized data. The obtained normalization looks adequate with 100 randomizations. Despite the adequate normalization, the later steps use larger numbers of randomizations, as they will enable better analysis of empirical pvalues for the observed results.
Evaluation of class scoring functions with artificial data
Next, we evaluate the class level scoring functions on very large set of various artificial datasets. We currently lack complete knowledge on what can be considered as a positive signal for a differentially expressed gene class. However, we propose a few signal types that can be intuitively justified. Case 1: The whole gene class shows a similar type of regulation, either up or down. Nevertheless, the complication in this simple case is that the signal for the gene class can be quite diffuse with larger STD, depending on the different measurement efficiencies for various RNAs. Case 2: only a subset of the gene class shows signal (up or downregulation) like in case 1. Potential explanations are: Only a subset of the gene class requires regulation or the gene class includes false positive genes that are not really members of the gene class. Case 3: Gene class shows both up and downregulation at the same time. Here gene class can comprise differently behaving subclasses or the members can be from a pathway with activities that either increase or decrease the activity of the pathway.
Here we have used an artificial data generation similar to previous work [7, 10]. Test covers an even larger number of various signal subtypes for positive classes than previous research. We vary the size of the gene class that represents the signal, the magnitude of positive signal, the standard deviation of the positive signal, and the percentage of gene class that shows positive signal, as in the previous work. The signal is also either set to be only in positive direction (only upregulation, case 1 and case 2), or the genes representing the signal are set to be randomly either positively or negatively regulated with equal probability (case 3). To our knowledge, this is the first work evaluating both of these cases to this extent. Details on the various parameters used in the artificial signal generation are shown in the Methods (supplementary text S1 [see additional file 1]).
First, the analysis was carried out with various gene class sizes, testing each class size separately. Thus, the compared results for positive and negative cases are obtained with the same gene class size. Next, the results for different class sizes and signal types are combined. This measures the performance of the method in the case where the various gene class sizes are separately normalized with mean and standard deviation estimates obtained from randomized data. This is currently the standard procedure of most programs. In another option, the positive and negative results for different class sizes are pooled and the separation is monitored within the pooled data. This measures the stability of the functions when the gene class size varies. This corresponds to the performance of the method without any class specific normalization with randomized data (proposed for example with iGA [11]).
In each testing situation with each parameter setting for positive datasets we generate 200 positive datasets and 600 negative datasets. Next, we calculate the test scores with each of these datasets. We used Area Under receiveroperating characteristics Curve (AUC) as a measure of separation. This was considered better than the analysis of power [7], as power analysis requires a selection of the pvalue cutoff, measuring separation only at that point, whereas AUC measures the separation across the whole list. We have first summarized the information from four testing situations by comparing every method pair and calculating two scores:
where i and j specify methods that are compared and l = 1, 2.. R refers to dataset, obtained with specific parameter values. Score A represents the mean of differences of AUC scores for two methods across all the datasets in the testing situation, whereas score B represents which of the methods, i or j, was more frequently better. We generated with both of these scores a difference matrix representing all the method pairs. Next the average for each method i was taken, by letting j go over all the methods. Now we have two measures for the tested method: One corresponds to average of difference to other methods and the other corresponds to frequency of better performance than with other methods. Results for these two scores are shown in table 1.
The results show the best overall performance for GSZscores. The only exception to the rule is the better performance of iGA in the testing with simultaneous up and downregulation with fixed class size. Even there GSZscores showed good performance in a larger subset of the datasets (shown by the score B). This discrepancy of two measures proposes that we frequently see slightly better performance from GSZscore and a small subset where iGA clearly surpasses GSZscores. Subsequent analysis links this to datasets showing signal around 1 STD (discussed later with fig. 1 and fig. 2). However, with real datasets this testing situation would require normalization of iGA results with several randomized runs, which is time consuming. Indeed, when the evaluation is done without class size normalization, iGA performance drops. Table 1 also suggests potentially best regularization parameter values for GSZscore. None of the parameter values is optimal in all cases. Yet w1= 0.1 or 0.2 and w2 = 0.3 or 0.5 show good performance across the tests. From the obtained results we selected w1 = 0.2 and w2 = 0.5 for later analysis, although we point out that other nearby parameter values yielded equally good results. ttest has been similarly evaluated with several parameter values (with w1 = 0 as the standard ttest), with the best performance always occuring with the largest value w1 = 3.
Table 1 summarizes a large amount of information, when in reality methods perform differently with different types of signal [7]. This was highlighted by visualizing the score for each combination of signal representing subset size and signal magnitude separately, like in previous work [7]. Each result is averaged over tested class sizes and standard deviations of the signal. We further modify the earlier visualization by marking the best ranking methods with colours (blue, cyan, green, yellow for the best, the 2nd, the 3rd and the 4th best). The following work concentrates on analysis done with each class size separately.
Fig. 1 and 2 reveal trends, pointing out that no test is optimal across all the signal types. GSZscore, with different parameter settings, is among the top methods across the most signal types. If we look at fig. 1, KS method shows the best performance with small signal magnitudes (signal between 0  1 STD. STD refers here to the standard deviation of the background distribution.) This was followed by iGA, showing best performance when signal is between 1.25  1.5 STD and signal representing portion between 0.6  0.9. GSZscore is optimal when signal is larger than 1.5 STD or the signal portion is smaller than 0.4. The only deviations from this trend are modKS, which is optimal when signal is between 1  1.75 STD and signal proportion is 0.1, and ttest, which is optimal when whole class shows signal that is between 0.75 STD and 2.25 STD. A simpler view is obtained by looking at results obtained with simultaneous up and downregulation (fig. 2). Here, good performance is only represented by KS test (when signal ranges 0.25  0.5 STD), iGA (when signal ranges from 0.75  1.5 STD) and GSZscore (when signal is larger than 1.5 STD). Similar results were obtained for cases where class sizes are allowed to vary. The only exception is that in these tests GSZscore and ttest clearly improve, and the performances of iGA and especially KS drop (data not shown).
Evaluation of the methods with real dataset
Previous analysis evaluated various methods using artificial datasets. The problem with such simulated data was that we do not have knowledge on how regularly one sees various signal types. Furthermore, the artificial data analysis assumed totally independent genes, which is a biologically unrealistic assumption. Therefore it is crucial to evaluate the performance also with real datasets. We use here ALL (Acute Lymphatic Leukaemia) dataset [20], due the reasons discussed in the Methods. These analysis steps are also represented as pseudo codes in supplementary text S1 [see additional file 1].
With a real dataset the true positive and true negative GO classes (also referred to as gold standard) are not known. Therefore we used the methods themselves to estimate the positive classes. This can omit some positive classes that none of the methods can discover and equally report some negative classes as positive. Nevertheless, we are not aware of any other realistic solutions for selecting positive GO classes, and minor errors should not corrupt our results. We also presume that the GO classes form a smooth transition with strongly positive GO classes at one end and strongly negative GO classes at the other end rather than a sharp binary classification to positive and negative classes. We therefore decided to use a varying rank threshold, selecting q best classes (q = 1, 2...200) from each method as positive (similar to [15]). This allows each method to report exactly the same number of positive GO classes and enables testing with a varying number of positive GO classes, concentrating to the best ranks in the selected positive GO class list. The competing idea, with selection of a fixed set of positive GO classes, would depend on the used pvalue threshold and on the used multiple testing correction. Furthermore, it would place methods depending on their stringency in unequal situation, as methods with a tendency to report many classes as positive would contribute a larger portion of positive GO classes. Also, the evaluation of sorted GO classes with decreasing significance resembles the popular analysis of GO classes as a sorted list.
We chose the AUC score to see how well the sorted list of GO classes separates different selected positive and negative classes. These results were also further confirmed with overall rank correlation between the ordered gold standard GO class list and ranked GO list from the evaluated method. Two methods have different aims: AUC with varying rank threshold focuses on most positive GO classes in control set, whereas the rank correlation monitors the overall correlation between the test and control set. Resulting evaluations are based on rank rather than the actual score values obtained for different methods. This was considered to be beneficial as different methods can be expected to follow different distributions, and including the actual score values would therefore put methods in unequal position. Only case where also the actual signal scores were used was the comparison of method with itself using pearson correlation (see table 2).
One of our evaluation setups requires the combination of the gold standards from all the evaluated methods. It would be natural to combine the rank results from various methods by averaging, but instead the best rank for each class was used here. This was done to minimize the effects of correlations between the methods, as with mean score two identically behaving methods would contribute twice in the positive class selection, whereas with the maximum score they would contribute only once.
The evaluation was done in two ways: i) Method itself is used to define k best GO classes from one half of the dataset. Next the method is run using the other half of the dataset to see how well it can predict its outcome from the correlating dataset. ii) All methods are allowed to select k best GO classes as positive GO classes using one half of the dataset. Next, each method is run on the other half of the data and it is tested how well they can predict the combined positive outcome of all the methods from the correlating dataset. Results from both cases were analyzed using AUC with varying rank threshold and also with standard correlations.
These evaluations have different evaluation principles. In case (i), we use each method itself to define gold standard from the separate data, selecting the method that is most robust to the variances between the biologically replicated datasets. Furthermore, as each method is tested against its own results, the occurrence or lack of correlation with other scoring functions does not affect the evaluation. However, this could favour methods that have biased selection, like a preference of large classes. Neither does it penalize methods that miss GO classes that other methods are able to find. In case (ii), the combined results of all the methods from the separate data represents the same gold standard for all methods. Here method has to be robust to variances between the datasets, and also be able to predict combined positive classes coming from all the methods.
These analysis steps were done with four splits of the dataset (8 replicates). Results, averaged over all the replications, are shown in fig. 3 and 4 and tables 2 and 3. Also each split was analyzed separately [see additional files 5 and 6]. Results vary considerably at small ranks between replicates, especially in case (ii), but they stabilize at larger ranks, especially when rank threshold is between 10 and 20. This can be explained by the increasing size of the positive class set, used in the AUC scoring, which is potentially less sensitive to variations and errors in its definition.
In case (i) GSZscore shows best performance, with an even performance to ttest at the very first ranks. The detailed analysis of differences in replicates (see supplementary figure S3, [additional file 5]) showed equal performance between ranks 140 and the performance of ttest drops to clear separation after rank 40. The observed difference in results suggests that GSZscore might see a larger set of classes with positive signal than ttest, and later analysis steps confirm this. Separation to other methods was obvious with iGA, modKS and KS showing consistently weaker AUC scores across all the ranks in all splits. These results are further confirmed with rank and Pearson correlations, which show the clear best performance for GSZscore (table 3).
In case (ii) GSZscore shows the best performance, and iGA is the close second best method. Their performances overlap in replicates across the top 10 ranks and after rank 20 GSZscore consistently tops iGA in all replicates (see supplementary figure S4, [additional file 6]). Ttest shows also somewhat equal performance at top 5 ranks but then it drops. ModKS shows also slightly better performance than GSZscore in two of the replicates at very smallest ranks. These are only a few rank positions, and the performance of modKS drops at these replicates extremely fast. Indeed, modKS is mainly 3^{rd}or 4^{th}method, and it is the weakest method at the smallest rank in one of the splits. After rank 20 the order of methods is consistently GSZscore, iGA, modKS or ttest and KS. Rank correlation analysis confirms also these results (table 2) showing the best performance for GSZscore.
The overall performance of GSZscore was best in these evaluations, as in case (i) it slightly surpassed ttest but is clearly better than iGA and in case (ii) it slightly surpassed iGA but was obviously better than ttest. Furthermore, it consistently topped modKS in case (ii) and clearly showed better overall performance in case (i). It was also observed that KS test shows the undoubtedly weakest performance in both tests. For a novel method, like GSZscore, it is also interesting to see which methods are most similar to it. This was analyzed with the whole ALL dataset, using rank correlation. Results are shown in table 4. The correlations that GSZscore shows with iGA and and ttest further underlines the GSZscore's mathematical similarities with these methods.
Comparison of empirical pvalue signals
Earlier evaluations omitted the number of significant classes reported by each scoring function. This is, however, a relevant feature as smaller pvalues and larger number of significant classes point out that the method is better in extracting signal. This evaluation was first done with ALL dataset [20]. Empirical pvalues were obtained by using 1000 row and column randomizations and selecting always the larger, less significant, of two pvalues. Only the scoring function is changed in this comparison and everything else is kept identical. We exclude the usage of asymptotic distributions as they are not available for all functions. All the results are shown here as log10(pvalues). We compared top 100 GO classes from each function by plotting the empirical pvalues for these classes. Cases where log cannot be calculated (empirical pvalue = 0) were modified by adding a value 0.5 to the calculus (see Methods, supplementary text S1 [see additional file 1]). We focus on the magnitude of reported pvalues, as the selection of significant classes would depend strongly on the used multiple testing correction.
The actual pvalue signals are represented as two separate results in the fig. 5. The fig. 5A shows pvalues from pooled data. Here each GO class is first normalized using the mean and the STD from the randomizations and next all the randomizations are pooled and used together to obtain pvalues. Fig. 5B shows pvalues from class data. These empirical pvalues are obtained by using only the class specific randomizations with each class. Our analysis showed that these two log(pvalue) results correlate strongly, and we further show in supplementary text S1 [see additional file 1] that our normalizations with randomized data was rather stable. Note that the class data has got a smaller intensity range and also discrete steps in the reported logpvalues, which is why it does not separate functions as well as pooled data.
Fig. 5 shows that the GSZscore obtains significantly larger log(pvalues) in the analysis of pooled log(pvalues) (fig. 5A). The difference is surprisingly large with the best log(pvalue) for GSZscore being 6.7, whereas the second best method (iGA) has best log(pvalue) 4.2. Furthermore, the difference is one to two orders of magnitude across the whole top class list (GSZscore pvalue is 10  100 times smaller than pvalue of any compared method). GSZscore also reports 23 classes with pvalue = 0 with class data, whereas other functions report only 1  3 classes. These are represented by largest values in the fig. 5B. These observed differences are larger than any of the differences between the other methods in the comparison.
ALL dataset represents very strong biological signals. However, it is equally important to evaluate the scoring functions with expression data with weaker signal levels. Therefore, we replicated the pvalue comparison with p53 dataset [13] and represent the results in fig. 6. Although the differences to other methods are now smaller, GSZ is clearly the best method in pooled data over top 15 classes. In class data the differences are smaller. GSZ produces again largest number of classes with pvalue = 0 (6 classes), but also KS produces 5 classes. In addition, GSZ drops to the 2^{nd}best after iGA around rank 10. This data proposes that GSZ is the best function across top ranks but it drops down later. These weaker ranks probably represent more random background signal.
Biological comparison of class scoring functions
We also evaluate the biological relevance of the reported classes. This monitors the ability of scoring functions to report the relevant biological classes. We generated a sorted list of GO classes from each scoring function using the pvalues from pooled data and class data in combination (see Methods, Supplementary text S1 [additional file 1]). This analysis is replicated on both real datasets, ALL and p53. The evaluation requires that we define functional classes that can be considered biologically positive, by selecting classes that can be linked with research setup. This selection is justified from the background biology, and it should not favour any of the evaluated scoring functions.
P53 dataset compares cancer cell lines with a mutation in p53 transcription factor to the ones without mutations. Here the classes with link to programmed cell death (apoptosis) were considered as positive. Selected classes were further confirmed by looking their regulatory relationship with p53 from the GeneGo Inc. MetaCoreTM database ([21], http://www.genego.com). As p53 data represents a small number of positive classes, we represent its results here in table 5. A more detailed presentation is in the supplementary table S1 [additional file 7]. Although most methods select the same apoptosis class (release of cytochrome c from mitochondria) as best, GSZ tops others by reporting 4 apoptosis related classes. GSZ reports also strongest pvalues to each of these apoptosis related GO classes. KS and, unexpectedly, also Ttest show here very weak performance.
As ALL dataset compared B and T cell lymphomas with each other, we considered classes linked to immune response as positive. These include classes like immunoglobulin complex, T cell receptor complex etc. Two groups of classes were considered as border cases: the generic MHC protein complex and interleukin proteins. Due to the large number of positive classes, these results are summarized in fig. 7. Each positive class is weighted here as one and each border case classes were weighted as 0.5. More detailed presentation of the results is in supplementary table S2 [see additional file 8].
Fig. 7 shows that our ALL dataset represents a very strong biological signal with GSZ, ttest, and iGA exhibiting 11  12 positive classes among the top 15 classes. ModKS and KS clearly represent weaker performance with 7 positives and one border case (modKS), and 3 positives and 2 border cases (KS) among the top 15 classes. Outside the top ranks, the density of positive classes in ttest and iGA test drops, whereas GSZ clearly shows the best performance. Results from both ALL and p53 dataset propose that GSZ is outstanding among the evaluated scoring functions in finding relevant biological signals.
Analysis of pairwise differences in class pvalues
In order to evaluate further the differences between the scoring functions we focused on the pairwise differences in the log10(pvalues) for each GO class between GSZscore and each of the other methods. This allows the analysis of the actual GO classes representing the strongest separation between the compared methods. Here, again, we monitor whether these classes are relevant to research setup and also the difference in the log(pvalues).
The outline of this analysis is represented in the supplementary text S1 [see additional file 1], and detailed results of the comparison in supplementary tables S3 and S4 [see additional file 9 and 10]. The top classes, where GSZ outweighs other methods are mostly all positive classes. Also the differences in pvalues, favouring GSZ, are very large. The opposite classes are more or less random selections from the GO classes.
Also pvalue differences are much smaller. There are also some interesting details in the comparisons:
GSZ vs. Ttest
When comparing GSZ with Ttest we observe large number of relevant GO classes with simultaneous up and downregulation (large STD signal, small mean signal) in ALL dataset. These were naturally missed by Ttest, whereas GSZ was able to report them. Furthermore, we observed surprisingly weak performance of Ttest on p53 dataset.
GSZ vs. iGA
The comparison with iGA showed many classes from ALL dataset that had a strong regulation. GSZ was able to report very strong pvalues for them, whereas iGA did not notice them since it discards the expression values. This underlines that inclusion of differential gene expression scores benefits the test statistic. GSZ clearly reported from p53 dataset better pvalues for 3 apoptosis classes.
GSZ vs. modKS
On ALL dataset modKS shows very strange behaviour, and it is clearly outperformed by GSZ. Our assumption is that weak performance of modKS is related to its sensitivity to outlier expression values in GO classes. This outlier sensitivity is further demonstrated by supplementary figures S1 and S2 [see additional file 3 and 4].
GSZ vs. KS
Standard KS showed weakest performance in these comparisons reporting mostly quite random GO classes.
These results bring forth GSZ as an outstanding method, when evaluating the biological relevance of the reported classes. There was only one exception on this rule in the p53 dataset, where iGA and KS reported stronger signal on ARF protein signalling classes. These apoptosis classes were not reported as linked to p53 according to GeneGO. Also their average regulation was very small.
Comparison of GSZ with program packages
We also compared our analysis pipeline with other actual software packages. This comparison is less clear to interpret, due to a large number of variables between different analysis pipelines. In addition, the programs report pvalues postprocessed in various ways. Therefore, the evaluation is focused on the order of the reported GO classes omitting the actual pvalues. These obtained GO class lists are evaluated using the selected biologically positive classes again. We selected three software packages for comparison. These were Gene Set Enrichment Analysis package (GSEA, R code version [13]), Signal Pathway (SP [9]) and Gene Set Analysis (GSA [15]).
All compared methods were run with the same set of GO classes, with no maximum size limit and the minimum size limit set to 3. Although these settings are suboptimal for the biological analysis, they allow a thorough evaluation of the packages against the variations in the GO class size. All the methods were also tested with 1000 randomizations. Both these settings correspond to the parameters used with the GSZ analysis in function comparison.
Two methods (GSEA and GSA) are onesided tests, generating two separate outputs. However, we need a single GO class list for comparison with SP and GSZscore. This was accomplished by combining the two GSEA output lists and sorting them using the absolute value of Normalized Enrichment Score (NES), used also to sort the result classes in the GSEA. With GSA we did a similar procedure, ordering in the first round with the normalized score (calculating Zscores with permutation results, returned by GSA) and in the next round using the pvalues. However, we also discuss the rankings in the separate lists. The following two chapters represent the results from two datasets.
ALL dataset
First we compare methods with ALL dataset. We used same biologically positive classes as in the earlier evaluation with the ALL dataset. The obtained results can be seen in the supplementary table S5 [see additional file 11] and the summary in the fig. 8. The first observation is that three out of four methods clearly report positive classes. These are GSZ, SP, and GSA. These results are in agreement with the strong signal that was reported from this dataset in the earlier comparison of scoring functions. GSEA, nevertheless, shows a strong deviation from this rule. None of the reported classes was able to obtain stronger signal than its default cutoff (False Discovery Rate ≤ 0.25). Furthermore, the reported classes are quite arbitrary with only 4 positive classes among the top70 classes. This is a very strong disagreement between GSEA and all the other evaluated methods.
There are also differences between the top three methods. GSZ, GSA, and SP all have only a few negative classes among top 20 classes. However, after the rank 20 their density drops in GSA and SP results. Indeed, there is only 7 negative classes among top 40 classes in GSZscore results and one border case, whereas SP results have 17 negative classes and one border case, and GSA has 18 negative. This makes their error rate among these top classes over two times bigger and highlights GSZ again as the best method. A detailed analysis shows that SP misses totally positive classes with large variance signal and small mean signal (thymic Tcell related classes and immunological synapse). Likewise, GSA represents a quite small signal for some classes emphasized by SP and GSZ, like 'MHC class II protein complex' and 'antigen prosessing and presentation' classes.
P53 dataset
Next, we compared the methods with p53 dataset. Here, again, we select the programmed cell death (apoptosis) related GO classes as positive. Obtained results are shown in table 6. More detailed results are shown in supplementary table S6 [see additional file 12]. Again, GSZ shows good performance with 4 apoptosis classes among its top15 results, showing strong ranking positions, such as 1 and 6.
SP reported one apoptosis related class at the rank 8, with empirical pvalue = 0. This was the same class rated as first class by modKS, iGA and GSZ. It had actually the first rank in the column randomization of SP, proposing that it is here more useful than the combined ranking from two randomizations. Thus, SP performs quite nicely, but reports only one cell death related class among the top50 classes, whereas GSZ was able to report 4 of them among top15 classes. Furthermore, the standard ranking of SP was not optimal here.
GSA generates separate lists for up and down regulated classes. Also, the outputs of GSA represent several GO classes with same pvalues. Therefore, we generated new ordering for GO classes, separately for upregulated, downregulated and combined set of GO classes. One apoptosis related class is seen in the downregulated set of GO classes with empirical pvalue = 0. Its original rank was 5^{th}and after our reordering with Zscores its rank is improved to 2^{nd}position. However, when the two GO class lists are pooled, the rank of the apoptosis class is only 6. Therefore, its performance seems to be weaker than GSZ.
GSEA generates two lists, instead of one. We again combined them into a single list. The apoptosis related classes were in the list of downregulated classes at the ranks 8, 14 and 24, and in the results of the pooled list at the ranks 9, 21 and 37. However, if we were to use the default GSEA FDR cutoff (0.25), none of these classes would have been reported. Also, all these apoptosis related classes were low in the result ranking. Therefore, the performance of the GSEA seems quite weak again.
Altogether, both datasets point out to GSZ as best method for detecting biologically relevant classes. SP and GSA show quite equal second best performance. Furthermore, both datasets pinpoint GSEA as the apparently the weakest method.
Biological analysis of results
GSZ reported 4 apoptosis related classes from p53 data: Release of cytochrome c from mitochondria, apoptotic mitochondrial changes, induction of apoptosis by intracellular signal and caspase activation via cytochrome c. What these classes represent is some of the key steps for programmed cell death: intracellular apoptotic signalling to mitochondria, removal of cytochrome c and its later complex formation with caspase. In addition, other methods (iGA and KS) highlight the ARF protein signal transduction. This group plays a relevant role in apoptotic signalling, although it shows here only very weak regulation. This suggests that a combination of two or three scoring functions could actually highlight all the possibly relevant biological signals in the dataset.
ALL dataset showed classes representing simultaneous up and downregulation, such as immunological synapse, positive thymic t cell selection and lymphocyte activation. Key trends here seem to be: Cell surface receptor signalling that activates immune response and developmental steps of T cells in thymus. Both these processes are central to immune processes. GSZ is able to recognize in particular the immunological synapse, a communication method between B cells and T cells. Note that part of the immunological synapse is expressed in T cell and part in B cell and it is complete when two cells interact with each other. These relevant immunological processes could not be observed from this dataset by methods that monitor only the mean expression of the GO class. Here ttest and SP missed them. Only GSZ and iGA performed satisfactorily on these classes. Also GSA from program comparison was able to detect a fraction of these.
ALL dataset also represented GO classes with strong mean up or downregulation. These classes link to T cell receptor complex, antigen processing via MHC class II protein complex, and to MHC class II protein complex. What is paradoxical is that some methods did not give good scores to these classes. Especially iGA and KS performed weakly with these. This was probably due to the limitation of the analysis only to the order of the genes. Only GSZ and ttest in the comparison of scoring functions reported good scores for these classes. Also SP gave very strong signals to these classes. Surprisingly, GSA obtained less strong results on some of these classes (MHC class II protein complex and antigen processing related classes). Altogether, these results point that across different biologically relevant signals GSZ shows better performance compared to the other, competing algorithms and programs.
Discussion
Current work proposes a novel class scoring function for threshold free gene set analysis. The function is a Zscore that is shown to have similarities to hypergeometric function and to correlation (see Results and supplementary text S2 in [see additional file 2]). Gene Set Zscore is shown to be stable under row randomization to variations of subset and class sizes. Furthermore, it can be stabilized with column randomizations against the biases caused by genelevel correlations. GSZscore had the best overall performance in the artificial data analysis. The second best class scoring function was iGA. It is, however, significantly heavier to compute than the GSZscore. In addition, we have performed a detailed evaluation on a real dataset, where the number of positive classes and the way how they were selected was altered. Here, GSZscore was consistently the best performing class scoring function, whereas the performance of other functions varied. A further proof favouring GSZscore was obtained when analyzing the empirical pvalues from two datasets. Here GSZscore was the top method, especially with ALL dataset showing 10  100 fold better pvalues across the top hundred classes. Altogether, these numerical results indicate that the GSZscore has an outstanding performance, when compared to the other represented class scoring functions.
We also monitored the biological performance of scoring functions by looking how they reported relevant classes. This evaluates the scoring functions for explorative biological comparison of sample groups. GSZscore showed best performance, reporting biologically relevant classes with different signal types from ALL dataset. Other scoring functions showed weaknesses, discussed in the background. From p53 dataset only GSZscore reports four relevant apoptosis related GO classes.
Biological analysis was also replicated against three widely used gene set analysis packages: Gene Set Enrichment Analysis, Gene Set Analysis and Signal Pathway. GSZscore clearly showed better performance in biological explorative analysis also against these tools. We omitted any detailed evaluation of pvalue signal between the program packages, since there is large number of variables between packages making comparison difficult. However, we point that GSZscore, GSA and SP show strong pvalues, whereas GSEA shows very weak signals.
It should be noted that we are excluding from the scoring function evaluation all the other (critical) analysis steps, and treat the functions in a consistent manner. This potentially weakens the performance of modKS, for which the original article [13] reported the need of separate evaluation of the upregulated and downregulated tail of the expression dataset. Also, the original work used a correlation test for differential expression, and used a different normalization with randomized data (for KS and modKS). These were excluded, as our aim was to keep the evaluation of the methods as constant as possible. Varying also other analysis steps would make it difficult to interpret what actually caused the differences in the results. Furthermore, our biological analysis shows better performance for modKS in our analysis pipeline than what is observed in GSEA, indicating that our analysis pipeline performs better.
Most earlier publications focus only on one of the two potential null permutations. We propose a separate normalization with row and column randomizations of a dataset, and the selection of the less significant outcome (normalized score value or pvalue) as the result. This represents a pessimistic but also intuitive perspective, where the more likely null model is allowed to explain the observed results. As a drawback, the obtained pvalue will naturally be conservatively biased. Nevertheless, our permutation evaluations show that GSZ is stable under row randomizations (see supplementary text S1 [additional file 1]). Also Efron and Tibshiriani [15] have proposed an alternative way to combine two nulls.
There are some details that were observed to cause problems when comparing gene expression data with functional classifications. These problems probably occur in many other biological data comparisons as well. A dataset should have a single, explicit representative for a single gene/transcript. Several measurements strongly violate the independence assumption, utilized by all the popular methods that we are aware of. Therefore, one can obtain seemingly highly significant GO class results when several measurements are associated with a single false positive gene. Similarly, it is unclear how to treat measurements (probe sets) that are associated with many genes. Fortunately, the current remapping projects (see for example [22]) aim to solve these types of problems, and turn out to be crucial for any gene expression data analysis. Furthermore, we have observed that genes not linked to classification should be filtered from the analysis. These are often classified as not belonging to the class, but such a solution caused false positive signal amongst the very largest classes of GO structure. We want to point out that our ALL dataset has been preprocessed to take these details into account.
Randomizations play critical role in analysis (see supplementary text S1, [additional file 1]). Row randomizations revealed substantial bias in some scoring functions (KS, modKS, normal ttest) and slight bias in iGA. Even more biased results were observed when doing the column randomization. Especially the iGA showed a dramatic bias here (largest log_{10}(pvalue)>70). Note that iGA represents a variation of the standard threshold based hypergeometric analysis, and these results suggest that one might need similar column randomization also with those methods [19]. Also some gene set analysis methods currently omit column randomization totally [7, 11, 14], which might cause bias in their outcome.
The evaluation of gene set enrichment type analysis methods is a nontrivial task, and its problems are not usually discussed. We can design testing with artificial data where we know the positive and negative classes. However, the relevance of different signal types is not known. With real dataset we have an opposite problem: We have a true signal across different classes, but the actual classification to positive and negative classes can be considered to be a subjective and vague decision. Here, we have first used the same methods that we are evaluating for the selection of the positive classes, but point out that this is not an optimal solution. Furthermore, we propose that it is reasonable to perform the testing so that the threshold for the positive class selection is allowed to vary. This enables the testing with varying number of positive classes for monitoring consistently good performance. Also the definition of the positive classes is varied. One can start with an assumption that the results from the method itself, from a similar dataset, represent the optimal solution that the method is then allowed to predict (referred sometimes to as concordance analysis. See similar work in [15]). The use of two datasets would correspond to a reallife situation where one laboratory is trying to confirm the results from another laboratory. Previous evaluation omits the positive results from other methods, and therefore we decided to repeat the analysis using results pooled from several methods. These test situations would correspond to a situation where the second lab is confirming results from the first lab (with one method), by using many methods in combination. A further improvement to the presented analysis would be the testing with smaller subsets of the original data, to potentially see further variations in the method performances.
Real data was evaluated even further by generating a 1000 row and column randomizations, and looking at the actual pvalues. We focused on the top 100classes from ALL dataset where GSZscore clearly represented a larger signal with 250  7100 times better signal for strongest class, and consistently over ten times better pvalues for top ranks. Obtained results were further analyzed by looking at the biologically positive classes. Analysis showed that GSZscore results show largest number of biologically relevant classes among their top ranks. Scoring functions were further compared using the pairwise differences of the GOclass pvalues. These differences from ALL dataset highlight a) classes where Ttest fails due to the heterogeneous regulation, b) classes where iGA represents weak signal although the class represents a clear regulation, c) strongly regulated classes where modKS fails totally. P53 dataset gave also similar results. Overall, our work represents one of the most detailed evaluations of GO analysis methods with artificial and real data, and we hope to inspire the field to implement similar detailed and comprehensive evaluations in GO analysis method comparisons.
The combination of artificial and real data analysis revealed some exciting insights. KS test had the third best performance in the artificial datasets, but with real data it showed the undisputedly weakest performance. This seems to confirm the hypothesis that the signal area where KS shows its best performance in artificial data (see fig. 1 and 2) is unimportant in the real datasets [7, 13]. Also modKS function showed very illogical behaviour with ALL dataset, preferring classes with smaller average regulation. Altogether, these results suggest that the use of KS and modKS should probably be avoided in GO class analysis. We observed even stronger differences in program package comparison. There GSEA package was indisputably the weakest program.
The introduced GSZscore could be further improved. One of its weaknesses is the requirement of the prior variance k, which corrects for signal when the subset from gene list is small. However, our artificial data evaluation shows set of parameter values that represent equally optimal performance (w 1 = 0.1, 0.2 and w 2 = 0.3, 0.5). These are used to add prior variance to calculus. Similar method is regularly used in gene expression data analysis and in various areas of data mining (like Ridge Regression, biasvariance tradeoff or addition of pseudo counts to the calculus). Yet the performance could potentially be improved with a more exact measure using an estimated cumulative distribution (for more exact pvalue) or Bayes Factor. In addition, the performance might be easily improved by using a larger number of score values from the calculated GSZscore profile, obtained over the gene list, than just a single maximum value. A separate variation on this theme is the connection to maxmean statistic [15], which was observed to correspond to GSZscore obtained when the threshold, selecting the analyzed subset, is set to correspond zero in differential expression scores (see supplementary text S2, [see additional file 2]). The presented GSZscore should be applicable also to other relevant research topics, such as feature extraction. In addition, the mathematical derivation of the GSZscore should provide an easy platform for further modifications.
Conclusion
With increasing demand of gene list data analysis for gene expression, the introduced Gene Set Zscore function should represent a significant addition to analyst's tool palette. It represents a novel addition of signallevels to ranked list analysis and gives reliable results.
Methods
Estimates for the expected value and the variance of GSZscore
To calculate the estimates for the mean and variance, we have to note that the analysis situation can be modelled with a composition of two functions. The first one defines the probability for the observed number of class members in the analyzed subset. The second one defines the probability for obtaining the observed sum for the class members. The exact analysis with Bayes factors or using exact cumulative distribution to calculate pvalues, would require the definition of a complex distribution with potentially too heavy calculations. Therefore, we prefer to use a Zscore with relatively simple estimates for mean and variance. Mean is simpler to derive from these two, and it is also required in the variance estimate. Thus we represent it here first. We start with the difference
Here pos refers to positive genes (class members) and neg refers to negative (class nonmembers). We are dealing here with a sum of N values S = Σ_{ n }X_{ n }selected from the pool of X_{1}, X_{2},..., X_{ M }values, which in turn is a subset of the total dataset (with size of L genes with K positive genes). The N = N_{ pos }is the number of positive genes included to our data subset (with M = N_{ neg }+ N_{ pos }). When the subset is randomly sampled from the total data pool, we do not know the N beforehand, but it can be expected to follow hypergeometric distribution, representing sample (without replacement) of size M from the pool of L genes including K genes classified as positive. Expected values for both terms in the eq. (8) can be defined using conditional probabilities
which is simply the probability weighted sum of expected values, conditional on the number of class members in the sum. Proof that the E(SN = i) = iE(X) is shown in supplementary text S1 [see additional file 1]. E(X) represents the expected value, or mean, of the data subset X_{1}, X_{2},..., X_{ M }. Expected value for eq. 8 can be represented as:
So the expected value is a simple function of expected value for hypergeometric distribution E(N), expected value for the test scores in the analyzed subset E(X) and the number of data points in the subset M. In practice we replace the E(X) with empirical mean calculated from the whole subset. This intuitive result is conditional on the selected data subset (M), on the size of the whole data pool (L) and on the number of positive genes in the whole pool (K). Estimate of the variance is somewhat harder to obtain. Here we start with the definition of variance for a single summation (S) in eq. (8).
By fixing N we can express the expectation E(S^{2}N) as a sum of variance and the squared expectation.
The first term in the latter equation represents the variance of the sum of N values selected from the pool of X_{1}, X_{2},.. X_{ M }values with variance D^{2}(X) (which is again replaced by empirical estimate obtained from the subset). Note that D^{2}(S) = 0, when N = 0, N = M or if X_{ i }= X_{ j }for all values of i and j (resulting to D^{2}(X) = 0). Derivation of D^{2}(SN), used in eq. 12 is represented in the supplementary text S1 [see additional file 1]. Substituting eq. 12 to the eq. (11) we get
This is simply an equation of the variance and the mean of hypergeometric distribution (E(N), D^{2}(N)), the mean and the variance for the subset of the data (E(X), D^{2}(X)) and the size of the data subset (M). The result stays same, whether we monitor the variance of sum for the negative or the positive data points. Equation is similar to the well known equation for the variance of the sum of N copies of identically independently distributed variables D^{2}(S_{ iid }) = D^{2}(X)E(N) + E(X)^{2}D^{2}(N). Here the difference is caused by various draws of X_{ i }not being independent.
The equation (13) represented the variance of the first summation in the equation (12). For the whole variance of the eq. (8), we have to multiply eq. (12) with squared constant 2^{2} = 4. Furthermore, we have to consider the variance of the latter term ME(X). A closer look reveals that this is simply a multiplication of two constants having no variance at all, no matter what outcome we observe in our test dataset.
Therefore it does not affect our variance estimate. So the final score for a data subset analysis becomes:
In contrast to KS and modKS, when OGL is divided to two subsets, our score gives a different score for upper and lower subsets. Therefore, with eq. 14 it is required to consider separately the lower and upper end of the list. We took the largest absolute score from these two lists as the final outcome. Another modification to our score function was made due to the observed instabilities with small subset and class sizes. This is simply a prior variance k added to the variance:
Note that similar procedure is regularly used with ttest in the expression data analysis. Selection of the prior variance will be discussed in the results section.
Summary of the remaining Methods
The BioConductor package [23], probe remapping [22], Robust Multiarray Average (RMA, [24]) and Intensity Based Modified Ttest (IBMT, [18]) were used in the gene expression data analysis. Gene IDConverter [25] and Protein Information Resource (PIR, [26]) were used in the generation of GO dataset. The rest of the Methods are shown in the supplementary text S1 [see additional file 1].
Abbreviations
 ALL:

Acute Lymphatic Leukaemia
 AUC:

Area Under receiveroperating characteristics Curve
 ECD:

Empirical Cumulative Density function
 EVD:

Extreme Value Distribution function
 GO:

Gene Ontology
 GEVD:

Generalized Extreme Value Distribution function
 GSA:

Gene Set Analysis
 GSEA:

Gene Set Enrichment Analysis
 IBMT:

Intensity Based Modified Ttest
 iGA:

iterative Group Analysis
 KS:

KolmogorovSmirnov test
 LR:

Likelihood Ratio
 ML:

Maximum Likelihood
 modKS:

modified KolmogorovSmirnov
 OGL:

Ordered Gene List
 RMA:

Robust Multiarray Average
 SP:

Signal Pathway.
References
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, IsselTarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nature genetics 2000, 25: 25–29. 10.1038/75556
Mewes HW, Dietmann S, Frishman D, Gregory R, Mannhaupt G, Mayer KF, Munsterkotter M, Ruepp A, Spannagl M, Stumpflen V, Rattei T: MIPS: analysis and annotation of genome information in 2007. Nucleic acids research 2008, (36 Database):D196–201.
Kanehisa M, Goto S: KEGG: kyoto encyclopedia of genes and genomes. Nucleic acids research 2000, 28: 27–30. 10.1093/nar/28.1.27
Hosack DA, Jr GD, Sherman BT, Lane HC, Lempicki RA: Identifying biological themes within lists of genes with EASE. Genome biology 2003, 4(10):R70. 10.1186/gb2003410r70
Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, Puigserver P, Carlsson E, Ridderstrale M, Laurila E, Houstis N, Daly MJ, Patterson N, Mesirov JP, Golub TR, Tamayo P, Spiegelman B, Lander ES, Hirschhorn JN, Altshuler D, Groop LC: PGC1alpharesponsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nature genetics 2003, 34(3):267–273. 10.1038/ng1180
Jr GD, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, Lempicki RA: DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome biology 2003, 4(5):P3. 10.1186/gb200345p3
Nilsson B, Hakansson P, Johansson M, Nelander S, Fioretos T: Thresholdfree highpower methods for the ontological analysis of genomewide geneexpression studies. Genome biology 2007, 8(5):R74. 10.1186/gb200785r74
Pavlidis P, Lewis DP, Noble WS: Exploring gene expression data with class scores. Pacific Symposium on Biocomputing 2002, 474–485.
Tian L, Greenberg SA, Kong SW, Altschuler J, Kohane IS, Park PJ: Discovering statistically significant pathways in expression profiling studies. Proceedings of the National Academy of Sciences of the United States of America 2005, 102(38):13544–13549. 10.1073/pnas.0506577102
BenShaul Y, Bergman H, Soreq H: Identifying subtle interrelated changes in functional gene categories using continuous measures of gene expression. Bioinformatics 2005, 21(7):1129–1137. 10.1093/bioinformatics/bti149
Breitling R, Amtmann A, Herzyk P: Iterative Group Analysis (iGA): a simple tool to enhance sensitivity and facilitate interpretation of microarray experiments. BMC bioinformatics 2004, 5: 34. 10.1186/14712105534
Berriz GF, King OD, Bryant B, Sander C, Roth FP: Characterizing gene sets with FuncAssociate. Bioinformatics 2003, 19(18):2502–2504. 10.1093/bioinformatics/btg363
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP: Gene set enrichment analysis: a knowledgebased approach for interpreting genomewide expression profiles. Proceedings of the National Academy of Sciences of the United States of America 2005, 102(43):15545–15550. 10.1073/pnas.0506580102
Newton MA, Quintana FA, Boon JAD, Sengupta S, Ahlquist P: Randomset Methods Identify Distinct Aspects of the enrichment signal in geneset analysis. Annals of Applied Statistics 2007, 1: 85–106. 10.1214/07AOAS104
Efron B, Tibshirani R: On testing the significance of sets of gene. Annals of Applied Statistics 2007, 1: 107–129. 10.1214/07AOAS101
Ackermann M, Strimmer K: A general modular framework for gene set enrichment analysis. BMC bioinformatics 2009, 10: 47. 10.1186/147121051047
Lönnstedt I, Speed TP: Replicated microarray data. Statistica Sinica 2002, 12: 31–46.
Sartor MA, Tomlinson CR, Wesselkamper SC, Sivaganesan S, Leikauf GD, Medvedovic M: Intensitybased hierarchical Bayes method improves testing for differentially expressed genes in microarray experiments. BMC bioinformatics 2006, 7: 538. 10.1186/147121057538
Goeman JJ, Buhlmann P: Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 2007, 23(8):980–987. 10.1093/bioinformatics/btm051
Chiaretti S, Li X, Gentleman R, Vitale A, Wang KS, Mandelli F, Foa R, Ritz J: Gene expression profiles of Blineage adult acute lymphocytic leukemia reveal genetic patterns that identify lineage derivation and distinct mechanisms of transformation. Clinical cancer research 2005, 11(20):7209–7219. 10.1158/10780432.CCR042165
Nikolsky Y, Ekins S, Nikolskaya T, Bugrim A: A novel method for generation of signature networks as biomarkers from complex high throughput data. Toxicology letters 2005, 158: 20–29. 10.1016/j.toxlet.2005.02.004
Dai M, Wang P, Boyd AD, Kostov G, Athey B, Jones EG, Bunney WE, Myers RM, Speed TP, Akil H, Watson SJ, Meng F: Evolving gene/transcript definitions significantly alter the interpretation of GeneChip data. Nucleic acids research 2005, 33(20):e175. 10.1093/nar/gni179
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, Leisch F, Li C, Maechler M, Rossini AJ, Sawitzki G, Smith C, Smyth G, Tierney L, Yang JY, Zhang J: Bioconductor: open software development for computational biology and bioinformatics. Genome biology 2004, 5(10):R80. 10.1186/gb2004510r80
Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP: Summaries of Affymetrix GeneChip probe level data. Nucleic acids research 2003, 31(4):e15. 10.1093/nar/gng015
Alibes A, Yankilevich P, Canada A, DiazUriarte R: IDconverter and IDClight: conversion and annotation of gene and protein IDs. BMC bioinformatics 2007, 8: 9. 10.1186/1471210589
Barker WC, Garavelli JS, Huang H, McGarvey PB, Orcutt BC, Srinivasarao GY, Xiao C, Yeh LS, Ledley RS, Janda JF, Pfeiffer F, Mewes HW, Tsugita A, Wu C: The protein information resource (PIR). Nucleic acids research 2000, 28: 41–44. 10.1093/nar/28.1.41
Acknowledgements
This work was enabled by the computing facilities of the CSC  IT Center for Science. PT would also like to thank Jarno Tuimala (CSC) and Dario Greco (Helsinki University) for endless R/Bioconductor tutoring and Jarmo Pirhonen (CSC) for Matlab visualization tips. Thanks to Jarkko Salojarvi for pointing Gene Set Analysis. A special thank's to Yrjo Leino (CSC) for series expansion to pvalue calculus. PO has been funded by Ida Montin and Niilo Helander funds, which are gratefully acknowledged.
Author information
Authors and Affiliations
Corresponding author
Additional information
Authors' contributions
PT developed, implemented and tested the methods. PO proposed the research topic and contributed to the datasets and visualization. PM contributed additional mathematical proofs. LH supervised this research. PT, PO and LH evaluated results. All contributed to the writing of the manuscript. All authors read and approved the final manuscript.
Electronic supplementary material
12859_2009_3037_MOESM1_ESM.PDF
Additional file 1: Supplementary text S1. This text shows in detail the stability of the GSZscore and competing methods with randomized data, as the threshold moves through the gene list. In addition, the stability of the methods is monitored as a function of the GO class size. Text also highlights the false positive signals seen with column randomizations. Furthermore, the performance of the used normalization with randomized datasets is shown, and the normalized scores and the empirical pvalues, obtained with row and column randomizations, are compared to each other. The text also represents the detailed comparison of pairwise differences in pvalues for different scoring functions. In addition, most of the Methods are introduced here. Mathematical supplement shows the derivation of E(SN) and D^{2}(SN), required for the derivation of the GSZscore. (PDF 6 MB)
12859_2009_3037_MOESM2_ESM.PDF
Additional file 2: Supplementary text S2. This text represents a mathematical comparison of the GSZscore with hypergeometric Zscore, with correlation between the differential expression scores and GO classification, with maxmean score (by Efron & Tibshiriani) and with Random Sets scoring function by Newton et al. (PDF 481 KB)
12859_2009_3037_MOESM3_ESM.PDF
Additional file 3: Supplementary figure S1: Stability of the GSZscore as the threshold goes through the gene list. Distribution of the GSZscore values as the threshold is moved along the gene list. Subset is smallest at the left and largest (the whole gene list) at the right end of the plot. Results are calculated using all the 4511 GO classes from diabetes dataset with randomized GO class matrix. Blue lines show seven percentiles (0, 5, 25, 50, 75, 95, and 100) at each position. For comparison, the red line shows minimum and maximum scores from the nonrandomized diabetes dataset. Notice the good stability with the regularized GSZscores. Figure is discussed more in a more detail in the supplementary text S1 [see additional file 1]. (PDF 179 KB)
12859_2009_3037_MOESM4_ESM.PDF
Additional file 4: Supplementary figure S2: Stability of KS and modKS as the threshold goes through the gene list. Behaviour of the KS test and the modKS test with the same randomized and positive dataset. Lines represent the same percentiles as in the supplementary figure S1 [see additional file 3] with the same colouring. Notice the biases between different threshold positions, especially when the results are compared with the earlier supplementary figure S1 [see additional file 3]. Notice also the less clear separation between the negative and the positive dataset for modKS in the lower figure. (PDF 341 KB)
12859_2009_3037_MOESM5_ESM.PDF
Additional file 5: Supplementary figure S3: Performance comparison in case (i) with each split analyzed separately. Performance of the methods in each split in case i. Figure represents the AUC score for each evaluated method as the rank limit of the positive GO classes is increased. Note that AUC is calculated here using the whole evaluated GO class list, and it is the size of the used positive GO class set that varies. Methods represented are GSZscore: blue line with circles, ttest: green line with cross, KS test: red line with box, modKS test: cyan line with diamond, iGA: magenta line with x. Notice that although the signal levels vary between the replicates, the differences between the methods are stable. GSZscore and ttest show equal performance among the smallest ranks, whereas the GSZscore is clearly the best among the larger ranks. Other methods show weaker signal. Figure is zoomed to the upper signal area so that most of the curve for the KS is left outside. (PDF 15 KB)
12859_2009_3037_MOESM6_ESM.PDF
Additional file 6: Supplementary figure S4: Performance comparison in case (ii) with each split analyzed separately. Performance of methods in each split in the case ii. The figure represents the AUC score for each evaluated method as the rank limit of the positive GO classes is increased. Methods are coloured identically to the earlier figure. Here, GSZscore and iGA show the best performance at very top ranks, with GSZscore slightly surpassing iGA. Note that performance of ttest and modKS varies considerably across the replicates. After rank 20 GSZscore is constantly best in all the replicates. (PDF 15 KB)
12859_2009_3037_MOESM7_ESM.XLS
Additional file 7: Supplementary table S1: Top 100 GO classes reported by each scoring function from p53 dataset. This table presents 100 best scoring classes for each scoring function from p53 dataset. Table shows class rank, two combined empirical logpvalues, mean of expression values of GO class genes, mean of absolute expression values for GO class genes, class size and class name. In addition, we show four different empirical logpvalues and 7 percentiles (0, 5, 25, 50, 75, 95, 100) for the expression values of the class members. Empirical logpvalues are the ones used to calculate the combined log(pvalue). Biologically relevant (positive) classes are shown in bold font and border cases are with underlined font. Results from each method is presented as a separate Excel worksheet. Furthermore, the first worksheet represents a summary of the results for each scoring function. (XLS 212 KB)
12859_2009_3037_MOESM8_ESM.XLS
Additional file 8: Supplementary table S2: Top 100 GO classes reported by each scoring function from ALL dataset. This table presents 100 best scoring classes for each scoring function from ALL dataset. Columns are similar to the columns in supplementary table S1 [see additional file 7]. Again each method's results are on separate sheet and the first sheet represents the summary. (XLS 210 KB)
12859_2009_3037_MOESM9_ESM.XLS
Additional file 9: Supplementary table S3: GO classes showing largest pairwise differences from ALL dataset. This table presents 50 GO classes with the strongest combined log(pvalue) differences in favour, as well as against the GSZscore, in a pairwise comparison with other scoring functions. Each comparison with each competing scoring function is on a separate sheet. Columns in the tables are identical with columns in tables supplementary table S1 and S2 [see additional files 7 and 8]. (XLS 99 KB)
12859_2009_3037_MOESM10_ESM.XLS
Additional file 10: Supplementary table S4: GO classes showing largest pairwise differences from p53 dataset. This table presents 50 GO classes with the strongest combined log(pvalue) differences in favour, as well as against the GSZscore, in a pairwise comparison with other scoring functions. Each comparison with each competing scoring function is on a separate sheet. Columns in the tables are identical with earlier tables. (XLS 86 KB)
12859_2009_3037_MOESM11_ESM.XLS
Additional file 11: Supplementary table S5: Top GO classes from GSZscore and compared program packages from ALL data. This table presents top classes obtained from GSA, SP, GSEA and GSZscore. Summary of the results is on the first sheet. Separate sheets show each programs results with additional data. Furthermore, we represent two result sheets for scoring functions with one sided test (GSA, GSEA). One presents the one sided test results and the other presents the combined results. (XLS 3 MB)
12859_2009_3037_MOESM12_ESM.XLS
Additional file 12: Supplementary table S6: Top GO classes from GSZscore and compared program packages from p53 data. This table presents top classes obtained from GSA, SP, GSEA and GSZscore. Summary of the results is on the first sheet. Separate sheets show each programs results with additional data. Furthermore, we represent two result sheets for scoring functions with one sided test (GSA, GSEA). One presents the one sided test results and the other presents the combined results. (XLS 1 MB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Törönen, P., Ojala, P.J., Marttinen, P. et al. Robust extraction of functional signals from gene set analysis using a generalized threshold free scoring function. BMC Bioinformatics 10, 307 (2009). https://doi.org/10.1186/1471210510307
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1471210510307