To generate reliable and biologically relevant gene clusters form expression data, we evaluated several available clustering algorithms (e.g. K-means, self-organizing maps) and selected agglomerative hierarchical clustering as the method of choice (Murtagh, 1985; Eisen et al., 1998; de Hoon et al., 2004; R Development Core Team, 2006). The hierarchical clustering method was chosen because of three main advantages: (1) the method requires no prior knowledge about the optimum number of the final clusters, (2) it is extremely robust in joining highly similar items into proper similarity groups, and (3) it provides an information-rich data output that represents the relative distances between all clustered items in a dendrogram (Becker et al., 1988). The main disadvantages of the approach are the complexity of its data output, the lack of predefined boundaries between clusters, and its weaker performance in identifying local expression similarities in a small subset of the samples (Prelic et al., 2006). However, most of these challenges can be overcome by applying efficient postprocessing methods of the obtained dendrograms, such as tree cutting methods (Gutierrez et al., 2007). Popular fuzzy clustering approaches (Krishnapuram et al., 2001) that allow memberships in several clusters—as opposed to strict clustering with unique memberships—were not considered for this study because it is difficult to efficiently prioritize and mine the complex cluster memberships from these methods in the downstream functional analysis steps. As an implementation of the hierarchical clustering algorithm, we used the hclust function (Murtagh, 1985) from the statistical programming environment R (R Development Core Team, 2006). As distance measurement we used correlation coefficients and as the cluster joining method complete linkage (see “Material and Methods” for more details). To obtain discrete clusters from the resulting dendrograms, we developed for this study a novel, to our knowledge, hierarchical threshold clustering (HTC) method. The corresponding R script is available in Supplemental Data S10. This method selects clusters in hierarchical clustering dendrograms based on a maximum tolerable distance between cluster members by applying an all-against-all distance test on all possible subtrees, while maintaining unique cluster memberships. As the threshold we chose for this step a minimum correlation coefficient of 0.6. This relatively conservative HTC setting ensures that all members of any given cluster share with all other members of the same cluster correlation coefficients between the selected cutoff of 0.6 and the highest possible value of 1.0. The exact cutoff value of 0.6 was chosen because it resulted in the highest enrichment of functionally related genes compared to alternative cutoff settings (Supplemental Data S4). Additionally, other gene expression correlation studies have used the same or very similar cutoff values (Haberer et al., 2006; Wei et al., 2006).

Applying the above strategy, we calculated four separate clustering data sets using both the Pearson correlation coefficients (PCC) and the Spearman correlation coefficients (SCC), in their signed and absolute forms as distance measures. The following text will refer to the four methods as PCC, SCC, PCCa, and SCCa, respectively (Supplemental Data S5). All four data sets were generated because of their complementary performance characteristics. The clustering with absolute correlation values allows the identification of positively and negatively correlated gene expressions, whereas the sign-specific approach joins only positively correlated items into similarity groups. The rank-based Spearman approach is limited to identifying global similarities in expression profiles, whereas the Pearson approach is very sensitive in detecting both global and local similarities. In particular, the latter detects local similarities with wide amplitude changes relative to the background, which can result in extreme cases in coclustering of outliers. A consensus approach between several or all methods was not considered because such a strategy would artificially deflate the cluster sizes and compromise the transparency of the results.

The distributions of the obtained numbers of clusters including their sizes from the four clustering methods are summarized in Figure 2. Because the sign removal increases the potential pool sizes of gene pairs with correlation values above a given cutoff, one would expect larger cluster sizes for the data sets with absolute correlation values compared to their signed counterparts. This trend can be observed in the many individual clusters in Supplemental Data S5, but the effect is not very pronounced in the global representation of Figure 2. These relative increases in cluster sizes are not as frequent as expected because of two main reasons. First, the number of highly negatively correlated gene pairs is much smaller than the number of positively correlated gene pairs (data not shown; compare Haberer et al. [2006]). Second, the assignment of a negatively correlated gene to a cluster at an earlier stage of the hierarchical clustering process can prevent other potential members from joining the same cluster at a given cutoff level, if they do not share the required degree of correlation with the existing members. This is particularly the case in combination with a complete linkage joining method that was chosen for this study to minimize the number of false positive members in the generated clusters.

Depending on the stringency of the applied prioritization filters listed in Table IV, our combined clustering and GO term enrichment strategy associated 277 to 1,541 PUF genes to overrepresented GO annotations. In comparison to the GO annotations currently available for these PUF genes, our method associated 216 to 1,050 of them to more specific GO terms in the MF category, 225 to 1,089 in the BP category, and 239 to 1,096 in the CC category (Supplemental Data S6). The large number of PUF genes associated with functionally informative annotations demonstrates the great potential of our approach for guiding future experimental studies on these genes.

Based on enrichment P values, the most highly overrepresented functional categories in the obtained cluster set are the BPs: ribosome assembly, photosynthesis pathways, and cell wall metabolism (Table V). This finding is largely in agreement with related gene coregulation studies in Arabidopsis (Haberer et al., 2006; Wei et al., 2006). With regard to ribosome assembly, 124 of the 410 GO annotated genes for cytosolic, plastidial, and mitochondrial ribosome components appear in seven clusters (see Table V; cluster identifiers 23, 32, 37, 39, 182, 239, and 299); 272 ribosomal genes appear in clusters with five or more members of the nonprioritized data set. Although cluster 23 consists exclusively of genes annotated as ribosomal genes (GO:0005840; P value 1.2 × 10⁻⁶⁴), the other six clusters are highly enriched in ribosomal genes and they contain among others 16 PUF genes. Equally interesting is the observation that photosynthesis-related annotations are highly overrepresented in five large clusters (cluster identifiers 4, 9, 45, 110, and 304). These clusters represent 51 of all the 121 genes that are currently annotated by the GO system as photosynthesis components (GO:0015979). Because both processes, photosynthesis as well as ribosomal activities, require the coordinated assembly of many proteins to large complexes and protein-protein interaction networks, it is not unexpected that their corresponding genes are tightly coregulated. In alignment with the association hypothesis of this study, several of the PUF members in these functionally extremely uniform clusters may be involved in processes that are connected to the enzymatic or regulatory networks of photosynthesis and ribosomal activities.

Interestingly, our method also identified a cluster (identifier 77) that is highly enriched in cell wall-related annotations (e.g. GO:0009834; P value 4.2 × 10⁻¹⁵), such as cellulase synthase genes. A very similar cluster of genes was recently described and experimentally verified by two groups (Brown et al., 2005; Persson et al., 2005) who specifically mined public expression data for genes that are coregulated with the cellulose synthase genes CESA4, CESA7, and CESA8. In addition, comparable results were described by Jen et al. (2006). This example demonstrates that our genome-wide expression clustering approach generates biologically meaningful data. An additional interesting cell wall-related cluster is cluster 349 that contains eight genes for Pro-rich extensin domain proteins.

The majority of the clusters in our data set contain one or more PUF genes (Table IV), but only a few of the larger clusters consist predominantly of PUF genes. Cluster 17 represents an exception to this rule. The 43 members of this cluster contain 26 PUF genes, and its characterized members show no clear enrichment of specific functions. Based on the high abundance of PUF genes in the entire data set (approximately 32%), PUF-gene-enriched clusters occur much less frequent than those enriched in PKF genes; clusters consisting exclusively of PUF genes are entirely absent (Table IV). One explanation for this difference could be that the expression of most PUF genes is controlled by the same regulatory networks as many PKF genes. If this is the case, PUF genes are more likely to appear in expression clusters together with PKF genes than without them.

Our method also identified clusters that are enriched in abiotic stress-response annotations. For instance, clusters 85 and 912 are highly enriched in heat stress-related genes (GO:0009408; P values 6.7 × 10⁻¹⁸ and 1.8 × 10⁻⁶). Interestingly, 10 of the 23 members in the cluster 85 were identified by the subsequent DEG analysis of this study as genes that respond specifically to heat stress and to a much lesser extent to other types of abiotic stresses (see Supplemental Data S7). Based on the available coexpression data, the nine PUF genes of this cluster are now excellent candidates for discovering novel gene functions involved in heat stress-response pathways. Additionally, this example illustrates that the two chosen approaches of this study, expression clustering and DEG analysis, complement and confirm each other. The hypoxia cluster 203 is another interesting abiotic stress cluster (Supplemental Data S6; Fukao and Bailey-Serres, 2004). This cluster does not appear in the most stringently prioritized data set (Table V), because it did not pass the applied uniformity filter. Nevertheless, it is enriched in hypoxia-responsive genes (cluster identifier 203; GO:0001666; P value 2.0 × 10⁻⁵), and it contains several members that are involved in cellular respiration processes, such as genes for the alcohol dehydrogenase ADH1 (AT1G77120), a pyruvate dehydrogenase (AT4G33070), and a hemoglobin-like oxygen binding protein that affects ATP levels under hypoxia (AT2G16060; Hebelstrup et al., 2007). Whether the five PUF genes of this cluster are also involved in hypoxia-response processes can be addressed in experimental studies.

One of the main challenges of performing systematic DEG analyses on large and diverse gene expression data sets from public sources is the identification of the given design parameters to determine for each experiment set its biologically most meaningful analysis strategy. This step is extremely crucial because every analysis needs to focus on the specific treatment factors of an experiment. The alternative of performing simply all possible comparisons will provide meaningless results for many experimental designs because it would generate a large number of illegitimate contrasts between biologically incomparable samples. To define reasonable analysis strategies for public GeneChip microarray expression data sets, all their replicates and the most useful sample comparisons need to be determined manually to provide the proper experimental design parameters to the downstream statistical methods for identifying DEGs. The MIAME and MGED ontology annotations (Brazma et al., 2001; Whetzel et al., 2006) of the public microarray depositories provide the essential information about the experiments, but efficient facilities to completely automate the DEG analyses on a large scale are not available at this point.

To perform large-scale DEG analysis of public expression data, we chose for this study a human-supervised analysis strategy, in which we determined for each experiment set its optimum analysis parameters. The goal of this analysis was to identify all PUF and PKF genes that respond to specific or a wide range of conditions by enumerating their significant expression modulations in the corresponding experiment classes. For this, the available experiment annotations were manually evaluated and the most reasonable set of sample comparisons were recorded in an experiment definition table that contained all the required input parameters to control the downstream statistical DEG analysis in an automated manner (Supplemental Data S2). Typically, we chose for each experiment set a design strategy that focused the analysis on the primary treatment as the main experimental factor. Multifactorial analysis strategies were avoided as much as possible. For instance, when an experiment contained a stress treatment as the primary experimental factor and time or different tissue types as secondary factors, then we compared only samples from identical tissues that were collected at the same time points. Additionally, comparisons between different experiment sets were not considered to exclude unknown variables, such as sample handling differences between laboratories (Hong et al., 2006). It is important to stress here that, depending on the design of a given experiment and its available annotations, it is often difficult to select a single most meaningful analysis strategy. Thus, our chosen strategy may not provide a perfect solution for every experiment set, but it represents a practical and reasonable compromise for performing systematic DEG analyses on large expression data sets from public databases.

In total our large-scale DEG analysis survey included 333 comparisons between samples with two to four technical or biological replicates from 41 experiment sets of six experiment categories. Table III provides an overview of the corresponding sample and experiment sets, and Supplemental Data S2 contains all detailed information including the chosen analysis strategies for these data sets. Because the abiotic stress category is by far the largest data set, containing 524 chip hybridizations of 254 biosamples (Kilian et al., 2007), the following description of our DEG results will be restricted to this most comprehensive treatment category (Table VI). The data for the other categories are provided in the online database of this project (see below). As the statistical method for identifying DEGs with the determined experiment analyses strategies, we used linear models for microarray data (LIMMA) from Smyth (2004, 2005), using in all cases as the confidence threshold a false discovery rate (FDR) of ≤0.01 in combination with a minimum fold-change filter of 2.

Applying the above DEG analysis strategy, we were able to identify 269 PKF and 104 PUF genes that showed expression changes in the majority of the 10 considered abiotic stress categories (Fig. 3; Supplemental Data S7). This set of a total of 373 generic stress DEGs was determined by filtering the generated DEG data set for members that showed one or more significant expression changes in at least 80% of all stress categories. Interestingly, 95% of these DEGs also appear in the generated gene expression clusters of the previous analysis (Supplemental Data S5). The subsequent GO term enrichment analysis revealed that stress-related annotations are highly overrepresented in this group of DEGs (see Supplemental Data S8). About 48 of its members (13%) are associated with the GO term “response to stress” from the BP ontology (GO:0006950; P value 2.0 × 10⁻¹³). This enrichment indicates that our strategy has a high selectivity for identifying stress-response genes. Therefore, many PUF encoding genes in this data set may be directly or indirectly involved in generic stress-response pathways. Among the different groups of identified stress responsive genes (see below; Fig. 3), the generic stress DEG set represents by far the largest group.

Similarly, other studies have shown that stress-regulated genes frequently exhibit expression changes to a wide range of different abiotic stress treatments rather than a refined subset of stresses (Rodriguez and Redman, 2005; Kilian et al., 2007). The group of generic stress DEGs contains 48 genes that are annotated as transcription regulators in the MF ontology (GO:0030528; P value 2.5 × 10⁻³; Supplemental Data S8). This enrichment emphasizes the central role of transcription factors for the control of many stress-response pathways. Moreover, it opens the possibility that several of the 104 PUF genes of this data set may be involved in similar transcription control processes.

We also used the generated abiotic stress DEG data set for identifying genes that respond predominantly to a specific type of stress. These specific stress DEGs were defined as follows. Firstly, they had to show in 25% of all comparisons of a given stress type significant changes. Secondly, they had to exhibit at the same time at least four times as many changes than in the other nine stresses (Fig. 3; Supplemental Data S7). This frequency-based filtering approach appeared to be more efficient for associating DEGs with specific stresses than overly strict filtering methods. This is the case because most stress-response genes are not highly specific for a single type of stress (Kilian et al., 2007). As a result, strict filtering for genes responding only to a single stress will fail to identify any candidate genes in our comprehensive data sets. It is important to emphasize here that the chosen filtering approach is a practical compromise, but not a perfect solution to the problem of assigning DEGs reliably to different stress types. Therefore, the complete DEG results are provided in Supplemental Data S7 where users can apply their own custom filters and prioritize strategies.

With the chosen frequency filter we were able to identify specific stress DEG sets within six of the 10 treatment types (Table VI; Fig. 3). The data sets for the stress treatments—light, oxidative, and wounding stress—did not contain any genes that meet our filtering criteria, and the drought data set contained only a single member. The lack of specific stress DEGs in these data sets indicates that the genome-wide expression response patterns to these four stresses widely overlap with those from other stresses. For the remaining six treatment categories we identified in total 608 PKF and 206 PUF genes that responded predominantly to single stresses. The functional analysis of these specific stress DEG sets with our GO term enrichment approach showed no outstanding enrichment of specific gene functions. Instead, the results contained mostly moderately enriched GO annotations from a wide spectrum of molecular processes and BPs (Supplemental Data S8). Similar to the generic stress data, the different groups of specific stress DEGs included various marker genes that are characteristic for stress-related gene sets. For instance, they contained many genes that are annotated with the GO term “response to stress” (see Fig. 3). This term is significantly enriched in the heat stress data set (P value 1.3 × 10⁻²), whereas the other five treatment sets contain it with considerable, but not significantly enriched frequencies (P values of ≥5 × 10⁻²). In addition, the heat stress and genotoxic stress data sets showed the expected enrichment of genes that are associated with heat response and DNA repair processes, respectively (GO:0009408; P value 4.9 × 10⁻³ and GO:0006281; P value 6.1 × 10⁻⁵).

The normalization of the raw data Cel files was performed in R using the MAS5 and RMA algorithms that are implemented in the affy package form the BioConductor project (Irizarry et al., 2003, 2006; Qin et al., 2006). To allow in the DEG analysis comparisons between the different samples of an experiment set, the RMA normalization was performed in batches for entire experiment sets (Table III). This batch normalization is only required for the quantile-based RMA approach, but not for the MAS5 scaling approach. The present call information of the nonparametric Wilcoxon signed rank test was computed with the affy package to estimate the amount of unexpressed genes (Liu et al., 2002; McClintick and Edenberg, 2006). The obtained expression values from both normalization methods were uploaded to the PED database.

We thank the community projects—R, BioConductor, TAIR, The Institute for Genomic Research, GEO, and AtGenExpress—for providing the excellent software and data resources that were used by this project. T.G. acknowledges support from the Bioinformatics Core Facility, the Center for Plant Cell Biology, and the Institute for Integrative Genome Biology at University of California, Riverside.