Mining gene expression data by interpreting principal components

Background There are many methods for analyzing microarray data that group together genes having similar patterns of expression over all conditions tested. However, in many instances the biologically important goal is to identify relatively small sets of genes that share coherent expression across only some conditions, rather than all or most conditions as required in traditional clustering; e.g. genes that are highly up-regulated and/or down-regulated similarly across only a subset of conditions. Equally important is the need to learn which conditions are the decisive ones in forming such gene sets of interest, and how they relate to diverse conditional covariates, such as disease diagnosis or prognosis. Results We present a method for automatically identifying such candidate sets of biologically relevant genes using a combination of principal components analysis and information theoretic metrics. To enable easy use of our methods, we have developed a data analysis package that facilitates visualization and subsequent data mining of the independent sources of significant variation present in gene microarray expression datasets (or in any other similarly structured high-dimensional dataset). We applied these tools to two public datasets, and highlight sets of genes most affected by specific subsets of conditions (e.g. tissues, treatments, samples, etc.). Statistically significant associations for highlighted gene sets were shown via global analysis for Gene Ontology term enrichment. Together with covariate associations, the tool provides a basis for building testable hypotheses about the biological or experimental causes of observed variation. Conclusion We provide an unsupervised data mining technique for diverse microarray expression datasets that is distinct from major methods now in routine use. In test uses, this method, based on publicly available gene annotations, appears to identify numerous sets of biologically relevant genes. It has proven especially valuable in instances where there are many diverse conditions (10's to hundreds of different tissues or cell types), a situation in which many clustering and ordering algorithms become problematic. This approach also shows promise in other topic domains such as multi-spectral imaging datasets.


Background
Bioinformatics has placed much emphasis on using various unsupervised clustering techniques as a means to understand the information present in gene microarray expression datasets. Clustering techniques produce a rich taxonomy of results by defining groups of genes that act more or less similarly across a number of experimental conditions. The diverse approaches to clustering genes by expression levels include k-means [1], self-organizing maps [2], hierarchical algorithms [3,4] and probabilistic models [5]. Some approaches permit clustering of the conditions as well [6][7][8]. Based on co-expression, genes that comprise individual expression clusters are often postulated to be co-regulated, and to the extent that this hypothesis is correct in any specific biological situation, the gene cluster definitions can offer key insights into gene regulatory network (GRN) structure and function.
Another common data mining task is to try to identify small sets of genes that can serve as effective predictors of disease diagnosis or prognosis. While clustering at its best is good at finding sets of genes that are similarly expressed across all conditions within a dataset, many issues (e.g. selection of K, stochastic effects, and "noise" from large numbers of genes that change little over most of the conditions) can prevent clustering from successfully highlighting small groups of interestingly co-expressed genes [9,10]. This often encountered problem is addressed in part by hierarchical phylogenetic ordering algorithms such as average linkage in Mike Eisen's cluster program [3], but the information biologists seek regarding shared sub-patterns of co-expression can be obscured by both algorithmic and visualization constraints. The algorithmic limitations in hierarchical clustering confound and "cover" the presence and organization of smaller and more specific gene groups that are similar across only a subset of conditions within the larger dataset. In any case, biologists generally subjectively define a cluster of genes from such phylogenetic trees based largely on human pattern recognition. Finally, nothing inherent in the clustering approach helps to direct a biologist to which cluster is interesting or relevant. Instead, biologists generally take the path of focussing on a group of genes exhibiting a pattern of expression that supports a specific hypothesis, or search for a known gene or genes of interest within a cluster to form an explanation for others in the cluster. Support vector machines have been shown to be useful for identifying small sets of related and predictive genes [11][12][13][14], but represent a supervised learning approach which requires one to first define a set of known positive examples, a set of known negative examples, and a specific covariate to predict. We wanted an unsupervised algorithm that would help us to find relationships and structure in the data that is more specific than what clustering algo-rithms usually deliver, yet is hypothesis independent. We found it efficient and useful to use as an independent starting tool a very direct approach based on principal components analysis (PCA, see Methods section). A virtue is that this approach is computationally efficient for very large datasets, especially compared with most clustering algorithms, but is also applicable to much smaller ones. It allows one to directly explore each of the independent and diverse sources of variation present within a gene expression dataset and to subsequently identify the specific genes that vary the most, together with the conditions in which they vary.
Unlike conventional clustering and ordering algorithms, this PCA based approach permits a gene to be highlighted and grouped as influential in multiple condition sets, whereas in cluster membership a gene is typically assigned to one unique cluster. The "single cluster assignment" quality of traditional clustering and ordering algorithms is problematical because it tends to hide commonality of expression that is restricted to a small, interesting, and often entirely unpredicted subgroups of tissues, cell types, treatments or other condition types. This situation, perhaps because of inherent properties of gene network structure, will arise increasingly as the number and diversity of conditions represented in expression datasets increases.
The clustering method of Barkai et al. [15,16] addresses this issue of multiple membership in a different way, by using randomly-selected gene sets to iteratively search for and refine self-consistent groups. Their approach, which is related to PCA through singular value decomposition (SVD), also permits genes to be assigned to multiple "expression modules." In contrast to the method presented here, there is no provision for correlating modules with covariate data.
The use of principal components analysis presented here differs from other recent applications in gene expression analysis. PCA is most commonly used in as a means of dimensionality reduction prior to clustering [7,17] or prior to classification [18,19]. It is also used to visualize or confirm clustering results [19][20][21]. In contrast, our use of PCA aims to find, then examine, and where possible, generate hypotheses to explain individual principal components. In this manner, we build on the observations by Hilsenbeck [22] Raychaudhuri et. al. [23] who used PCA to gain insight into the underlying factors present in the Chu et. al. yeast sporulation experiments [24]. Wall et. al. [25] introduced a novel use of singular value decomposition (SVD) for gene expression analysis that identifies non-exclusive gene groups, and Selaru et. al. [26] illustrated the potential of PCA to detect molecular phenotypic bases that correspond to relevant clinical or biological features of human tumors. Their approach identifies a subset of principal components that correlate well with known covariates. Here we introduce methods that extend beyond producing gene groups and observing a few principal components. Our methods provide a path for systematically analyzing each principal component by identifying the genes most influential in defining a particular principal component and the conditions in which those influential genes vary significantly. Finally, we describe methods which aim to explain each principal component's observed variance in terms of the condition variables deemed most likely to be driving the variance. We also introduce a software package that implements these methods. These methods and our software package provide automated and objective way of doing what a biologist naturally tries to do through inspection and pattern recognition.

Implementation
We have developed a Python package to implement the PCA interpretation capability described in detail in the Methodology section. This PCA analysis package has been added to CompClust developed previously [9,27]. The combined packages allow one to cluster, classify and visualize numeric datasets that have discrete or numeric annotations (referred to as labelings, or labelled datasets), and to compare labelings with confusion matrices and metrics such as normalized mutual information (NMI) [28]. This PCA analysis tool (including the complete results for the dataset analysis described in the Results section) has also been made accessible through the CompClustWeb webbased interface [29]. Our software makes use of data manipulation and graphical plotting using the matplotlib package [30], and the statistics are generated using the rPy package [31] and Gary Strangman's Python stats package [32].
The web-based front-end permits users to get a complete report on the interpretation of each principal component, including interactive PCA projection plots with the principal component's extreme genes (PCEGs) highlighted; ranked lists of the PCEGs with detailed annotations; interactive significance-ordered gene set trajectory plots that permit users to drill down to the individual gene level; similarly ordered condition reports ordered by expression difference and grouped by significance, including covariate info (with significantly correlated covariates highlighted); and finally a report of any suggestive covariates that are well correlated with the significant column grouping, including the confusion matrices and/or plots of statistics scores to back up the conclusions. All principal component analysis and results generation is implemented in a Python package so that analyses of large datasets can be executed in a batch mode rather than through the graphical interface. Further, the software that implements the CompClustWeb interface is provided within the CompClust package, so a software developer can create his or her own CompClustWeb server to review results of their PCA interpretation.

Application to microarray expression data, Case 1: GNF human data
We obtained gene expression data from the Genomics Institute of the Novartis Research Foundation ("GNF") Gene Expression Database via their SymAtlas web site [33,34]. The dataset is a challenge for most clustering algorithms because it contains 158 tissue samples hybridized to two Affymetrix microarray chips: U133A and GNF1H. The dataset combines the measurements of these chips to provide a total of 33,689 unique probe identities across the 158 tissue samples. Expression data are signal intensities estimated by Affymetrix Microarray Suite v5. For our analysis we used the log base 2 of the expression signal, and included data for all tissues and probes (noting that absent and present calls were not provided with the signal intensities). We applied our principal components analysis tool to generate interpretations for each of this dataset's 158 principal components.
As detailed in the Methodology section, for each principal component we identified a set of gene probes occupying the high and low extremes of that principal component's axis (we refer to these as principal component extreme genes, or PCEGs). One can adjust parameters to recover smaller and larger numbers of PCEGs per component by specifying either a likelihood threshold or an explicit number of PCEGs. The PCEGs are those probes having the most highly weighted values for that principal component, selected because they stand out from the others, they influence the principal component's direction, and thus they warrant further investigation. We selected probes with likelihoods less than extremeThresh = 0.00001, which yielded on average 20 low and 20 high extreme genes per principal component, though the sets sizes do vary considerably (µ = 18.9, σ = 17.2). Next we identified the tissues in which the high PCEGs showed significantly different expression than the low PCEGs. Visualizations produced include scatter plots of the extreme genes in PCA sub-spaces (PC N-1vs. PC N), and extreme gene trajectories in original tissue order as well as with tissues ordered by decreasing difference of mean of high PCEGs and mean of low PCEGs. The latter trajectory plot emphasizes how the extreme genes for a principal component show a pattern of expression that imposes a partitioning of tissues. It is left to human interpretation to examine the extreme genes and the tissue partitioning exposed by each principal component, and thus to build hypotheses that attach meaning to the sources of variation. The percentage variance explained by the top 50 principal components is provided in Table 1. Example results for two illustrative  principal components, PC7 and PC21, are shown in Figures 1 through 4 and in Tables 2 through 5, and are [35] contains the complete collection of PCA interpretation results generated for these dataset for all principal components, as well as results of a comparable analysis done at extremeThresh = 0.001 which yielded larger PCEG sets.
We addressed the question of biological and statistical significance of PCs and the sets of extreme genes identified. Each set of high and low extreme genes from each principal component was tested for Gene Ontology (GO) statistical enrichment when compared to the human GO annotations from NCBI's loc2go dataset using the hypergeometric to calculate the p-value of each GO term. Terms that were enriched for a particular PC at 1% significance threshold and that were still significant following a Bonferroni correction for multiple hypothesis testing as described in [36] are reported as enriched (see Table 6). 26 of the top 42 PCEG lists, derived using a stringent cutoff of 0.00001, produced significant GO enrichments; no PCEG sets beyond PC42 showed significant enrichment. As discussed below, many of the significant results showed obvious biological coherence and relationships to the specific samples associated with the PC of origin. This argues that PCs containing less than 1% of the total variation in this dataset are still relevant and point to coherent and important gene sets and their related samples.
Relationships between extreme gene sets and the corresponding sets of driving samples could be discerned for many PCs. The Methodology section presents a way to additionally correlate each principal component's sample partitioning with any available sample covariates. Although some human sample covariate information is provided in our test case, the GNF human expression dataset is not amenable to this additional layer of analysis because multiple subject's RNA samples were pooled prior to amplification and array hybridization. However, a second publicly available dataset with rich covariate information is presented below.

Application to microarray expression data, case 2: human diabetes
We acquired a Human diabetes expression dataset [37] from the Broad Institute Cancer Program dataset repository [38] along with the corresponding phenotype covariate data, and applied the filtering step as they described to produce a set containing 10,983 probes across 43 sam-         [35]).
This dataset contains more than 50 covariates, which provides the opportunity to interpret each principal component by searching for covariates that correlate well with expression patterns in the PCEG sets. As described in the Methodology section, we asked if any of the covariate annotations are well correlated with the partitioning of samples into UP N , FLAT N and DOWN N sets. Covariate distributions were compared across different partitions (when sufficient data was available) and any significant trends identified were recorded (see Table 7). For covariates identified as significantly correlated with a principal component's sample partitioning, covariate distribution plots were generated to further investigate and evaluate the apparent relationship. For example, Figure 5 illustrates that PC14's UP 14 , FLAT 14 and DOWN 14 sample partitions appear to be significantly related to two covariate measurements: Insulin_0 (sig = 0.0010); and Type2b_(%) (sig= 0.0077). The Pearson's correlation between the mean expression for the PC14EG-high set and Insulin_0 and Type2b_(%) covariates are r = 0.411 and r = 0.467 respectively.

Discussion
This PCA-based data-mining tool highlights specific patterns of expression and associates them in a convenient way with the genes and samples responsible for those patterns. Some associations in the first few principal components (PCs) of the GNF set reflect major features in the data that are expected. This includes the global high and low constitutive expression profiles of PC1 (67% of variance in the GNF dataset). A component similar to this is often the first or second PC in Affymetrix array datasets. GNF PC3, in contrast, highlighted brain/neuronal tissues, which we expected in this dataset because there are many more samples from brain regions than from any other tissue, and there are thousands of genes that are expressed in  H 7 )-mean(L 7 )), which has the effect of placing the most significantly affected tissues at the top of each list. The most significant tissues in UP 7 are at the left of Figure 2, and the most significant conditions in DOWN 7 are at the right of Figure 2.

UP7 DOWN7
SmoothMuscle   21 , ordered most extreme first, including common name, description, and associated Gene Ontology terms (provided with the GNF dataset). These extreme genes are shown as blue points in Figure 3.   21 , ordered most extreme first, including common name, description, and associated Gene Ontology terms (provided with the GNF dataset). These extreme genes are shown as blue points in Figure 3. (Continued) a general brain pattern. The GO enrichment analysis associated with PC3lowEGs confirmed this impression by identifying neurogenesis, central nervous system, and synaptic terms as significantly enriched for PC3EG (Table 6).

PC-21 Value
We asked if principal components that account individually for small fractions of variation in the data are likely to be significant. Conventional practice generally ignores principal components accounting for a few percent, or less, of total variation, on the assumption that such minor components are most likely dominated by noise. We believe this assumption, that all of the PC's accounting for small fractions of data should be ignored because they are artefact, to be wrong in the context of our analysis. We believe this, in part, because in our analysis we find many PCEGs for minor components are statistically enriched (Table 6). Further, computational experiments using randomized data fail to produce any significant enrichments. This GO enrichment analysis, it should be noted, tends to underestimate the fraction of gene sets that are significant. This is a known artefact in instances where the input gene number (here those passing the selected rank sum threshold) is small, and its effects are exacerbated by the fact that GO annotations for human are still very much in a building phase. Many genes that will eventually be associated with GO terms are not yet entered. This means that reducing the threshold modestly, and therefore increasing the gene number, can uncover additional significant GO term enrichment in some of these PCs. For complete results of GNF analysis at p < = 0.001 see our supplemental materials web site [35].
Viewed from a biological perspective, this PCA mining revealed several different classes of relationships. GNF PC21 is a good example of a component that highlights a coherent gene set and its corresponding tissues, many of which would also be grouped together by conventional clustering algorithms. This is true even though PC21 accounts for just 0.25% of total variation. The PC21EGlow set (defined at p < = 0.00001) was enriched in a statistically significant way for five different GO terms (Table  6), and these terms (myogenesis, muscle contraction, etc.) tell a simple and internally consistent story about muscle development and function. The top tissue samples driving PC21-low are skeletal muscle, tongue (also composed largely of striated muscle), heart and thyroid (which includes a population of myoid cells). Most top-ranked genes in this PCEG set are so specific for striated muscle that they would also appear together in conventional clusterings, although many clustering approaches become technically problematic with datasets of this size. How-  Table 2.

UP21 DOWN21
Pancreas ever the PC21EG-low list differs from a conventional muscle cluster because it also includes some genes that are partly associated with muscle and partly associated with other tissues. PP1R1A, a regulatory subunit of protein phosphatase1, is such a gene. A role for it in striated muscle is suggested based on its coherent presence in tongue, skeletal and cardiac samples, even though it might well not have been seen in this light by standard clustering.
The second example is GNF PC7, which accounts for 0.81% of variation. It illustrates a different kind of biological relationship that more strongly distinguishes results of PCA mining from classical clustering. Top extreme genes associated with PC7 by inspection turn out to be a "who's who" of extracellular matrix components (a specific subset of fibronectins, collagens, laminins plus matrix associated proteins like MFAP5, MGP, LUM; regulatory molecules that mediate stability and function of those matrix components (thrombospondin, SPARC, ADAMTS1, Plod2); and matrix associated signalling and matrix associated signal modulators (insulin like growth factor binding proteins 7, 8 and 10; Sema3c). GO analysis confirms what inspection of the top PC7EGs suggested: namely that a set of extracellular matrix components are expressed in these driving tissues. It is instructive to look at the individual expression profiles for these genes directly at the GNF website and also in aggregate, as represented in the tissue (conditions) list in Table 2. The most prominent contributing tissues associated with high expression of these genes are informative specifically because they do not constitute a group that would have been selected a priori as a coherent set based on known tissue function or shared developmental origin.
This is useful because a biologist interrogating the GNF database would not likely have constructed a query combining adipocytes, smooth muscle, bronchial epithelium, nor would one expect traditional clustering algorithms to place these genes so close to each other as to catalyze the same observation. Similarly, conventional ordering algorithms would not have placed them adjacent to each other because other parts of their expression profiles, containing different genes than those in PC7, would dominate their positions. And concerning genes, they have, in addition to commonalities of expression highlighted by PC7, differences from each other in additional diverse tissues not highlighted by PC7. The PCA grouping gives impetus and a necessary starting gene list to search for one or more factors or regulatory RNAs with a similar expression pattern, or to search for a shared and perhaps evolutionarily conserved cis-acting DNA sequence motifs. It is unlikely that these working hypotheses would have been arrived at easily by widely used methods of gene expression analysis.  Table 4. , and thus grouped by significance (UP 7 group at left, FLAT 7 group in middle and LOW 7 group at right) at test1Thresh = 0.05. Table 3 lists the tissues within UP 7 and DOWN 7 that occur at each end of this plot.

GNF PC7 Extreme Gene Trajectory Plot with Tissues Ordered by Significance
The diabetes dataset offered us an opportunity to add more value to principal component interpretations by searching for covariates that appear correlated at some significance level. The relationships highlighted between covariates and principal components are suggestive, but not conclusive by themselves. Rather, they provide hypotheses that a researcher may wish to further investigate. While we have not delved deeply into this dataset, we believe that a number of the principal components are highlighting a number of meaningful sources of variation present. It is not clear whether this set should exhibit the same proportion of meaningful principal components as the GNF dataset, as by design the NGT vs. GM2 dataset does not contain substantial diversity of samples. Likewise, the selection of covariates was focused on a narrow set of measurements selected to be indicative of diabetes status, and so many covariates are redundant. We anticipate this tool will be maximally useful in cases where datasets are rich in both sample complexity and diversity of covariates.

Conclusion
Results presented above show that this PCA-mining approach can guide a user to biologically significant observations that both complement and reinforce those from conventional clustering analysis. The software package and web interface make this style of microarray analysis straightforward and accessible.
We have applied this to four additional microarray datasets (as yet unpublished) and to one multi-spectral imaging dataset. In each case we found the interpretations that the tool presented to be useful. In general, it seems that the top few principal components identify very broad characteristics of the data. Digging to the deeper components that comprise smaller but more particular substructure in the data leads to more subtle but often meaningful observations, many being complementary to standard clustering.
With respect to the top few components, PC1 is usually the approximate diagonal through the sample/condition space, explaining the overall variation in absolute expression level. For some other datasets we have noticed that the top few PCs can also highlight effects of preprocessing normalization steps or global data quality issues. This means they do not necessarily expose the most important biological variation. Thus, in one microarray dataset not shown here, PC2 was found to be extremely well correlated with a measure of quality of samples, as reflected by the percent of Affymetrix probes called present. Given this evidence of data quality effects comprising a major source of variation over the entire dataset, one might be motivated to remove the major offending conditions, and then repeat the PCA interpretation on the remaining conditions (columns). The idea is that an independent source of variation might be obscured by more dominant signals or noise present in the data from the offending condition.
Our experience thus far leads us to think that this PCA interpretation method will contribute to microarray expression analysis, as one part of a panel of methods that are sensitive to different features in a dataset, such as sample number, gene number, and distribution of variation across the samples. The PCA method should be especially useful for large, complex datasets that offer rich variation among many samples. What is certain is that there are almost always multiple sources of variation in a dataset and that in any specific study their nature and relative strength is informative, whether the origin is an easilyunderstood biologic one, a technical one, or a poorlyunderstood but nonetheless biologically pertinent one.
We are continuing to explore ways to improve our methodology and software package. We anticipate further advances will come with software infrastructure improvements to permit covariate analyses of both column (sample) covariates and row (gene) covariates. The CompClust dataset labeling capability [27] allows a user to attach diverse and numerous labelings to rows or columns. For example we can pull in additional row (gene probe) annotations such as Gene Ontology (GO) functional groups. Beyond explicitly comparing the NMI significance of specific row partitionings for discrete covariates, we GNF PC21 Extreme Gene Trajectory Plot with Tissues Ordered by Significance  Table 5 lists the tissues within UP 21 and DOWN 21 that occur at each end of this plot.
plan to add routines to CompClust to automatically indicate when a group of genes are found to be enriched in specific GO categories (as was done in our analyses above), and more generally to handle large, discrete, multi-valued distributions of values. Our use of NMI treats discrete covariates as discrete random variables that can have at most a single value per condition, and so does not optimally address issue of multi-valued discrete random variables (e.g. GNF data has covariate "concomitant medications" with values like "aspirin", "tylenol", and "aspirin & tylenol"). We are considering more elaborate extensions of mutual information or alternatives that might be able to take further advantage of such multi-valued entries.

Availability and requirements
The PCA interpretation software is implemented as one component of the CompClust Python package [9,27], which is freely available for non-commercial use. The software capability is also accessible through the CompClust-Diabetes PC14 Sample Partitioning is Correlated with Certain Covariates Figure 5 Diabetes PC14 Sample Partitioning is Correlated with Certain Covariates. When sufficient covariate data is available (a number of the measurements are missing for certain covariates), covariate distributions are compared across partitions and significant differences are reported (as in Table 7). When a covariate is identified as significantly correlated with a principal component's sample partitioning, covariate distribution plots can be generated to further investigate and evaluate the apparent relationship. For example, the diabetes dataset's PC14 extreme genes partition the samples into UP 14

PCEG Set GO Description (p-value)
Web web-based interface [29]. The software that implements the web application is also included within the CompClust software distribution.

Methodology
We have developed the following algorithm for identifying and analyzing multiple independent sources of variance present within multi-dimensional sample datasets, in particular those that are produced by gene microarray expression experiments. The overall approach can be summarized as follows: 1) perform principal components analysis of the dataset; for each principal component: 2) identify the most extreme gene probes (those with the highest or lowest weighting) for that principal component; 3) identify and group any conditions in which those extreme probes vary significantly; 4) identify any condition covariates that correlate well with the condition grouping. By extending the interpretation of each principal component from extreme genes (rows) to ordered groups of significant conditions (columns) and further to identifying statistically significant correlations with column covariates, we attempt to make full use of the available data, in an objective and data-driven way, to analyze and provide meaningful interpretations of the diverse sources of variation present within the dataset. The PCEGs can be identified and ranked in one of two ways: by identifying points having a low probability (p < = extremeThresh) of belonging to a Gaussian fit to the distribution of points along the PC n axis, or by taking a fixed number of nExtreme points at each tail of the distribution. H n is the resulting set of data points having the highest coordinate values for PC n , and L n is the resulting set of data points having the lowest coordinate values for PC n . These high and low extreme gene point sets are informative in and of themselves because they represent the most extreme of the data points along a principal axis of variation. As such, the high and low PCEG sets are some of the primary outputs generated by our procedure. We use the term "extreme" in a very general sense, in that the points stand out because they are far from the main distribution. We do not mean to imply that such points are  Figure 5).

Determine the principal components of the dataset
Principal Component Number 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 either biologically relevant or nuisance data that should be removed; rather we are interested in these points in an unbiased way. By further analyzing their pattern of expression in the original axes we hope to gain a better understanding of their possible biologic significance.

Identify significant conditions for each principal component
The extreme gene points comprising H n are located near one edge of the high-dimensional cloud of points, and points L n are near the opposite edge. Thus, points H n are likely to have coordinate values that are maximally different from points L n in a subset of the original column space coordinate system. Our procedure next seeks to identify in which of the original columns (the original axes or dimensions) we find the greatest difference of values for points H n versus points L n . We do this by comparing distributions of values in H nj versus L nj for each of the original columns j (j ∈ [1, nc]). A two-sided Wilcoxon rank sum test is used to estimate the likelihood that these two sets of values are drawn from the same distribution [39]; the resulting p-values for each column are used to rank order and group the columns, rather than as actual probabilities. Columns having a likelihood less than a user-defined significance level test1Thresh are identified and placed into one of two column sets: UP n for those where column j has mean(H nj ) > mean(L nj ), and DOWN n for those columns j where mean(H nj ) <mean(L nj ). Remaining columns that do not show significant variation are placed in the column set FLAT n . The column sets are also meaningful outputs of our procedure, as UP n and DOWN n describe the groups of columns in which the extreme genes H n and L n vary significantly. Our procedure can output these columns and column sets in various orders simply as an aid to human interpretation, including: original column order; grouped by set and within set ordered by the Wilcoxon p-value significance; ordered by mean difference, mean(H nj )mean(L nj ); or ordered by the eigenvector column loading. Taken together, the PCEG point sets and significant column sets should provide valuable insight to researchers wishing to interpret each of the sources of variation identified by the principal components procedure.

Interpret each principal component using covariate annotations
When provided additional covariate annotations C, the procedure seeks to determine which, if any of the annotations C k are well correlated with the partitioning of columns into the sets {UP n , FLAT n , DOWN n }. A discrete annotation C k containing m unique values V 1 , V 2 , ... V m also defines a partitioning of the columns {KV 1 , KV 2 , ..., KV m } where KV 1 is the set of columns that share the value V 1 , KV 2 are those that share value V 2 , and so on. An information theoretic measure known as normalized mutual information (NMI) [28] describes the degree to which two discrete random variables share information. When there is high mutual information, knowing the value of one of the variables should be useful predictor of the other variable. (See [9] for a description of the merits of NMI in terms of clustering and understanding the relationships between clusterings.) We construct the 3 x m confusion matrix to compare the {UP n , FLAT n , DOWN n } column partitioning with the {KV 1 , KV 2 , ..., KV m } partitioning and calculate an NMI score between the partitionings. Because the usual NMI score is not symmetric (i.e. NMI(r,c) ≠ NMI(c,r)), we use a variant that we refer to as the average NMI score, which is simply the average of the NMI of the confusion matrix and the NMI of the transpose of the confusion matrix. Those covariates C k having an average NMI greater than a user defined threshold nmiThresh are added to the set of significant covariate annotations A n .
We apply a different approach when evaluating the C k that are continuous covariates. We need to assess whether each principal component's {UP n , FLAT n , DOWN n } column partitioning correlates with each C k distribution of values. We can separately score three different partitioning schemes, UP n vs. DOWN n , {UP n +FLAT n } vs. DOWN n , and UP n vs.{FLAT n +DOWN n }, by determining if C k 's value distributions differ significantly across the partition. E.g. does C k within UP n have a different distribution than C k within DOWN n ? Next, does C k within {UP n + FLAT n } have a different distribution than C k within DOWN n ? Finally, does C k within UP n have a different distribution than C k within {FLAT n +DOWN n }? For each partition scheme we again use a two-tailed Wilcoxon rank sum test, including the small sample adjustments when sample size is less than 10, to determine whether the covariate's value distributions on each side of the partitioning differ significantly from each other. A minSetSize parameter can be specified as desired to reduce false positives when set sizes are very small, e.g. when comparing a distribution of 2 values vs. 7 values. Thus we calculate the Wilcoxon p-value for three partitionings of columns: UP n vs. DOWN n , {UP n +FLAT n } vs. DOWN n , and UP n vs. {FLAT n +DOWN n }. Those covariates having a p-value less than the user defined threshold test2Thresh for any of the three partitionings are added to the set of significant covariate annotations A n .
Upon completion of the covariate analysis, covariates in the set A n , that previously met user-controlled significance thresholds are reported by the software. Covariate reports provide the following information: For discrete-valued covariates, corresponding confusion matrices and average NMI scores are reported; for continuous-valued covariates, the three Wilcoxon p-values are reported together with supporting plots illustrating the covariate distributions. There is no guarantee that any covariates will be significantly related to a principal component. Conversely, spurious relationships might be reported, especially in the case of small numbers of samples due to small column partitions. The tool simply points to those covariates related to a principal component that also satisfy a usercontrolled significance threshold. It is up to the investigator to consider these hypotheses and to confirm the interesting ones through further investigation.

Terminating condition
We have shown that some large-scale expression datasets have biologically pertinent structure that is revealed by deep PC analysis that goes well beyond the first few principal components. However there are limits to the depth of mining and these limits depend on both size and character of the dataset. In all cases, the last principal component is not free to seek a source of variation because it must be orthogonal to all prior nc-1 components. To some degree that is also true of some portion of latter principal components that explain ever-diminishing fractions of the variance. We suggest that a natural terminating condition exists: When a principal component cannot find any columns in which the extreme gene sets show significant differences, there is no need to proceed to subsequent principal components. We observe that this condition is often not met because the extreme genes are typically differentially expressed in at least a few of the original columns (the original axes or dimensions), even for the most minor principal components. We also observe that variants of a dataset (e.g. representative column subsets) can affect the relative ordering, but not the existence, of multiple factors or sources of variation that are reflected in the minor principal component regime. We may therefore choose to investigate all principal components, but do so with the expectation that minor principal components will describe increasingly subtle sources of variation, which can, and often do, include noisy processes inherent in the data source.