# Threshold-free high-power methods for the ontological analysis of genome-wide gene-expression studies

- Björn Nilsson
^{1}Email author, - Petra Håkansson
^{1, 2}, - Mikael Johansson
^{2}, - Sven Nelander
^{3}and - Thoas Fioretos
^{1}

**8**:R74

**DOI: **10.1186/gb-2007-8-5-r74

© Nilsson et al.; licensee BioMed Central Ltd. 2007

**Received: **6 January 2007

**Accepted: **8 May 2007

**Published: **08 May 2007

## Abstract

Ontological analysis facilitates the interpretation of microarray data. Here we describe new ontological analysis methods which, unlike existing approaches, are threshold-free and statistically powerful. We perform extensive evaluations and introduce a new concept, detection spectra, to characterize methods. We show that different ontological analysis methods exhibit distinct detection spectra, and that it is critical to account for this diversity. Our results argue strongly against the continued use of existing methods, and provide directions towards an enhanced approach.

## Background

A fundamental challenge in genome-wide gene-expression studies is to translate complex microarray data into an understanding of the biological conditions being studied. A widely used approach to this problem is ontological analysis - or functional gene-category analysis - the aim of which is to enable data interpretation in the light of known functional relationships between genes. In essence, the methodology seeks to identify categories of functionally associated genes - predefined in external ontologies such as the Gene Ontology Consortium taxonomy [1] - that show deviating expression patterns compared to the general gene population. The underlying motivation is that such categories are presumably likelier to be biologically relevant than gene categories whose expression patterns do not exhibit distinctive features.

Most ontological analysis approaches published so far rely on discrete statistical procedures (binomial, hypergeometric, chi-square or Fisher's exact test) to test for relative enrichments of gene categories within lists of significant genes [2]. These methods are widely used and numerous software packages exist. Nevertheless, discrete methods suffer from a drawback in that the results fundamentally depend on an (essentially arbitrary) threshold for calling genes differentially or non-differentially expressed [3, 4].

To overcome this problem, threshold-free methods for identifying potentially relevant gene categories were recently proposed. Most of these are based on the Kolmogorov-Smirnov (KS) goodness-of-fit test [3–6], although rank-based approaches have also been suggested [7, 8]. The important conceptual advantage of threshold freedom is that the expression data for all genes are considered simultaneously, without the uncertainty associated with previous gene list extraction.

In the study reported here, we enhance the ontological analysis methodology in several important respects. Particularly, we first consider enhanced methods for detecting potentially relevant gene categories. These methods are based on classical and recent examples of a particular class of goodness-of-fit techniques - empirical distribution function (EDF) statistics - that are threshold-free and can be expected to have high statistical power: that is, the chance of detecting a relevant gene category, given it is there, is increased. We carefully assess each method using extensive simulations and by application to multiple real microarray datasets. Second, we develop a new concept, 'detection spectra', which serves to map the prototypic gene categories that are preferentially detected by a given method. We show that different ontological analysis methods exhibit distinct detection spectra, and that it is critical to be aware of this diversity. We also show that, in terms of detection spectrum, the methods represent a continuum ranging from KS on the one extreme to the discrete methods on the other, whereas the remaining methods exhibit intermediate properties. In particular, one method based on the Zhang C (ZC) statistic qualifies as an effective, threshold-free replacement for discrete methods, something that has been previously lacking. Third, to simplify the characterization of detected categories in terms of underlying enrichments of over- or underexpressed genes, we equip each method with an indicator function. These functions indicate the direction of transcriptional deviation, and support the biological interpretation of the ontological analysis results. Finally, we develop a fast significance computation scheme that allows EDF-based analyses to be performed in acceptable time. In conclusion, we introduce attractive alternatives to existing methods for the ontological analysis of microarray experiments, and give directions for the choice of method in practice.

## Results

### Evaluation by simulation

*N*); the proportion of modulated genes (

*π*); and the mean and standard deviation of the modulated gene scores (

*μ*and

*σ*). By varying these parameters, we could artificially recreate gene categories with a broad range of score-distribution dissimilarities.

### Detection spectra

To achieve near-exhaustive testing, we selected 1,800 parameter configurations from wide and relevant intervals, and determined the method powers for each one (see Materials and methods). Hence, for each method, we obtain an 1,800-dimensional performance profile, or detection spectrum, indicating the category types that can be detected.

*σ*= 0.1), intermediate (

*σ*= 0.5) and diffuse (

*σ*= 1.0) effect spreads, respectively. This is consistent with the fact that narrow effects spreads cause discrepancies near the center of the category gene score distribution (KS optimal), whereas intermediate and diffuse spreads lead to dissimilarities which, to greater extents, engage the tails (AD better suited). Fourth, we estimated the coverages of the detection spectra by computing overall (average) powers. The highest values were observed for AD, ZA, ZC, and ZK, implying that these methods are able to detect a broader range of categories than discrete and previous threshold-free (KS-based) methods (Additional data file 2). Taken together, these simulations clearly show that different ontological analysis methods focus on different types of categories, and provide an exact map of the method performances under varying circumstances.

### Method-method relationships

### Application to real data

We proceeded to apply all methods to real microarray data, starting with a dataset from a recent study of ours (P.H., B.N., A Andersson, C Lassen, U Gullberg, and T.F., unpublished work). The aim of this study was to map the transcriptional response of cells to the activity of the fusion oncogene *BCR/ABL1*, associated with chronic myeloid leukemia (CML), in a reverse way by blocking the activity of the fusion protein using the tyrosine-kinase inhibitory drug imatinib mesylate [9]. Essentially, expression profiles of imatinib-treated and non-treated CML cell lines were acquired, and gene scores quantifying the imatinib response were computed (see Materials and methods).

*BCR/ABL1*-mediated leukemogenesis or in the effects of imatinib. For example, consistent with data in the literature, significant enrichments of overexpression were observed in the categories 'heme biosynthesis', and enrichments of underexpression in the 'interferon-gamma signaling pathway', the 'MAPKKK cascade', and in categories related to apoptosis regulation. Also identified was the 'EGF receptor pathway', individual members of which are known to become phosphorylated/activated by

*BCR/ABL1*. Taken together, these findings support the validity of the ontological analysis methodology. Finally, to point at important connections between gene score distributions and detection spectra, and to exemplify the utility of the indicator functions, we selected four illustrative categories, which are discussed in Figure 5.

Functional profile of the imatinib-induced transcriptional response

KS | CM | AD | ZA | ZK | ZC | D1 | D2 | D3 | D4 | D5 | D6 | Δ | Δ | Size* | Biological process (GO) |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

+ | + | -1.00 | -1 | 3 | 'de novo' IMP biosynthesis | ||||||||||

+ | + | + | + | + | + | + | + | + | + | + | +1.00 | +1 | 5 | Heme biosynthesis | |

+ | -1.00 | -1 | 3 | Inactivation of MAPK activity | |||||||||||

+ | + | -1.00 | -1 | 5 | Intracellular transport | ||||||||||

+ | -1.00 | -1 | 5 | Negative regulation of apoptosis | |||||||||||

+ | -0.87 | -1 | 22 | Regulation of transcription | |||||||||||

+ | + | -0.97 | -1 | 7 | Regulation of translational initiation | ||||||||||

+ | + | -1.00 | -1 | 7 | Translational initiation | ||||||||||

KS | CM | AD | ZA | ZK | ZC | D1 | D2 | D3 | D4 | D5 | D6 | Δ | Δ | Size* | Biological process (ABI) |

+ | + | + | +1.00 | +1 | 17 | Hematopoiesis | |||||||||

+ | + | + | + | -1.00 | -1 | 25 | Inhibition of apoptosis | ||||||||

+ | + | + | + | + | -0.89 | -1 | 9 | Macrophage-mediated immunity | |||||||

+ | + | + | + | + | -0.99 | -1 | 37 | MAPKKK cascade | |||||||

+ | + | + | + | -1.00 | -1 | 21 | Protein complex assembly | ||||||||

+ | -1.00 | -1 | 4 | rRNA metabolism | |||||||||||

KS | CM | AD | ZA | ZK | ZC | D1 | D2 | D3 | D4 | D5 | D6 | Δ | Δ | Size* | Molecular function (GO) |

+ | + | + | + | -0.90 | -1 | 257 | ATP binding | ||||||||

+ | -0.61 | -1 | 17 | ATP-dependent helicase activity | |||||||||||

+ | + | + | + | + | -0.91 | -1 | 57 | GTPase activity | |||||||

+ | + | + | + | + | -1.00 | -1 | 3 | Protein kinase C activity | |||||||

+ | -1.00 | -1 | 6 | Protein tyrosine phosphatase activity | |||||||||||

+ | + | + | + | + | + | -0.93 | -1 | 95 | RNA binding | ||||||

+ | + | + | + | + | -0.91 | -1 | 14 | Translation initiation factor activity | |||||||

KS | CM | AD | ZA | ZK | ZC | D1 | D2 | D3 | D4 | D5 | D6 | Δ | Δ | Size* | Molecular function (ABI) |

+ | -0.87 | -1 | 18 | Protein kinase | |||||||||||

KS | CM | AD | ZA | ZK | ZC | D1 | D2 | D3 | D4 | D5 | D6 | Δ | Δ | Size* | Molecular pathway (ABI) |

+ | -1.00 | -1 | 23 | EGR receptor signaling pathway | |||||||||||

+ | + | + | + | + | + | + | + | -1.00 | -1 | 8 | Interferon-gamma signaling pathway | ||||

+ | + | + | + | -1.00 | -1 | 3 | Metabotropic glutamate receptor group I pathway | ||||||||

KS | CM | AD | ZA | ZK | ZC | D1 | D2 | D3 | D4 | D5 | D6 | Δ | Δ | Size* | Cellular component (GO) |

+ | + | + | + | -1.00 | -1 | 147 | Cytoplasm | ||||||||

+ | + | + | + | +1.00 | +1 | 4 | Kinesin complex | ||||||||

+ | + | + | + | + | + | + | +0.73 | +1 | 9 | Microtubule associated complex | |||||

+ | + | + | + | -0.93 | -1 | 23 | Nuclear pore | ||||||||

+ | + | + | -1.00 | -1 | 24 | Nucleolus | |||||||||

+ | + | -1.93 | -1 | 21 | Ribonucleoprotein complex |

### Application to other studies

To verify the generality of our results, we applied all methods to 25 other differential expression comparisons (Additional data file 4) based on seven publicly available microarray datasets (see Materials and methods). While a detailed description of the vast amount of resulting data is beyond the scope of this paper, we point out the following recurrent observations. First, as in Table 1, the choice of method strongly impacted on the results in accordance with the simulated method-method agreements. Second, the numbers of detected categories also approximately followed the simulated method relationships. In broad outline, the threshold-free methods detected many more categories than the discrete methods, because of the noticeable difference in overall power between these two groups. In contrast, the differences within the threshold-free group were less pronounced and more variable, which is explained by the fact that these methods have more similar overall powers, implying that the number of detected categories will, to a greater extent, be determined by the match between the detection spectrum and the set of deviating categories that are actually present in the data. Taken together, these findings further underscore the fact that different methods focus on different category types, and, hence, that it is important to be aware of this in practice.

### Computational efficiency

The total time required for analyzing the 25 studies (all methods and ontologies) was 32 seconds (C++ implementation; 2 GHz Core2Duo PC), illustrating the benefit of the fast significance computation scheme (see Materials and methods).

## Discussion

The ontological analysis of genome-wide studies relies fundamentally on the validity and continued growth of ontologies providing annotations of gene function. However, efficient computational methods are needed to integrate these annotations with data in an optimal way. We have addressed the latter problem by considering gene-category identification methods based on high-power EDF statistics.

We have shown that the value of these methods lies in their higher overall powers - implying an ability to detect a broader range of potentially biologically relevant gene categories - and in their detection spectra, which are distinct from those of existing methods. Previously, KS-based and discrete approaches have focused on high-proportion-low-effects and low-proportion-high-effects deviations, respectively [4], whereas methods with intermediate detection spectra and threshold-free methods for detecting low-proportion-high-effects deviations have been lacking. The methods described fill these gaps. In particular, our data suggest ZC to be a new method of choice for low-proportion-high-effects-oriented analysis. Offering excellent low-proportion-high-effects coverage, high overall power and the obvious advantages of threshold freedom, ZC virtually removes the need for discrete methods. Regarding the remaining methods, ZK is slightly less tail oriented than ZC; ZA and AD focus on intermediate-proportion- moderate-effects categories; Finally, CM resembles KS.

As shown, the choice of category-detection method has a profound impact on the results of the ontological analysis. However, the question of what prototypic categories are most biologically relevant, and hence should be the primary target in ontological analyses, is an open problem. In the absence of solid evidence supporting that one category type is generally more biologically relevant than others, the choice of method must be partly guided by the investigator's preferences and project-specific considerations. For example, in data containing strongly differentially expressed genes (for example, the imatinib study presented here), it may be natural to optimize the analysis for low-proportion-high-effects categories, making the tail-sensitive methods (ZC, ZK and discrete) the methods of choice. On the other hand, in datasets where differentially expressed genes display predominantly low-to-moderate effect sizes, it seems more reasonable to focus on intermediate-proportion-moderate-effects and high-proportion-low-effects categories, motivating the choice of methods with more center-oriented detection spectra in this case.

Alternatively, multiple methods can be used in concert, provided that appropriate statistical corrections are made. Such an approach would yield results tables similar to Table 1, and offers the benefit of allowing the user to make indirect conclusions about the characteristics of the distributional deviations of the detected categories. For example, if a category is detected as significant by all methods (for example, the heme biosynthesis category in the imatinib experiment), then, quite clearly, its gene-score distribution must be highly aberrant, most probably because of a high-proportion-high-effects size enrichment. In contrast, if a category is called significant by one method (for example, AD with the EGF receptor signaling pathway in the imatinib data), and not by the others, then the distributional deviations must fall within the detection spectra of that method but outside the detection spectra of the other methods. In the EGF receptor signaling pathway example, a reasonable conclusion - given the detection spectra established in Additional data file 1 - would be that underlying deviation is likely to be of intermediate-proportion-moderate-effects-size type, as such categories represent the detection optimum of AD.

Regardless of category-detection statistic, the reference gene population, null model, gene score, and ontology also influence the results and must be chosen judiciously [2, 3, 7]. We recognize that a shared limitation of many ontological analysis methods, ours included, is that dependencies between genes are not taken into account when computing significances, something that may lead to underestimated *p* values. First steps have been taken to develop dependency-modeling schemes, for example SAFE [3] or CatMap [7]. While the methods described can be adopted into those frameworks if desired, additional efforts are needed to address the problem of modeling dependencies in detail.

Other features introduced are indicator functions and a fast significance computation scheme. The indicator functions facilitate the interpretation of the results of the ontological analysis. Their advantage compared to existing approaches is that the need for two separate tests per category, one to detect enrichment of overexpression and one to detect enrichment of underexpression, is removed. A limitation is that enrichments cannot be distinguished from (contralateral) depletions. Such ambiguities can be resolved graphically (Figure 5), or through the development of improved versions in future studies. The fast significance computations are not crucial to the ontological analysis as such, but are valuable in that they allow the procedure to be performed within an acceptable time-frame.

Finally, we recognize the limitations of the gene-category model used for computing the detection spectra. First, as already discussed, the model assumes independence between genes. Second, for tractability, we have limited our treatment to gene categories with only one group of the modulated genes. While the model could be extended to multiple modulated groups, this would obviously increase the complexity of the study at the expense of presentational clarity and understandability. Third, we have only considered gene scores with approximately normal distributions, which is a minor limitation as the most frequently used gene scores are based on *t*-statistics. Nevertheless, the properties of the described methods for scores with distinctly different distributions (for example, scores based on the *F*-statistic) remain to be established.

## Conclusion

We have presented novel ontological analysis methods constituting attractive alternatives to existing approaches. Hence, this work contributes to the repertoire of useful methods aiding the interpretation of genome-wide gene expression studies.

## Materials and methods

### Threshold-free category-detection methods

Let *F*_{
N
}(*x*) and *F*(*x*) denote the empirical (cumulative) distribution functions for the gene-specific differential expression scores *x*_{1},...,*x*_{
N
}for an *N*-gene category and the scores *x*'_{1},...,*x*'_{
M
}for an *M*-gene reference gene population, respectively. We consider six EDF statistics to measure discrepancy between *F*_{
N
}(*x*) and *F*(*x*), that is, to detect gene categories with deviating gene-score distributions. A technicality that arises is that, normally when using EDF statistics, *F*(*x*) is specified by a continuous function. This is obviously not the case here as *F*(*x*) is an EDF, jumping by 1/*M* at each *x*'_{
i
}. However, we note that, in this application, this issue can be ignored because *M* is large (on the order of 5,000 to 40,000), making *F*(*x*) sufficiently smooth to be regarded as continuous.

*F*(

*x*) and

*F*

_{ N }(

*x*)

*y*

_{ i }=

*F*(

*x*

_{ i }). While intuitively straightforward and capable of detecting discrepancies near the center of the distribution, KS fails to notice subtle discrepancies in the tails as well as small but consistent deviations. Second and third, we use two members of the Cramér-von Mises family of quadratic EDF statistics, defined by

*ψ*(

*x*) is a suitable weight function. We consider

*ψ*(

*x*) = 1, which generates the Cramér-von Mises (CM) statistic itself [12]

*ψ*(

*x*) =

*F*(

*x*)

^{-1}(1-

*F*(

*x*))

^{-1}which gives more weight to the tails and generates the Anderson-Darling (AD) statistic [13]

### Discrete category detection methods

For comparison, we also included a discrete category-detection method in the study. As a representative of this class of methods, we used the binomial test, which is routinely used as an approximation to hypergeometric procedures, such as Fisher's exact test, when the population is large. The technical details can be found in standard statistics textbooks or in work on ontological analysis (see [2] and references therein). Throughout, we used six thresholds for calling genes differentially expressed (1.5, 1.8 to 3.0; denoted D1 to D6) and used gene scores that would always be compatible with these values (see Microarray datasets).

### Category characterization methods

While the EDF statistics effectively detect deviating - and thus presumably biologically relevant - gene categories, they do not, in their basic form, indicate whether the deviations are caused by enrichments of overexpressed genes, underexpressed genes, or a mixture of both. In previous approaches, this problem has been addressed by performing two separate tests for each category, one to detect enrichments of overexpression and one to detect enrichments of underexpression. Here, we proceed differently and instead derive an indicator for each detection method. The advantage of these functions is that the direction of transcriptional deviation is determined in a continuous manner, removing the need for double testing.

*δ*

_{ i }are defined separately for each statistic and depend on

*x*

_{ i }and

*y*

_{ i }(calculations not shown). We let the indicators Δ be the raw correlation between $\delta ={\left\{{\delta}_{i}\right\}}_{1}^{N}$ and ${\delta}^{\prime}={\{sign({F}_{N}({x}_{i})-F({x}_{i}))\cdot {\delta}_{i}\}}_{1}^{N}$, that is,

When the denominator is zero, we let Δ = 0. For KS and ZK, which are based on max operators instead of sums, we let

Δ = *sign*(*F*_{
N
}(*x*_{
i
}) - *F*(*x*_{
i
})),

where *i* is the index used when computing the KS or ZK statistic. The Δ indicators characterize gene categories by considering the distributional dissimilarities that led to their detection. In categories with unexpectedly many overexpressed genes, we have *F*_{
N
}(*x*_{
i
}) ≤ *F*(*x*_{
i
}) for all *i*, implying Δ = 1. Conversely, categories with unexpectedly many underexpressed genes, will receive Δ = -1. Moreover, for CM, AD, ZA and ZC, Δ will attain intermediate values depending on the balance between the two types of genes. The KS and ZK indicators are less informative, evaluating to either -1 or 1 depending on the predominant direction of deviation.

### Significance computations

The null distributions for the EDF statistics are unknown, and, in some cases (ZA, ZC and ZK), asymptotic theory is lacking. To compute significances, we therefore used a procedure based on Monte Carlo simulation by gene permutations, which is currently a standard scheme in ontological analysis although it does not account for dependencies between genes. More elaborate schemes seeking to model dependencies using sample label permutations have been suggested [3, 7], and the methods above can be adopted into those frameworks if needed.

In principle, the null distributions could be simulated from scratch for every category. However, that approach turned out to be exceedingly time-consuming. We instead note that the assumed continuity of *F*(*x*) implies that the EDF statistics are distribution-free. Hence, their null distributions can be pre-computed by drawing *y*_{
i
}'s from a uniform distribution, a procedure that is essentially equivalent to permuting genes when the population is large (assumed). This strategy completely avoids simulations at runtime, allowing entire ontological analyses to be performed in instants. Throughout, the distributions were pre-simulated using 10^{8} Monte Carlo replicates (per category size and statistic), and were compressed to tractable sizes using a recent algorithm (B.N. unpublished work).

### Simulation model

To simulate the reference gene score population, we drew 10,000 scores from a standard normal distribution (zero mean, unit variance). This choice is motivated by the fact that differential expression is frequently assessed using the *t*-statistic or variance-moderated versions thereof [15–17], in which cases the population scores will be approximately normally distributed as most genes are non-differentially expressed. Furthermore, to simulate deviating gene categories, we used the mixture model previously proposed in [4], in which a proportion of the category genes are given scores from a modulated normal distribution (non-zero mean, non-unit variance) whereas the remaining genes are given scores from a standard normal like the reference population (Figure 1). The model parameters are: the number of category genes (*N*), the proportion of modulated genes (*π*), the mean (effect size) of the modulated gene scores (*μ*), and the standard error (effects spread) of the modulated gene scores (*σ*). The parameter values were: *N* = 10, 30 and 100 genes, which are typical category sizes; *π* = 0.1, 0.2 to 1.0, which is essentially exhaustive; *μ* = 0.2, 0.4 to 4.0, covering very weak to very strong effects. Because the EDF statistics are distribution-symmetric, negative and positive *μ* values will yield identical results. Hence, the evaluation can be restricted to positive values without loss of generality. Finally, *σ* = 0.1, 0.5 and 1.0, corresponding to narrow, intermediate and diffuse effects spreads, respectively. Thus, the total number of four-parameter combinations was 3 × 10 × 20 × 3 = 1,800. For each combination, 100,000 random categories were generated and tested for conformity with the population distribution. The parameter-configuration-specific statistical powers were estimated as the proportions of categories called significant at the *p* < 0.001 level (the full set of raw data is in Additional data file 5). To verify robustness, the experiments were repeated with numerous other cutoff levels, yielding results in broad agreement with those presented.

To quantify the diversity of category types detected, we computed overall (average) powers across all parameter configurations and across *π* and *μ* for fixed *N* and *σ*. To quantify method-method agreements, we computed the Spearman rank correlation and the Jaccard similarity coefficient (or Jaccard index) for all pairs of methods. The Spearman metric, the correlation between the rank-transformed *p* values, measures similarity between category rankings. The Jaccard similarity coefficient, the proportion of categories called significant by both methods, reflects similarity between sets of detected categories.

### Microarray datasets

The imatinib data (P.H., B.N., A Andersson, C Lassen, U Gullberg, and T.F. unpublished work) were generated at our lab by culturing five CML cell lines in the presence or absence of imatinib mesylate. Expression profiles were obtained at 3 and 12 h after drug exposure using 27 K cDNA arrays. For each time point and treatment group, two technical replicates were obtained, yielding 2 × 2 × 2 = 8 arrays per cell line. After filtering and probe merging, 5,532 unique Entrez Gene entries remained. The full set of microarray data will be made available upon acceptance of the original work.

In addition to the dataset from the imatinib experiment, we included seven publicly available expression array datasets. The Valk dataset [18] and the Radich dataset [19] were obtained from the NCBI Gene Expression Omnibus repository [20], accessions GSE1159 and GSE4170, respectively. The Zheng dataset [21] was obtained from the ArrayExpress repository [22]. The Bhattacharjee dataset [23] was obtained from the Broad Institute website [24]. The Ross dataset [25] was obtained from the St Jude Children's Research Hospital website [26]. The Andersson dataset [27] was obtained by personal communication with the corresponding author. The West dataset [28] was obtained from the Duke University website [29].

As a score of differential expression, we used Smyth's moderated *t*-statistic [15], which follows an approximate *t*-distribution under the null hypothesis whenever the data are reasonably normal. Hence, in any given study, the variance of the population scores will be near one, guaranteeing that the thresholds used with the discrete method are meaningful.

### Ontologies

A total of six ontologies from the Gene Ontology (GO) Consortium [30] and the Applied Biosystems Panther Gene Classification System (ABI) [31] were used: Biological Process (GO+ABI), Molecular Function (GO+ABI), Cellular Component (GO), and Molecular Pathway (ABI). The ontology versions used in the analyses were those available in December 2006.

### Software availability

To allow readers to readily apply the described methods to their own data, we provide a software package called RenderCat (stand-alone Windows executable). This software is publicly and freely available on request from B.N. The source code is open and can be downloaded from the SourceForge repository [32]. The package includes implementations of all the category-detection methods described, including the indicator functions and the fast significance computations. For the convenience of the user, we have also included functionality for creating tables similar to Table 1 (tab-delimited text or LaTeX format) and capability for rendering gene-category score-distribution plots similar to Figure 5 (bitmap format). To correct for multiple testing, the program uses the false-discovery rate [33, 34].

## Additional data files

Additional data are available online with this paper. Additional data file 1 is a figure representing the complete results of the simulation study. Additional data file 2 is a table listing overall powers. Additional data file 3 is a table containing the complete data from the method-method agreement assessment study. Additional data file 4 contains a list of additional differential expression studies. Additional data file 5 contains the raw data used for generating the detection spectra.

## Declarations

### Acknowledgements

This work was supported by research grants from the Swedish Cancer Society, the Swedish Children's Cancer Foundation, and the Medical Faculty of Lund University. We thank Jill Storry for superb proofreading of the manuscript.

## Authors’ Affiliations

## References

- Ashburner M, Ball C, Blake J, Botstein D, Butler H, Cherry J, Davis A, Dolinski K, Dwight S, Eppig J, et al: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000, 25: 25-29. 10.1038/75556.PubMedPubMed CentralView ArticleGoogle Scholar
- Khatri P, Draghici S: Ontological analysis of gene expression data: current tools, limitations, and open problems. Bioinformatics. 2005, 21: 3587-3595. 10.1093/bioinformatics/bti565.PubMedPubMed CentralView ArticleGoogle Scholar
- Barry W, Nobel A, Wright F: Significance analysis of functional categories in gene expression studies: a structured permutation approach. Bioinformatics. 2005, 21: 1943-1949. 10.1093/bioinformatics/bti260.PubMedView ArticleGoogle Scholar
- Ben-Shaul Y, Bergman H, Soreq H: Identifying subtle interrelated changes in functional gene categories using continuous measures of gene expression. Bioinformatics. 2005, 21: 1129-1137. 10.1093/bioinformatics/bti149.PubMedView ArticleGoogle Scholar
- Mootha V, Lindgren C, Eriksson K, Subramanian A, Sihag S, Lehar J, Puigserver P, Carlsson E, Ridderstrale M, Laurila E, et al: PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003, 34: 267-273. 10.1038/ng1180.PubMedView ArticleGoogle Scholar
- Lamb J, Ramaswamy S, Ford H, Contreras B, Martinez R, Kittrell F, Zahnow C, Patterson N, Golub T, Ewen M: A mechanism of cyclin D1 action encoded in the patterns of gene expression in human cancer. Cell. 2003, 114: 323-334. 10.1016/S0092-8674(03)00570-1.PubMedView ArticleGoogle Scholar
- Breslin T, Eden P, Krogh M: Comparing functional annotation analyses with Catmap. BMC Bioinformatics. 2004, 5: 193-10.1186/1471-2105-5-193.PubMedPubMed CentralView ArticleGoogle Scholar
- Lee H, Braynen W, Keshav K, Pavlidis P: ErmineJ: tool for functional analysis of gene expression data sets. BMC Bioinformatics. 2005, 9: 269-10.1186/1471-2105-6-269.View ArticleGoogle Scholar
- Deininger M, Buchdunger E, Druker B: The development of imatinib as a therapeutic agent for chronic myeloid leukemia. Blood. 2005, 105: 2640-2653. 10.1182/blood-2004-08-3097.PubMedView ArticleGoogle Scholar
- Kolmogorov A: Sulla determinazione empirica di una legge di distibuziane. Giorna Ist Attuari. 1933, 4: 83-91.Google Scholar
- Smirnov N: Estimate of deviation between empirical distribution functions in two independent samples. Bull Mosk Univ. 1939, 2: 3-16.Google Scholar
- Cramér H: On the composition of elementary errors: II, Statistical applications. Skand Akt. 1928, 11: 141-180.Google Scholar
- Anderson T, Darling D: A test of goodness of fit. J Am Stat Ass. 1954, 49: 765-769. 10.2307/2281537.View ArticleGoogle Scholar
- Zhang J: Powerful goodness-of-fit tests based on the likelihood ratio. J R Stat Soc B. 2002, 64: 281-294. 10.1111/1467-9868.00337.View ArticleGoogle Scholar
- Smyth G: Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004, 3: Article 3-Google Scholar
- Cui X, Hwang JG, Qiu J, Blades N, Churchhill G: Improved statistical tests for differential gene expression by shrinking variance components estimates. Biostatistics. 2005, 6: 59-75. 10.1093/biostatistics/kxh018.PubMedView ArticleGoogle Scholar
- Storey J, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci USA. 2003, 100: 9440-9445. 10.1073/pnas.1530509100.PubMedPubMed CentralView ArticleGoogle Scholar
- Valk P, Verhaak R, Beijen M, Erpelinck C, van Waalwijk , van Doorn-Khosrovani SB, Boer J, Beverloo H, Moorhouse M, van der Spek P, Lowenberg B, Delwel R: Prognostically useful gene-expression profiles in acute myeloid leukemia. N Engl J Med. 2004, 350: 1617-1628. 10.1056/NEJMoa040465.PubMedView ArticleGoogle Scholar
- Radich J, Dai H, Mao M, Oehler V, Schelter J, Druker B, Sawyers C, Shah N, Stock W, Willman C, Friend S, Lindsey P: Gene expression changes associated with progression and response in chronic myeloid leukemia. Proc Natl Acad Sci USA. 2006, 103: 2794-2799. 10.1073/pnas.0510423103.PubMedPubMed CentralView ArticleGoogle Scholar
- Gene Expression Omnibus repository. [http://www.ncbi.nlm.nih.gov/geo/]
- Zheng C, Li L, Haak M, Brors B, Frank O, Giehl M, Fabarius A, Schatz M, Weisser A, Lorentz C, Gretz N, Hehlmann R, Hochhaus A, Seifarth W: Gene expression profiling of CD34
^{+}cells identifies a molecular signature of chronic myeloid leukemia blast crisis. Leukemia. 2006, 20: 1028-1034. 10.1038/sj.leu.2404227.PubMedView ArticleGoogle Scholar - ArrayExpress repository. [http://www.ebi.ac.uk/arrayexpress]
- Bhattacharjee A, Richards W, Staunton J, Li C, Monti S, Vasa P, Ladd C, Beheshti J, Bueno R, Gillette M, et al: Classification of human lung carcinomas by mRNA expression profiling reveals distinct adenocarcinoma subclasses. Proc Natl Acad Sci USA. 2001, 98: 13790-13795. 10.1073/pnas.191502998.PubMedPubMed CentralView ArticleGoogle Scholar
- Broad Institute: cancer genomics publications. [http://www.broad.mit.edu/mpr/lung]
- Ross M, Zhou X, Song G, Shurtleff S, Girtman K, Williams W, Liu H, Mahfouz R, Raimondi S, Lenny N, Patel A, Downing J: Classification of pediatric acute lymphoblastic leukemia by gene expression profiling. Blood. 2003, 102: 2951-2959. 10.1182/blood-2003-01-0338.PubMedView ArticleGoogle Scholar
- St Jude Research. [http://www.stjuderesearch.org]
- Andersson A, Olofsson T, Lindgren D, Nilsson B, Ritz C, Eden P, Lassen C, Rade J, Fontes M, Morse H, et al: Molecular signatures in childhood acute leukemia and their correlations to expression patterns in normal hematopoietic subpopulations. Proc Natl Acad Sci USA. 2005, 102: 19069-19074. 10.1073/pnas.0506637102.PubMedPubMed CentralView ArticleGoogle Scholar
- West M, Blanchette C, Dressman H, Huang E, Ishida S, Spang R, Zuzan H, Olson J, Marks J, Nevins J: Predicting the clinical status of human breast cancer by using gene expression profiles. Proc Natl Acad Sci USA. 2001, 98: 11462-11467. 10.1073/pnas.201162998.PubMedPubMed CentralView ArticleGoogle Scholar
- Duke University IGSP - supplemental data. [http://data.cgt.duke.edu/west.php]
- The Gene Ontology. [http://www.geneontology.org]
- PANTHER - classification of genes and proteins. [http://www.pantherdb.org]
- SourceForge.net: RenderCat. [http://sourceforge.net/projects/rendercat]
- Bejamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995, 57: 289-300.Google Scholar
- Bejamini Y, Yekutieli D: The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001, 29: 1165-1188. 10.1214/aos/1013699998.View ArticleGoogle Scholar

## Copyright

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.