&CHAPTER 4 Robust Methods for Microarray Analysis GEORGE S. DAVIDSON, SHAWN MARTIN, KEVIN W. BOYACK, BRIAN N. WYLIE, JUANITA MARTINEZ, ANTHONY ARAGON, MARGARET WERNER-WASHBURNE, MÓNICA MOSQUERA-CARO, and CHERYL WILLMAN 4.1. INTRODUCTION The analysis of a complex system within an environment that is only subject to incomplete control is nearly impossible without some way to measure a large fraction of the system’s internal state information. As a result, it is only with the recent advent of high-throughput measurement technologies (able to simultaneously measure tens of thousands of molecular concentrations) that systems biology is really a possibility. As an example of the scope of this problem, consider that eukaryotic cells typically have on the order of tens of thousands of genes, each of which is likely to have several alternative splicing variants coding for the protein building blocks of the cell. These proteins undergo posttranslational modifications and have multiple phosphorylations such that there are likely to be hundreds of thousands, possibly millions, of variants. Hence the future of systems biology relies critically on high-throughput instruments, such as microarrays and dual mass spectrometers. The research reported here addresses the need for improved informatics to deal with the large volume of information from these techniques. At present, microarray data are notoriously noisy. Hopefully this technology will improve in the future, but for the immediate present, it is important that the analysis tools and informatics systems developed for microarray analysis admit and adjust to significant levels of noise. In particular, such methods should be stable in the presence of noise and should assign measures of confidence to their own results. In this chapter we describe our efforts to implement and assess reliable methods for microarray analysis. We begin with the structure of the typical microarray experiment, normalization of the resulting data, and ways to find relationships. We then Genomics and Proteomics Engineering in Medicine and Biology. Edited by Metin Akay Copyright # 2007 the Institute of Electrical and Electronics Engineers, Inc. 99 100 ROBUST METHODS FOR MICROARRAY ANALYSIS proceed to discuss the tools themselves as well as methods to assess the output of such tools. Much of the work in this chapter is based on the VxInsight visualization software [1]. However, we have tried to minimize overlap with existing published results and have emphasized new work. Throughout each section we will present results and examples from our research to motivate the specific approaches, algorithms, and analysis methods we have developed. The data sets we have analyzed will not be emphasized, although they were of course critical to the work. Most of the data sets we have used have been published, and these will be cited in the course of the chapter. The data sets so far unpublished have been provided by the University of New Mexico Cancer Research and Treatment Center. These data sets include an infant leukemia data set, a precursor-B childhood leukemia data set, and an adult leukemia data set, all of which have been presented at the American Society of Hematology conferences in 2002 and/or 2003. Publications concerning the biological implications of these data sets are forthcoming. Here we discuss methods only. Finally, we emphasize that this chapter is by no means a survey of techniques for microarray analysis. Indeed, the literature on microarray analysis is vast, and new papers appear frequently. We will mention, however, some of the seminal papers, including [2, 3], which describe the original technology, as well as [4], which gives an overview of the earlier work on microarray technology. Of course, microarrays would not be terribly useful without computational analysis, and some of the original work on microarray analysis includes hierarchical clustering of the yeast cell cycle [5, 6], an analysis of leukemia microarray data using an original method of gene selection [7], and an application of the singular-value decomposition to microarray data [8]. In addition to innumerable additional papers, there are a variety of books available which emphasize different aspects of the microarrays and microarray analysis. A few of the more recent books include [9– 12]. 4.2. MICROARRAY EXPERIMENTS AND ANALYSIS METHODS 4.2.1. Microarray Experiments Microarray experiments and their analyses seek to detect effects in gene expression under different treatments or natural conditions with the goal of clarifying the cellular mechanisms involved in the cells’ differential responses. Uncertainty is the rule rather than the exception in these analyses. First, the underlying systems (cells and/or tissues) are incredibly complex whether viewed from the dynamic process perspective or their physical realizations in space and time. Second, there is abundant variability between cells experiencing exactly the same conditions as a result of genetic polymorphisms as well as the stochastic nature of these chemical systems. Third, the collection and initial preservation of these cells are not a precisely controlled process. For example, when a surgeon removes a cancer, the tissue may not be frozen or otherwise processed for minutes to hours, all the 4.2. MICROARRAY EXPERIMENTS AND ANALYSIS METHODS 101 while the cells continue to respond to these unnatural conditions. Further, the tissues, or partially processed extracts, are often sent to another laboratory several hours or even days away from the original collection site, all of which offers opportunities for chemical changes. Fourth, these measurements are not easy to make; they involve many processing steps with a wide variety of chemicals, and at every step variability arises. The processing is often (necessarily) divided across several days and among several technicians, with inherently different skills and training. Further, the chemicals are never exactly the same; they are created in different batches and age at different rates. All of these affect the laboratory yields and the quality of the processing. Finally, the arrays themselves are technological objects subject to all sorts of variability in their creation, storage, and final use. In effect, the simple measurement of mRNA concentrations that we would like to make is confounded by huge uncertainties. To be able to make good measurements, it is essential that all of the mentioned steps be subject to careful statistical process control, monitoring, and systematic improvement. Further, the actual experiments should be designed to avoid, randomize, or otherwise balance the confounding effects for the most important experimental measurements [13, 14]. These are best practices. Unfortunately, they are not often followed in the laboratory. By the time the data are ready for analysis, they are typically presented in a numeric table recording a measurement for each gene across several microarrays. For example, one might be analyzing 400 arrays each with 20,000 gene measurements, which would be presented in a table with 20,000 rows and 400 columns. Often the table will have missing values resulting from scratched arrays, poor hybridizations, or scanner misalignments, to name just a few (from among a host) of the possible problems. The analysis methods should be able to gracefully deal with these incomplete data sets, while the analyst should approach these data with great skepticism and humility considering the complexity of the cellular processes and the error-prone nature of our microarray technology. Despite all of these issues, statisticians and informaticians, unlike mathematicians, are expected to say something about the structure and meaning of these data. Because, as Thompson has said, “[statisticians] should be concerned with a reality beyond data-free formalism. That is why statistics is not simply a subset of mathematics” [15, p. xv]. Here, we attempt to follow Thompson in discussing implications of, as well as our approaches to, analyzing these experiments, including considerable detail about the algorithms and the way the data are handled. 4.2.2. Preprocessing and Normalization As discussed in the previous section, microarray data typically have a large number of missing data, or values otherwise deemed to be nonpresent. We typically drop genes with too many missing values, assigning a threshold for this purpose that is under the control of the analyst. The raw values are then scaled to help with the processing. The distributions of microarray measurement values typically have extremely long tails. There are a few genes with very large expressions, while most of the 102 ROBUST METHODS FOR MICROARRAY ANALYSIS others are quite small. Tukey and Mosteller suggested a number of transforms to make data from such distributions less extreme and more like the normal Gaussian distribution [16]. In particular, taking logarithms of the raw data is a common practice to make microarray data more symmetric and to shorten the extreme tail (see Fig. 4.1a). However, we frequently use another transform to compress the extreme values, which is due to [17]. This rank-order-based score is an increasing function of expression level, for which the smaller values are compressed to be very nearly the same. This is particularly useful with array data, as many of these smaller values are due purely to noise. The Savage score is computed as follows: If the absolute expression levels are rank ordered from smallest to largest, X(1)  X(2)      X(n) , then the score for X(k) is given by SS(X(k) ) ¼ k X 1 n þ 1j j¼1 Figure 4.1b shows how this score affects the resulting distribution. In particular, Savage scoring compresses the extreme tail and emphasizes the intermediate and large counts. In this case, about 60% of the savage scores are below 1. The use of this scoring has two advantages over correlations using raw counts. First, because it is based on rank ordering, data from arrays processed with very different scales can still be compared. Second, because the noisiest fraction of the measurements is aggressively forced toward zero, the effect of the noise is suppressed without completely ignoring the information (it has been taken into account during the sorting phase). Finally, large differences in rank order will still be strong enough to be detected. The normalization of array data is controversial, although some form of centering and variance normalization is often applied [18, 19]. However, it has been argued FIGURE 4.1. (a) Distribution of log-transformed data from typical Affymetrix U94A microarray. (b) Distribution of the same data after Savage scores. 4.3. UNSUPERVISED METHODS 103 that for many experiments there is no intrinsic reason to expect the underlying mRNA concentrations to have the same mean (or median) values, and hence variance adjustment is suspect. Nevertheless, the analyst has the option to do such normalizations, if desired. In general, we avoid this issue by working with order statistics and Savage scores. 4.2.3. Overview of Basic Analysis Method After preprocessing the measurements with thresholding, rescaling, and various other transforms, we often perform our analysis as follows. First, we compare the genes or arrays under investigation by computing pairwise similarities with several techniques. These similarities are used to cluster (or assign spatial coordinates to) the genes in ways that bring similar genes closer together. These clusters are then visualized with VxInsight [1, 20]. Although we typically use genes in our analyses, we often use arrays as well. In fact, the analysis is the same, as all our techniques use either the gene matrix initially described, or the transpose of the gene matrix. Hence, throughout the remainder of this chapter, we will use the terms genes and arrays interchangeably. Next, the clusters are tested with statistical methods to identify genes and groups of genes that are differentially expressed in the identified clusters or genes otherwise identified with respect to experimental questions and hypotheses. The expression values for the identified genes are plotted and tested for stability. Those genes which seem particularly diagnostic or differentiating are studied in detail by reading the available information in online databases and in the original literature. Each of these analysis steps will be presented in an order approximately following the analysis order we use in practice. This type of analysis is quite typical for microarrays and is usually divided into two categories. The method first described for clustering genes is known as cluster analysis, or unsupervised learning. The term unsupervised is used because we are attempting to divine groups within the data set while avoiding the use of a priori knowledge. In contrast, the second method described for the identification of gene lists is usually called supervised learning. The term supervised refers to the fact that we are using prior information in our analysis, in this case attempting to discover which genes are differentially expressed between known groups. An unsupervised method asks, in essence: Are there groups in the data set and, if so, what are they? A supervised method asks: Given groups, can we predict group membership and can we learn what is most important in distinguishing the groups? 4.3. UNSUPERVISED METHODS 4.3.1. Overview of Clustering for Microarray Data Organizing large groups of data into clusters is a standard and powerful practice in exploratory data analysis. It has also become a standard practice in microarray 104 ROBUST METHODS FOR MICROARRAY ANALYSIS analysis. Typically, the first step after the initial data transformations involves the pairwise comparison of the data elements to determine their relative similarity or dissimilarity. A single comparison usually results in a single number. Thus when comparing the expressions of n genes across multiple experimental conditions, one might compute n(n 2 1)/2 correlations using the similarity measure between each possible pair of genes. After the data pairs have been assigned similarities, various grouping algorithms can be used. In the case of microarray analysis, there are a variety of methods that have been applied and/or developed. These include hierarchical clustering [5], the singular value decomposition/principal-component analysis [8], and self-organizing maps [21]. Our method (discussed next) belongs to the class of graph-theoretic clustering methods. Other such methods include those presented in [22 –24]. Finally, a general overview of clustering methods can be found in [25]. In this chapter we focus on a tool known as VxOrd [20], which uses a forcedirected graph layout algorithm. Our method requires, as do most clustering algorithms, a single similarity value for each of the data pairs. In the next section we show one way to compute the similarities and then consider the actual clustering algorithm. 4.3.2. Clustering Using VxOrd 4.3.2.1. Choosing a Similarity Measure One obvious candidate for measuring similarities is the simple correlation coefficient R due to Pearson [26], Pd  )(yi  y ) i¼1 (xi  x ffi Rxy ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P Pd  )2 di¼1 (yi  y )2 i¼1 (xi  x where, of course, 1  Rxy  1. Pearson’s R is just a dot product of two meancentered and normalized d-dimensional vectors. Thus, one can think of R as a measure of the extent to which the two vectors point in the same direction in the d-dimensional space. Of course, the vectors might lie along the same line, but point in opposite directions, which is the meaning of Rxy ¼ 1. If the vectors are completely orthogonal, the correlation will be zero. In fact, Pearson’s correlation is the measure of similarity that we, and the rest of the microarray community, use most often. It is, however, not the only possibility and in fact has some features that do not recommend it under certain situations. For instance, Pearson’s correlation is known to be sensitive to outliers, especially when n is small. Technically, R has a breakdown point of 1/n, meaning that as few as one extreme outlier in the n points can make the statistic completely different from the true measure of correlation for two random but correlated variables [26]. In practice, we have not found this to be a real problem, since we typically use hundreds of arrays. However, early in the development of microarray technology, many data sets were published with order tens of arrays. In these cases, it was deemed valuable to apply more computationally expensive but more robust 4.3. UNSUPERVISED METHODS 105 measures of correlation. Details about robust measures of correlation, including the percentage-bend correlation coefficient, can be found in [27]. We have also found occasion to use very different similarity measures. For example, we clustered genes based on the textual similarity of their annotations in the Online Mendelian Inheritance in Man (OMIM) [28] (http://www.ncbi.nlm. nih.gov/omim). In this case, the text was processed with the RetrievalWare search and retrieval package from Convera. RetrievalWare computes the similarity between two text documents with a proprietary algorithm. For each gene annotation the strongest 19 other annotations were accumulated to create a similarity file and then processed to produce a clustering. The more typical processing for microarrays is discussed in the following section. 4.3.2.2. Postprocessing for Similarity Measures While Pearson’s correlation has a breakdown point of 1/n (a single outlier can distort the statistic from one extreme to the other [26]), it is easy to compute and has been widely accepted in the microarray community. Because Savage-scored expression values are bounded, the influence of outliers is less important. As a result, the correlation coefficient is usually the basis of our similarity computations. When too few arrays are available to have confidence in R, the percentage-bend coefficient [27] is to be used instead. It is common to cluster directly with these coefficients. However, doing so ignores much of the available information because R is such a nonlinear function. For example, there is a slight change in significance when comparing two pairs of genes that have R ¼ 0.5 and R ¼ 0.55, respectively, but the relative significance between R ¼ 0.90 and R ¼ 0.95 can be quite large. Consequently, it is better to transform these correlations with a measure of their relative significance. This transform can be done by converting to the t-statistic for the observed correlation R between the pairs of values [26]: pffiffiffiffiffiffiffiffiffiffiffi R d2 t ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  R2 Both R and t have computational issues that should be addressed. In particular, R is undefined when the variance of either x or y vanishes, and hence a minimum, acceptable variance must be determined. We typically require that d X i¼1 (xi  x )2 d X i¼1 ( yi  y )2 . 0:0001 Otherwise, no correlation is computed for the pair of expression vectors. A related issue occurs with t when R approaches +1.0. In this case, the t-statistic becomes arbitrarily large. Because clustering will be comparing similarities, the strength of an extreme outlier will distort the clustering. Hence t should be clipped to avoid such extremes (we typically truncate values greater than 300, though even this value may be extreme). 106 ROBUST METHODS FOR MICROARRAY ANALYSIS Missing data also present concerns in computing R. Certainly, if too many values are missing, any computed similarity would be suspect. Recourse to the analyst’s experience and judgment is the best way to choose how many values can be missing before the comparison is not attempted. For large collections of arrays, requiring that at least 70 measurements be simultaneously present for both expression vectors has been acceptable in our experience. Computing all of these correlations produces a huge file of similarity comparisons. For example, the computation for an experiment [29] around Clostridium elegans, which has about 20,000 genes, required the computation of about 2  108 correlations. Using all of the correlations for clustering is neither necessary nor desirable. Most of the correlations will be modest and including them slows the clustering computation and introduces a great deal of resistance to the process that separates the mass of genes into clusters. If a particular gene has strong correlations with a few tens of other genes, they should eventually be clustered together. However, if there are hundreds or thousands of correlations weakly linking that particular gene to others, then the net sum of these weak correlations may overwhelm the few strong correlations. If only a few of the correlations will be used for clustering, some method of choice is required. The analyst can use all correlations above some threshold or just the strongest few correlations for each gene. We have found the latter approach to be sufficient. We have found that using the 20 strongest correlations is often a good starting point. However, even fewer correlations may suffice, especially with the methods discussed next for finding the most central ordination from a distribution of stochastic clustering results. 4.3.2.3. Clustering by Graph Layout Once the similarities have been computed, the VxOrd algorithm is used to cluster the data. VxOrd considers the data set as an abstract graph. In the case of microarrays, we usually think of the genes as nodes in a graph, and the edges as similarities between genes. The relationship between the data elements and their similarity values can be visualized as an abstract, weighted graph G(Vi , Ei, j , Wi, j ) consisting of a set of vertices V (the genes) and a set of edges E with weights W (the similarities between the genes). This graph is only topologically defined, as the vertices have not been assigned spatial locations. Spatial coordinates, called ordinations, are computed from the weighted graph using VxOrd, which places vertices into clusters on a two-dimensional plane, as shown in Figure 4.2. The ordinations are computed such that the sum of two opposing forces is minimized. One of these forces is repulsive and pushes pairs of vertices away from each other as a function of the density of vertices in the local area. The other force pulls pairs of similar vertices together based on their degree of similarity. The clustering algorithm stops when these forces are in equilibrium. Although the algorithm has been discussed in detail in a previous paper [20], we provide here a brief overview for convenience. VxOrd is based on the approach in [30]. Fruchtermann and Rheingold compute a force term for both attraction and repulsion. These terms are then used to generate new positions for the graph vertices. The algorithm combines the attraction and repulsion terms into one potential energy equation, shown below. The first term, in brackets, is due to the attraction 4.3. UNSUPERVISED METHODS 107 FIGURE 4.2. Data elements are nodes and similarities are edge values, which are clustered and assigned (x, y) coordinates by VxOrd. These coordinates are then used to visualize clusters, as shown on the far right. (From [20].) between connected vertices; the second term is a repulsion term. The equation is given by " # ni X (wi, j  l2i, j ) þ Dx, y Ki(x, y) ¼ j¼1 where Ki(x, y) is the energy of a vertex at a specific (x, y) location, ni is the number of edges connected to vertex i, wi, j is the edge weight between vertex i and the vertex connected by edge j, l 2i, j is the squared distance between vertex i and the vertex at the other end of edge j, and Dx,y is a force term proportional to the density of vertices near (x, y). In our ordinations, the energy equation is minimized iteratively in three phases. The first phase reduces the free energy in the system by expanding vertices toward the general area where they will ultimately belong. The next phase is similar to the quenching step that occurs in simulated annealing algorithms [31], where the nodes take smaller and smaller random jumps to minimize their energy equations. The last phase slowly allows detailed local corrections while avoiding any large, global adjustments. VxOrd also includes additional improvements. These improvements include barrier jumping, which keeps the algorithm from getting trapped in local minima; a grid-based method for computing Dx, y , which reduces the computation of the repulsion term from Q(jVj2) to Q(jVj); and edge cutting, which encourages exposure of clusters in the final stages of the optimization. 4.3.2.4. Clustering Parameters the VxOrd algorithm: The analyst has two important controls over 1. The number of similarities used for the clustering 2. The degree of edge cutting permitted, where edge cutting refers to removing key edges in order to expose cliques in the graph The first control concerns how many similarities are passed to the clustering algorithm. Every gene has some correlation with every other gene; however, most of 108 ROBUST METHODS FOR MICROARRAY ANALYSIS these are not strong correlations and may only reflect random fluctuations. By using only the top few genes most similar to a particular gene, we obtain two benefits: The algorithm runs much faster and, as the number of similar genes is reduced, the average influence of the other, mostly uncorrelated genes diminishes. This change allows the formation of clusters even when the signals are quite weak. However, when too few genes are used in the process, the clusters break up into tiny random islands containing only two or three very similar genes, so selecting this parameter is an iterative process. One trades off confidence in the reliability of the cluster against refinement into subclusters that may suggest biologically important hypotheses. These clusters are only interpreted as suggestions and require further laboratory and literature work before we assign them any biological importance. However, without accepting this trade-off, it may be impossible to uncover any suggestive structure in the collected data. The second control tells the algorithm how many edges may be removed so that cliques in the graph, or clusters in the ordination, may be exposed. This parameter must also be balanced for reliability of the clique against actual exposure of the clique. As an example of the impact of these parameters, consider Figure 4.3. Here we are clustering a set of 126 arrays, each with about 12,000 genes. First consider FIGURE 4.3. Effects of using different numbers of similarity links and different parameters for edge cutting. (a) Using too many similarity links, in this case 30, and only a single undifferentiated group is formed. (b) Using 15 similarity links and the data are no longer completely undifferentiated; some stronger similarities are beginning to force the emergence of structure. (c) Using 30 links but with the maximum edge cutting enabled, and clusters are still not apparent. (d ) Data can be separated the with 15 similarity links and aggressive cutting. 4.3. UNSUPERVISED METHODS 109 the effect of using too many similarities. Figure 4.3a shows the result when 30 similarities per array are used. However, when only the top 15 strongest similarities are used, as in Figure 4.3b, three groups begin to be apparent. When a set of elements has a relatively uniform set of similarities, it can be very difficult to separate them in to subclusters. However, there may be subsets of stronger similarities that would divide the data into clusters if allowed to express their influence in the absence of the other, mostly homogeneous similarities. In other words, we can reveal small cliques of vertices by cutting similarity relationships that have been constraining the vertices to remain an undifferentiated agglomeration. Figure 4.3c shows that no cliques are apparent when using 30 similarities per vertex, even with extremely aggressive edge cutting. On the other hand, the suggestive clusters seen in Figure 4.3b readily break into more detailed cliques when only 15 similarities per vertex are used and when aggressive edge cutting is enabled, as shown in Figure 4.3d. 4.3.2.5. Evaluating Utility and Significance of Clustering Clustering algorithms are designed to find clusters. However, one’s initial stance should be that there is no reason to suppose that the clusters found are more than artifacts. We have expended significant effort devising methods for evaluating the clusters produced by VxOrd. These efforts are described in detail in our previous publications, but for completeness we provide a short overview here. The first, most intuitive approach is to check that gene expressions are correlated within clusters and to investigate the biological meaning of the clusters. One of our first efforts in this direction was an analysis of the Spellman yeast data [32]. In this paper we compared the typical expression histories of the genes in each cluster to assure ourselves that genes in the cluster had, generally, uniform expression patterns and that these patterns were different in the various clusters. Figure 4.4 shows Spellman’s yeast cellcycle data clustered with VxOrd, overlaid with expression traces for typical genes in the various clusters. Not only do these traces seem homogeneous within the clusters and different between clusters, but they also have biological significance as the cells move through their replication cycle. Surprisingly, the various states in the cell cycle correspond to a clockwise progression around the roughly circular layout of gene clusters in this figure. This visual inspection was also recast in the same paper in a statistically more rigorous manner. Although we do not provide the details, we concluded that the VxOrd clusters were not artifacts. A statistical test used in [32] showed that a subset of genes associated with cell cycle phase G1 were collocated with p , 0.001 and further that CLB6, RNR1, CLN2, TOS4, and SVS1 collocated with p , 0.0001 for cells exiting from long-term stationary phase. Another test falling into the category of intuitive verification was performed in [33]. This work tested VxOrd clusters of C. elegans genes for clusters enriched in genes known to be involved in various biological processes. Stuart et al. [33] found significant statistical enrichments. These enrichments suggested that other genes in the same clusters could be expected to be involved in the indicated processes. This hypothesis was confirmed with laboratory experiments. 110 ROBUST METHODS FOR MICROARRAY ANALYSIS FIGURE 4.4. Cell cycle data with typical expression traces from each cluster. Interestingly, the clusters lay out in a circle corresponding to the temporal cell cycle phases. Another evaluation method to investigate the clustering algorithm uses exactly the same processing parameters but with randomly permuted versions of the measurements. If the clustering algorithm finds clusters or structures in this randomized data, then the results with the original data should be suspect. The processing methods discussed above have been tested in this way, and randomized data do not exhibit any organized structure; see, for example, the analysis in [29], where the randomized data resulted in a single, symmetric, and otherwise unorganized group of genes, which revealed structure in the data as well as lack of structure in the randomized data. If randomized data show no structure, then the structures in the actual data become more interesting and may possibly be useful. These methods have been useful in showing that the clusterings are not chance occurrences and have led to scientific insights. However, these approaches have not addressed two other important issues related to clustering: (1) how stable these clusters are with respect to variations in the measurements and (2) how stable they are with respect to different random initializations of the VxOrd clustering algorithm, which has an inherently stochastic nature. We investigated these two issues in [20]. To test the stability of the algorithm with respect to random starting points, we ran 100 reordinations of the Spellman cell cycle data [6], which had about 6000 genes. Each reordination was started with a different seed for the random-number generator. We then visually marked the elements of a cluster in one ordination and looked to see if they were still visually clustered together in the other ordinations. The results of this analysis were generally positive and are shown in 4.3. UNSUPERVISED METHODS 111 FIGURE 4.5. (a) Ordinations with different random starting conditions. (b) Effect of increasing edge noise on cluster stability. (From [20].) Figure 4.5a. We also performed a more rigorous neighborhood analysis (see discussion below) with good results. To determine if small changes or noise in the similarities would give small changes in the ordination results, we ran 80 reordinations where we added noise drawn from a Gaussian distribution with mean zero and standard deviations 0.001, 0.010, 0.050, and 0.100 and recomputed the ordinations. These different ordinations were also compared visually and statistically. These results generally showed that our the clusterings held together remarkably well, even when a large fraction of noise was added. The visual results are shown in Figure 4.5b. 4.3.2.6. Finding Most Representative Clustering Each randomly restarted ordination by VxOrd represents a sample from a distribution of possible ordinations arising from a particular similarity file. From this perspective, one might want to identify the best ordination, which is particularly hard because it is an extreme and further because the concept of best cluster or best ordination is not particularly well defined. However, the idea of a most representative ordination or most central ordination (MCO) can be defined with respect to the sample of observed randomly restarted ordinations. In this case, two ordinations are compared by neighborhood analysis to create a single measure of overall similarity between the two ordinations. With this method for comparing two ordinations, one can make all possible comparisons of the available randomly restarted ordinations and then select the ordination that is, on average, most like all the remaining reordinations. This idea of centrality of the distribution of ordinations might be further extended to the second moment to compute some measure of dispersion, which perhaps could eventually be extended to allow some sort of hypothesis testing about these ordinations. However, we have only investigated the centrality issue. As an example, we used massively parallel computers to calculate hundreds or in some cases thousands of reclustered ordinations with different seeds for the random-number generator. We compared pairs of ordinations by counting, for 112 ROBUST METHODS FOR MICROARRAY ANALYSIS every gene, the number of common neighbors found in each ordination. Typically, we looked in a region containing the 20 nearest neighbors around each gene, in which case one could find a minimum of 0 common neighbors in the two ordinations or a maximum of 20 common neighbors. By summing across every one of the genes, an overall comparison of similarity of the two ordinations was computed. We computed all pairwise comparisons between the randomly restarted ordinations and used those comparisons to find the ordination that had the largest count of similar neighbors. Note that this corresponds to finding the ordination whose comparison with all the others has minimal entropy and in a general sense represents the MCO of the entire set. Although not shown here, plots of the entropies, intersection counts, and cross plots of entropy showed that there were central ordinations, even for a data set that we found very difficult to break into stable clusters [34]. It is possible to use these comparison counts (or entropies) as a derived similarity measure in order to compute another round of ordinations. For example, given that 200 random reordinations have been computed, one can compute the total number of times gene gj turns up in the neighborhood of gene gk in the available 200 ordinations. This count, or the average number of times the two genes are near each other, will be high when the two genes are generally collocated (which should be a reflection of similar expression profiles for gj and gk). The clusters from this recursive use of the ordination algorithm are generally smaller, much tighter, and generally more stable with respect to random starting conditions than any single ordination. In fact, we found that most single ordinations were more similar to the MCO when using the derived similarity than when using Pearson’s R. We typically use all of these methods (computing the MCO from among about 100 to 200 random reordinations and computing neighborhood set sizes ranging from 10 to 30 by steps of 5) during exploratory data analysis in order to develop intuition about the data. As an interesting aside, the process of comparing pairs of ordinations results in similarity values between every ordination. These similarities can be used to create clusters of the clusterings! Figure 4.6 shows an example where we found that the random reclusterings seemed to fall into two different attractor basins, which may be interpreted as a sign that there were two different but perhaps equally valuable ways to look at the data and that no single cluster was able to represent both of these views. 4.3.3. Enhancing VxOrd In addition to developing the MCO algorithm, which is layered on top of the VxOrd algorithm, we put some effort into enhancing VxOrd itself. For clarity, we will hereby denote the enhanced version by VxOrd 2.0 and the original verision by VxOrd 1.5. The original motivation for developing VxOrd 2.0 was to cluster larger data sets, although we also found the algorithm to be more stable and more accurate on certain data sets. 4.3. UNSUPERVISED METHODS 113 FIGURE 4.6. Process of comparing individual clusters results in similarity values between each ordination. These similarities can be used to create a cluster of clusters. In this case, there seem to be two attractors, suggesting the data may have two useful visualizations. 4.3.3.1. Adding Graph Coarsening to VxOrd VxOrd 2.0 is based on graph coarsening, which has previously been applied to other graph layout algorithms in order to draw larger graphs more efficiently [35, 36]. Graph coarsening, also known as multilevel mesh refinement, works by replacing the original graph with a smaller but still representative graph. In the context of graph layout, we draw the smaller graph and use this initial layout as a starting point for drawing the original graph. In fact, we use an iterative approach which typically provides a succession of smaller graphs and hence a succession of graph layouts. Suppose the initial graph is denoted by G0. A coarsened version G1 of G0 is obtained by randomly visiting and subsequently merging nodes as follows: 1. Pick a node vi at random. 2. Choose a neighbor vj of vi such that the edge eij has maximal weight. If no such vj exists or vj has been previously visited/merged, go to step 4. 3. Merge vj into vi by adding edge weights of common neighbors or creating new edges if there are no common neighbors. 4. Repeat until all nodes have been visited/merged. 114 ROBUST METHODS FOR MICROARRAY ANALYSIS This algorithm is fast and has been observed [37] to produce much smaller but still representative versions of the original graph. Once the initial graph G0 has been coarsened, we can draw the smaller version G1 using VxOrd 1.5 and then reverse our coarsening to obtain a layout of G0. In fact, we repeatedly apply the coarsening algorithm to get an even smaller graph. The algorithm that we use is as follows: 1. Apply the coarsening algorithm until we obtain a suitably small graph (say 50 nodes) or until the algorithm fails to make any more progress. We obtain a sequence of graphs G0, G1, . . . , Gm from this process. 2. Draw Gm using VxOrd. 3. Refine Gm to obtain Gm21. Place the additional nodes obtained by this refinement in the same positions as their counterparts (merged nodes in Gm) and adjust with VxOrd. 4. Repeat step 3 using Gm22, . . . , G0. This algorithm requires additional adjustments in terms of the grid-based density calculations and the various parameters involved in the simulated annealing. With proper adjustments, however, we obtain better accuracy and stability with this algorithm (VxOrd 2.0) than with the previous version (VxOrd 1.5), as will be demonstrated in the following section. 4.3.3.2. Benchmarking VxOrd 2.0 We first benchmarked VxOrd 2.0 on the so-called swiss roll data set. Although this is not microarray data, it provides a useful example of the essential difference between VxOrd 2.0 and VxOrd 1.5. This data set was used in [38, 39] to test two different nonlinear dimensionality reduction algorithms. The data set is provided as points in three dimensions, which give a spiral embedding of a two-dimensional ribbon (see Fig. 4.7a). To test VxOrd, we considered each point to be a node in an abstract graph, and FIGURE 4.7. The swiss roll is a two-dimensional manifold embedded nonlinearly in three dimensions: (a) actual data set; layouts of associated graph using (b) VxOrd 1.5 and (c) 2.0 VxOrd. 4.3. UNSUPERVISED METHODS 115 we connected each node to its 20 nearest neighbors. In principle, VxOrd should draw the graph as a two-dimensional ribbon. The swiss roll is a useful benchmark because it is very easy to visualize but hard enough so that it will confound (at a minimum) any linear algorithm (such as principal – component analysis). It also confounded the original VxOrd 1.5, as can be seen in Figure 4.7b. Looking closely at Figure 4.7b, we can see that VxOrd 1.5 did well on a local scale (i.e., the colors are together and in the correct order) but failed on a global scale (the ribbon is twisted). Once graph coarsening was added to VxOrd, however, the global structure was also ascertained correctly, as shown in Figure 4.7c. In our subsequent analysis we found a similar story: VxOrd 2.0 does well on large data sets on a global scale but otherwise does not improve VxOrd 1.5. For our next benchmark, we again used the swiss roll data set, an adult leukemia data set provided by the University of New Mexico Cancer Research and Treatment Center, and the yeast microarray data set [6] used previously. In order to test the stability of our modification to VxOrd, we performed multiple runs of both VxOrd 1.5 and VxOrd 2.0 with different random starting conditions. We then compared the ordinations using a similarity metric obtained as a modification of a metric discussed in [40]. Our similarity metric is a function se (U, V ), where U, V are two VxOrd layouts of the same m-node data set x1, . . . , xn. The metric is computed by first constructing neighborhood incidence matrices NU,e and NV,e, where N*,e is an n  n matrix N*,* ¼ (nij), with nij ¼  1 if kxi  xj k , e 0 otherwise Now se (U, V) ¼ NU,e  NV,e kNU,e kkNV,e k where NU,e  NV,e is the dot product of NU,e and NV,e when both matrices are considered to be vectors of length n 2. Finally, we note that in order to make sure we can form reasonable e neighborhoods of a given node xi, we first scale both U and V to lie in the area [21,1]  [21,1]. To see how we can use this metric to assess the stability of VxOrd 2.0, we first revisit the swiss roll data set. In the top row of Figure 4.8, we show an all-versusall similarity matrix for 10 different random runs of both VxOrd 1.5 and VxOrd 2.0. This figure confirms the results from Figure 4.7 and shows that the metric is valid. In particular, we see that VxOrd 2.0 arrives at a consistent solution (indicated by higher similarities) while VxOrd 1.5 is less consistent. Next we computed the same similarity matrices for the adult leukemia data set and for the yeast time-series data, as shown in the second and third rows of 116 ROBUST METHODS FOR MICROARRAY ANALYSIS FIGURE 4.8. Comparison of the similarity matrices for random runs of VxOrd 1.5 and 2.0 using swiss roll, adult leukemia, and yeast data sets. On the left we show the runs produces by VxOrd 1.5, and on the right we see the runs produced by VxOrd 2.0. From top to bottom we have the swiss roll with 10 runs, the adult leukemia data set with 100 runs (10 runs for 10 different edge cutting parameters—more aggressive up and left), and the yeast data set with 10 runs. The color scale is shown on the right and is the same for all images. If we look at this figure as a whole, we see that the right-hand column has more red than the left-hand column and hence that VxOrd 2.0 (right column) is generally more stable than VxOrd 1.5 (left column). 4.4. SUPERVISED METHODS 117 TABLE 4.1 Average Values (Excluding the Diagonal) of the Similarity Matrices Shown in Figure 17 Swiss Roll AML Yeast 0.43 0.87 0.33 0.80 0.62 0.69 VxOrd 1.5 VxOrd 2.0 Figure 4.8. We arrived at similar results in each case. We also computed the average similarity across the matrices (excluding the diagonal) for the different cases, shown in Table 4.1. In the case of the adult leukemia data set, we also experimented with the edge cutting feature of VxOrd. In particular, we computed 10 random runs for each of the 10 most aggressive edge cuts. We found that even though VxOrd 2.0 was more consistent overall, it still was not consistent with the most aggressive cut. The yeast data set was larger (6147 nodes compared to 170 nodes in the adult leukemia data set) and the results of both VxOrd 1.5 and 2.0 were fairly consistent. In particular this suggests VxOrd 1.5 still does well on larger data sets without an inherently gridlike structure (as in the swiss roll). 4.4. SUPERVISED METHODS Clustering a microarray data set is typically only the first step in an analysis. While we have presented tools and techniques to help assure a reasonable and stable ordination, we have not yet discussed the most important part of the analysis: the steps necessary to determine the actual biological meaning of the proposed clusters. While we do not typically address the actual biology in a given analysis, we often provide the appropriate tools for such analysis. These tools must be both informative and accessible to the average biologist. 4.4.1. Using VxInsight to Analyze Microarray Data As stated earlier, most of our work is built upon a database with a visual interface known as VxInsight [1]. VxInsight was originally developed for text mining but has been extended for the purpose of microarray analysis. Once an ordination has been obtained (using the methods described previously), the ordination is imported into VxInsight, along with any annotation or clinical information. VxInsight uses a terrain metaphor for the data, which helps the analyst find and memorize many large-scale features in the data. The user can navigate through the terrain by zooming in and out, labeling peaks, displaying the underlying graph structure, and making queries into the annotation or clinical data. In short, VxInsight provides an intuitive visual interface that allows the user to quickly investigate and propose any number of hypotheses. Details about and applications of VxInsight can be found in [1, 20, 29, 32, 33]. 118 ROBUST METHODS FOR MICROARRAY ANALYSIS 4.4.1.1. Typical Steps in Analysis Using VxInsight VxInsight is very useful for an initial sanity check of a data set. We will typically cluster the arrays to look for mistakes in the scanning or data processing which might have duplicated an array. A duplication will often be apparent in the experiment because the pair of duplicated arrays will cluster directly on top of each other and will typically be far from the other clusters. We have discovered that many data sets cluster more by the day the samples were processed, or even by the technician processing the samples, than because of biologically relevant factors. Further investigation is needed, for example, if almost 100% of a particular processing set clusters by itself. In one case we found a very stable ordination consisting of two groups. After much confusion we discovered that the groups were divided by experimental batch and that one of the groups consisted of patients whose samples contained only dead or dying cells (perhaps due to bad reagents or problems with the freezing process). When the experiments were redone, the original clusters dissolved into more biologically meaningful clusters. One can often see the effect of confounding experimental conditions using this same method. For example, if a set of arrays is processed by the date they were collected and the date corresponds to separate individual studies, then the processing set (date) will be confounded with the study number. Well-designed studies control such confounding by randomizing the processing order or by carefully balancing the processing order. However, it is always wise to use these exploratory analysis methods to ensure that your main effect has not, somehow, been confounded. A more interesting phase of analysis begins after obviously bad data have been culled and the remaining data have been reclustered. The data may be clustered in either of two ways. In one approach, the genes are clustered in an effort to identify possible functions for unstudied genes. See, for example, [29, 32]. In the other approach, which is often seen in clinical studies, we cluster the arrays (the patients) by their overall expression patterns. These clusters will hopefully correspond to some important differentiating characteristic, say, something in the clinical information. As the analysis proceeds, various hypotheses are created and tested. VxInsight has plotting features that are helpful here, including a browser page with various plots as well as links to external, Web-based information. Although useful information can be gleaned by simply labeling different peaks in VxInsight, a more systematic method is even more informative. At the highest level, one may wish to select two clusters of arrays and ask: Which genes have significantly differential expressions between these two clusters? Given any method for identifying such genes, it is useful to display them within the context of the cluster-bygenes map. Sometimes the most strongly differentiating genes for the clusters of arrays may not have been previously studied. In this case, it can be very useful to find known genes that cluster around these unstudied genes using the cluster-bygenes map. This process is illustrated in Figure 4.9, which shows the original table of array data, clustered both by arrays and by genes. The lower map represents the result after clustering by arrays and shows two highlighted clusters (colored white and green, respectively). The genes with strongly differential expressions between the groups 4.4. SUPERVISED METHODS 119 FIGURE 4.9. Array of expression data for large number of experiments shown clustered by genes and by arrays. A list of genes with different expressions between two groups of arrays is shown. This list includes a short annotation and links to more extensive, Web-based annotations. of arrays are shown to the right of this map. Note that the list is sorted by a statistical score and also contains links to the available Web-based annotations. A curved arrow in the figure suggests the path between the gene list and the cluster-bygenes image. That connection is implemented with sockets and forms the basis of a more general analysis tool, which allows an arbitrary gene list to be sent from the analysis of the arrays to the analysis of the genes. 4.4.1.2. Generating Gene Lists There are many methods for generating gene lists or finding genes which are expressed differently in different groups. As stated in the introduction, this process is known as supervised learning, since we are using known groups to learn about (and make predictions about) our data set. Finding gene lists in particular is known as feature or variable selection, where the features in this case are genes. There are a wide variety of methods for feature selection, and we do not provide here an extensive survey of this area. We do, however, mention some of the methods developed specifically to select genes for microarray analysis. The method in [7] was one of the first gene selection methods proposed, the method in [41] applied feature selection for support vector machines to microarray data, and [42] discusses a variety of feature selection methods applied to microarray data. 120 ROBUST METHODS FOR MICROARRAY ANALYSIS For our purposes, we use a simple statistical method for gene selection. A geneby-gene comparison between two groups (1 and 2) can be accomplished with a simple t-test. However, we wanted to eventually support comparisons between more than two groups at a time, so we actually used analysis of variance (ANOVA). This processing results in an F-statistic for each gene. The list of genes is sorted to have decreasing F-scores, and then the top 0.01% of the entire list are reported in a Web page format, with links to the associated OMIM pages. The OMIM pages are then examined manually to hypothesize biological differences between the clusters. 4.4.1.3. Gene List Stability An analysis using the gene list feature of VxInsight typically progresses as follows. First, a question is posed within the VxInsight framework and a statistical contrast is computed for that question. The gene list is initially examined to see if any genes are recognized by their short descriptions, which, if available, are included with the genes. The plots are examined, and the OMIM annotations are read. If the gene appears to be important, the literature links and other relevant National Center for Biotechnology Information (NCBI) resources are studied. This analysis step is very labor and knowledge intensive; it requires the bulk of the time needed to make an analysis. As such, it is very important to not waste time following leads that are only weakly indicated. That is to say, before one invests a great deal of time studying the top genes on a list, it is important to know that those highly ranked genes would likely remain highly ranked if the experiment could be repeated or if slight variations or perturbations of the data had occurred. The critical issue about any ordered list of genes is whether we can have any confidence that this list reflects a nonrandom trend. To be very concrete, suppose that My Favorite Gene (MFG) is at the top of the list in our ANOVA calculations, that is, MFG had the largest observed F-statistic from the ANOVA. What can we conclude about the observed ranking for MFG? Certainly, a naive use of the F-statistic has no support because we tested, say, 10,000 genes and found the very largest statistic from all of those tests. So, an F-value for p ¼ 0.001 would likely be exceeded about 10 times in our process, even if all the numbers were random. Hence, the reported F-statistic should only be considered to be an index for ordering the values. However, if we could repeat the experiment and if MFG was truly important, it should, on average, sort into order somewhere near the top of the gene list. We cannot actually repeat the experiment, but we can treat the values collected for a gene as a representative empirical distribution. If we accept that this distribution is representative, then we can draw a new set of values for each of the two groups by resampling the corresponding empirical distributions repeatedly (with replacement), as shown in Figure 4.10. This process is due to Efron and is known as bootstrapping [43]. Now consider Figure 4.11, where we resample for every gene across all of the arrays in the two groups to create, say, 100 new experiments. These experiments are then processed exactly the same way as the original measurements were processed. We compute ANOVA for each gene and then sort the genes by their F-value. As we construct these bootstrapped experiments, we accumulate the 4.4. SUPERVISED METHODS 121 FIGURE 4.10. A bootstrap method uses the actual measured data as estimates for the underlying distribution from which the data were drawn. One can then sample from that estimated underlying distribution by resampling (with replacement) from the actual measurements. FIGURE 4.11. Actual data are processed to create the gene list, shown at the bottom left. The actual data are then resampled to create several bootstrapped data sets. These data sets are processed exactly the same way as the real data to produce a set of gene lists. The average order and the confidence bands for that order can be estimated from this ensemble of bootstrapped gene lists. 122 ROBUST METHODS FOR MICROARRAY ANALYSIS distribution of the location in the list where each gene is likely to appear. Using these bootstrap results one can determine, for each gene, its average order in the gene lists. Although the distributions for such order statistics are known, they are complex. On the other hand, the bootstrapped distributions are easily accumulated, and they are acceptable for our needs. In addition to the average ranking, we count the 95% confidence bands for each gene’s ranking as estimated by the bootstraps. We report both the upper 95% confidence band and the centered 95% confidence interval for each of the genes. The lower limit of this upper 95% confidence band (LLUCB) is recorded for later use (note that 5% of the time we would observe a ranking below LLUCB by random chance, even when our hypothesis is false, given the two empirical distributions). Next, we investigate the p-values for the observed rankings of these genes under the null hypothesis, H0, that there is no difference in gene expression between the two groups (1 and 2). In this case (when H0 is in fact true), the best empirical distribution would be the unordered combination of all the values without respect to their group labels. To test this hypothesis, we create, for example, 10,000 synthetic distributions by bootstrapping from this combined empirical distribution and process them exactly as we did the original data. We are interested in what fraction of the time we observed a particular gene ranking higher in the bootstrapped results than the appropriate critical value. There are several reasonable choices for this critical value. We could use the actual observed ranking or the average ranking from the bootstraps under the assumption that H0 was false. Instead, we take an even more conservative stance and choose a critical value using a power analysis to control our chance of a type II error. We set b ¼ 0.05, or 5%. If H0 were false (i.e., if the groups do have different means), then the earlier bootstrapping experiments suggest that one might randomly observe a ranking as low as LLUCB about 5% of the time. Hence, we examine the later bootstrap experiments (under H0 assumed true and thus no group differences) and find the fraction of the times that we observe a ranking at or above LLUCB. This value is reported, gene by gene, as the p-value for the actual rankings. In essence, we are saying that if H0 is true, then by random chance we would have seen the gene ranking above LLUCB with probability p. As LLUCB is much lower than the actual ranking, this p-value is very conservative for the actual ranking. To investigate the meaning of the actual F-statistics used to index these gene lists, we computed another bootstrap experiment. We were interested in the effect of scaling the original expression values by their Savage-scored order statistics. As previously discussed, this scoring is felt to be more robust than taking logarithms. However, we were concerned that this might influence our p-values, so we developed a code to estimate the expected F-statistic for the mth ranked gene in a gene list from two groups (1 and 2) respectively having j and k arrays. This code computes a large bootstrap after randomizing the Savage scores within each of the j þ k arrays. The code then computes the ANOVA for each gene and eventually sorts the resulting genes into decreasing order by F-statistics. The final result is a p-value (by bootstrap) for the two groups with the specific number of arrays. This 4.4. SUPERVISED METHODS 123 computation is rather intensive and should either be fully tabulated or run only as needed for genes uncovered by the earlier methods. We have not run extensive simulations of this code against the p-values or the list order distributions, but the limited checks did suggest that genes which always ranked near the top of the differentiating gene lists do have rare F-statistics based on the Savage-scored orders relative to the expected random distributions (data not shown). 4.4.1.4. Comparing Gene Lists As mentioned previously, the ANOVA plus bootstrap approach described above is only one way to find genes which may have important roles with respect to particular biological questions. Our collaborators, for example, have used support vector machine recursive feature elimination (SVM RFE) [41], a Bayesian network approach using a feature selection method known as TNoM [44], and a technique based on fuzzy-set theory as well as more classical techniques, such as discriminant analysis. By using several of these methods, one might hope to find a consensus list of genes. Our experience has shown that this is possible. While the lists from different methods are usually not exactly the same, they often have large intersections. However, the simultaneous comparison of multiple lists has been a difficult problem. We have developed a number of methods which have helped us understand that the lists may be different in the details but still very similar biologically. This makes sense considering that different methods might identify different but closely related elements of regulation or interaction networks. In that case, the methods suggest the importance of the network and the particular region in that network, even though they do not identify exactly the same elements. This relatedness suggests something similar to the kind of “guilty-by-association” method that has been used to impute gene functions for unstudied genes that cluster near genes with known function, as in [29]. Indeed, something similar can be used to evaluate the similarity of multiple gene lists. Figure 4.12a shows a VxInsight screen for clusters of genes. Highlighted across the clusters are genes identified by different methods (shown in different colors). In this particular case, one can see that the various methods do identify genes that are generally collocated, which suggests that gene regulations and interacting networks probably do play a strong role with respect to the question under consideration. Here, for example, the question was, “which genes are differentially expressed in two types of cancers [acute lymphoblastic/myeloid leukemia (ALL/AML)]?”. However, multiple methods do not always produce such strong agreement, as shown in Figure 4.12b. In this case the question was, “which genes are predictive for patients who will ultimately have successful treatment outcomes (remission/ failure)?” Unfortunately, this question had no clear answer. Interestingly, the ANOVA-plus-bootstrap method suggests a very stable set of genes for the first question, while the list for the second question is not stable and has confidence bands spanning hundreds of rank-order positions (data not shown). Finally, we have discovered that when two methods produce similar gene lists, the coherence may be due to the underlying similarity of the methods more than to any true biological significance. We discovered this fact using a visualization of the gene lists using principal-component analysis (PCA), a common technique used for 124 ROBUST METHODS FOR MICROARRAY ANALYSIS FIGURE 4.12. (a) General collocation of genes identified by different algorithms (shown with different colors). This collocations suggests that the different methods are in reasonable agreement. (b) Genes selected by each method are widely separated and show no coherence, suggesting that there is a lack of consensus among the methods. dimensionality reduction that involves projections of the original data onto the principal components. These components are ordered according to the amount of data captured by each component. The first component is the most informative, the second is the next most informative, and so on. Further information on PCA can be found in [45, 46]. An example of PCA applied to microarray data can be found in [8]. In our PCA-based visualization of multiple gene lists, each gene is considered to be a point in patient space, where each dimension corresponds to a different patient. Since, in this case, there were 12,000 genes and 126 patients, the spatial representation had 12,000 points (samples) in a 126-dimensional space. Of the 12,000 genes we only considered about 600 that occurred in the different gene lists, reducing our problem to 600 genes in 126 dimensions. Furthermore, because we were mainly interested in how the genes compared as discriminators, and not how their actual expression levels compared, we projected the genes onto the 126-dimensional unit sphere in patient space, as suggested in Figure 4.13a. Geometrically, this corresponds to comparing the directions of the genes in the various gene lists as opposed to their magnitudes. In order to understand this visualization, is it useful to imagine a sphere with a plane passing through the origin. The sphere corresponds to the unit sphere (the sphere with radius 1 centered at the origin) in the patient space and the plane corresponds to the plane determined by the first two principal components. The first principal component points in the radial direction of the sphere and the second principal component is tangential to the sphere at the sphere’s intersection with the first principal component. The vector representing a particular gene and it will intersect the unit sphere, and it will be near the equator of the sphere (unit circle) if it lies in the plane of the first two principal components. To the extent that the gene lies above or below the plane of the first two principal components, the projection of the intersection back down onto the plane will lie further inside the equator. The distribution of these projections onto the principal-component plane suggests how a given method of gene selection identifies important genes. 4.4. SUPERVISED METHODS 125 FIGURE 4.13. (a) A few genes from three different methods are shown intersecting the unit sphere, along with the projections of those intersections down onto the plane of the first two principal components. Note that genes near that plane will have projections that fall close to the arc of the sphere, while those above or below the plane will have intersections that fall well within the equator of the sphere. (b) The ALL-vs-AML gene list comparison. The gene lists that characterize ALL vs. AML are shown, with a different color for each of the methods used to obtain them. In distinguishing ALL from AML we found that most of the genes in the list were colocalized in our representative visualization. (c) Gene lists that characterize remission vs. failure are shown, with a different color for each of the methods used to obtain them. It can be seen in this figure that distinguishing remission from failure is a difficult task. 126 ROBUST METHODS FOR MICROARRAY ANALYSIS For instance, discriminant analysis and the ANOVA methods are much more similar to each other than to the Bayesian network approach. If we use PCA (see [46] for an introduction to PCA), we see further that many methods will be heavily influenced by differences in the first few principal components of the gene expression data. On the other hand, methods such as SVM RFE [41] are able to examine the simultaneous efficacy of groups of genes, some of which, individually, may not be discriminatory in the first or second principal component. One way to understand these differences is by considering where selected genes project onto the plane of the first two principal components: see Figure 4.13a, which schematically represents a few genes from three methods, identified by different colors. It is evident from Figure 4.13b that the gene lists selected for the ALL/AML problem are related. Unfortunately, it is equally obvious that the gene lists selected for the remission/failure problem are unrelated, as shown using the same analysis in Figure 4.13c. When distinguishing ALL from AML, we found that most of the lists were colocalized in our representative visualization (see Figs. 4.12a and 4.13b). When distinguishing remission from failure, on the other hand, we could not arrive at a satisfactory conclusion (see, Figs. 4.12b and 4.13c), which is also consistent with the results from ANOVA plus bootstrapping (data not shown). 4.4.2. Unifying Gene Lists Although it is useful to compare gene lists, the task of sifting through five or six such lists can be very time consuming for the biologist. For this reason, we also developed a quick-and-easy way to combine multiple lists into a single master list. In order to collate and compare the different gene lists we used a weighted voting scheme. In this scheme, we consider genes to be candidates and gene lists to be votes. In other words, each method suggests, in order of preference, which genes should be elected. Our method for combining the gene lists ranks the candidate genes according to the geometric mean of the voting order in each list, where TABLE 4.2 Rank 1 2 3 4 5 6 7 8 9 10 A Simple Combination of Many Gene Lists Geo Mean SVM Stepwise ROC ANOVA TNoM 1 6.17 6.49 8.57 10.14 10.24 10.9 13.34 13.66 13.91 1 4 3 2 5 6 10 8 11 12 1 2 1 6 2 8 3 5 4 11 15 14 1 6 2 3 10 11 4 5 1 3 4 23 11 25 Note: The overall rank was obtained by using the geometric mean of the ranks provided in columns 3 –7. An empty cell in the array indicates a value of 31. ACKNOWLEDGMENTS 127 each method is allowed only 30 votes (the length of our shortest list) and all other genes are given a vote of 31. An example is shown in Table 4.2. 4.5. CONCLUSION At this point in the analysis it may seem that the biology has dissolved into a sea of numbers and statistical methods. However, these methods are our only guideposts when we begin reading the known information about the indicated genes. Without them we could easily waste very valuable time and people in the study of genes which are only weakly, if at all, related to the central questions of the research. Guided by these methods, we can approach the literature with greater confidence and are much more likely to see the important biology reemerge in the gene annotations and the cited literature. However, even after these statistical filters, this literature is vast and is not organized to make our searching particularly easy. We have come to recognize that this step (where very knowledgeable scientists must read extensively) is the critical, rate-limiting step for our research. As a result, we (and many others) have begun work with the natural language processing (NLP) community to build tools that find, summarize, and reorder important parts of the available online literature to make that reading process simpler and more focused toward our research needs. Although we do not discuss such work here, a demonstration of a preliminary automatic Gene List Exploration Environment (GLEE) can be found at http://aiaia.nmsu.edu/. See also [34] Regardless, gene expression studies are providing new insights into molecular mechanism and hold the promise of deeper biological understanding. However, the speed at which groups of genes generated by microarray analysis can be put together in pathways is one of the limiting steps in the translation of these discoveries to applications. Mistakes and dead ends due to faulty microarray analysis tools are a particularly frustrating way to slow this analysis. The methods presented here are potentially useful in uncovering groups of genes that serve to fingerprint biologically important subtypes; further aiding biological discoveries; and refining diagnosis and improving assessment of prognosis. To provide greater confidence in our tools, we have also benchmarked our methods extensively for reliability. In fact, we believe that both of these factors (usefulness and reliability) are equally important, particularly for the analysis of microarray data. ACKNOWLEDGMENTS The authors of this chapter gratefully acknowledge the immense contributions of our collaborators, without which this work could not have been accomplished. Here we reported on the methods and techniques developed at Sandia National Laboratories. However, much of the work discussed sprang from our collaborators; we thank them for their generous time and patience in teaching us enough biology to be helpful and for including us in their laboratories. While developing these methods we have worked directly with and have had fruitful discussions with dozens of researchers from these laboratories and have benefited 128 ROBUST METHODS FOR MICROARRAY ANALYSIS by many generous introductions to other researchers. We owe each of these people our thanks for their critiques of our work and encouragements to continue. Our collaborations with the University of New Mexico have been very fruitful. The informatics work reported here stems directly from close interactions with Vickie Peck, who opened the world of genomics to us. Certainly, without the continuing encouragements and contributions of Maggie Werner-Washburne we would not have created useful tools or have been able to explain them to life scientists. One of the most important of these was Stuart Kim from Stanford University, who was able to see through the primitive nature of our early efforts and recognize what they could become. Our collaboration with Stuart stretched and improved all of our tools and led to an important, joint publication that continues to have impact. The microarray work with Cheryl Willman’s laboratory at the University of New Mexico Cancer Center drove the development of many of the statistical techniques presented here. Cheryl included us in her weekly laboratory meetings, which must surely have been made more tedious by our presence and the continuing need to teach us the rudiments of leukemia biology. Each of these groups is large and we have learned something from every one of you, for which we say thank you. We would especially like to acknowledge the help and collaborations of Moni Kiraly, Jim Lund, Kyle Duke, Min Jiang, Joshua M. Stuart, and Andreas Eizinger from the Kim Laboratory and, of course, Edwina Fuge, Jose Weber, Juanita Martinez, Anthony Aragon, and Angela Rodriguez from the Werner-Washburne Laboratory. From the Willman Laboratory, we must certainly acknowledge Susan Atlas, Paul Helman, Robert Veroff, Erik Andries, Kerem Ar, Yuexian Xu, Huining Kang, Xuefei Wang, Fred Schultz, Maurice Murphy, and particularly Mónica Mosquera-Caro and Jeffrey Potter, who have been immensely helpful to our research. We would particularly like to thank Jon C. Helton for suggesting the use of Savage scoring as a means to normalize microarray data. If we have unfortunately omitted someone to whom we owe our thanks and gratitude, please accept our apologies and recognize that we do value your help. As always, any misrepresentation or error with respect to our collaborators’ work is purely our own fault. Finally, we would like to thank the W. M. Keck foundation for funding the W. M. Keck Genomics Center at the University of New Mexico, which was particularly important in our work. This work was supported by grants from the National Science Foundation, (NSF) (MCB-0092374) to M.W.W., an NSF Minority Post-doctoral fellowship to M.J.M, and U.S. Department of Agriculture (99-38422-8034) to A. D. A. Development of the algorithms described here was funded in part by the W. M. Keck foundation and a Laboratory Directed Research and Development program from Sandia National Laboratories. In addition, a portion of this work was funded by the U.S. Department of Energy’s Genomics: GTL program (www.doegenomestolife.org) under project “Carbon Sequestration in Synechococcus Sp.: From Molecular Machines to Hierarchical Modeling” (www.genomes-to-life.org). Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the U.S. Department of Energy under Contract DE-AC04-94AL85000. REFERENCES 1. G. S. Davidson et al., “Knowledge mining with VxInsight: Discovery through interaction,” J. Intell. Inform. Sys., 11(3): 259 – 285, 1998. 2. M. Schena et al., “Quantitative monitoring of gene expression patterns with a complementary DNA microarray,” Science, 270(5235): 467 – 470, 1995. 3. J. L. DeRisi, V. R. Iyer, and P. O. Brown, “Exploring the metabolic and genetic conrol of gene expression on a genomic scale,” Science, 278(25): 14863– 14868, 1997. REFERENCES 129 4. P. O. Brown and D. Botstein, “Exploring the new world of the genome with DNA microarrays,” Nature Genet., 21: 33– 37, 1999. 5. M. Eisen et al., “Cluster analysis and display of genome-wide expression patterns,” Proc. Natl. Acad. Sci., 95(25): 14863– 14868, 1998. 6. P. Spellman et al., “Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization,” Mol. Biol. Cell, 9(12): 3273– 3297, 1998. 7. T. Golub et al., “Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring,” Science, 286(5439): 531 – 537, 1999. 8. O. Alter, P. O. Brown, and D. Botstein, “Singular value decompositon for genome-wide expression data processing and modeling,” Proc. Natl. Acad. Sci., 97(18): 10101– 10106, 2000. 9. H. C. Causton, J. Quackenbush, and A. Brazma, Microarray Gene Expression Data Analysis: A Beginner’s Guide, Blackwell Publishers, Malden, MA, 2003. 10. S. Draghici, Data Analysis Tools for DNA Microarrays, Chapman & Hall, Boca Raton, FL, 2003. 11. M.-L. T. Lee, Analysis of Microarray Gene Expression Data., Kluwer Academic, Boston, MA, 2004. 12. M. Schena, Microarray Analysis, Wiley, New York, 2002. 13. G. E. P. Box, J. S. Hunter, and W. G. Hunter, Statistics for Experimenters, Design, Innovation, and Discovery, Wiley, Hoboken, NJ, 2005. 14. J. F. Zolman, Biostatistics, Oxford University Press, New York, 1993. 15. J. R. Thompson, Simulation: A Modeler’s approach, Wiley, New York, 2000. 16. J. W. Tukey and F. Mosteller, Data analysis and regression, in F. Mosteller (Ed.), Addison-Wesley Series in Behavioral Science: Quantitative Methods, Addison-Wesley, Reading, MA, 1977. 17. I. R. Savage, “Contributions to the theory of rank order statistics—the two-sample case,” Ann. Math. Stat., 27: 590 – 615, 1956. 18. M. Bilban et al., “Normalizing DNA microarray data,” Curr. Issues Mol. Biol., 4(2): 57– 64, 2002. 19. J. Quackenbush, “Microarray data normalization and transformation,” Nature Genet, 32: 496 – 501, 2002. 20. G. Davidson, B. Wylie, and K. Boyack. “Cluster stability and the use of noise in interpretation of clustering”, in 7th IEEE Symposium on Information Visualization (InfoVis 2001), San Diego, CA, 2001. 21. P. Tamayo et al., “Interpreting patterns of gene expression with self-organizing maps: Methods and application to hematopoietic differentiation,” Proc. Natl. Acad. Sci., IEEE Computer Society, Washington, D.C., 96(6): 2907– 2912, 1999. 22. A. J. Enright and C. A. Ouzounis, “BioLayout—an automatic graph layout algorithm for similarity visualization,” Bioinformatics, 17(9): 853 – 854, 2001. 23. R. Shamir and R. Sharan. “CLICK: A clustering algorithm with applications to gene expression analysis,” in Proc. Int. Conf. Intell. Syst. Mol. Biol. (ISMB), 8:307– 316, AAAI Press, Menlo Park, CA, 2000. 24. Y. Xu, V. Olman, and D. Xu, “Clustering gene expression data using a graph-theoretic approach: an application of minimum spanning trees,” Bioinformatics, 18(2): 536–545, 2002. 130 ROBUST METHODS FOR MICROARRAY ANALYSIS 25. A. K. Jain, M. N. Murty, and P. J. Flynn, “Data clustering: A review,” ACM Comput. Surv., 31(3): 264 – 323, 1999. 26. R. R. Wilcox, Fundamentals of Modern Statistical Methods: Substantially Improving Power and Accuracy, Springer-Verlag, New York, 2001, p. 110. 27. R. R. Wilcox, “Introduction to robust estimation and hypothesis testing”, in G. J. L. a.I. Olkin (Ed.), Statistical Modeling and Decision Science, Academic, San Diego, CA, 1977, p. 188. 28. Online Mendelian Inheritance in Man (OMIM), McKusick-Nathans Institue for Genetic Medicine, John Hopkins University (Baltimore, MD), National Center for Biotechnology Information, and National Library of Medicine (Bethesda, MD), 2000. 29. S. Kim et al., “A gene expression map for Caenorhabditis elegans,” Science, 293(5537): 2087– 2092, 2001. 30. T. Fruchtermann and E. Rheingold, Graph Drawing by Force-Directed Placement, University of Illinois, Urbana-Champagne, IL, 1990. 31. S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, 220(4598): 671 – 680, 1983. 32. M. Werner-Washburne et al., “Comparative analysis of multiple genome-scale data sets,” Genome Res., 12(10): 1564– 1573, 2002. 33. J. M. Stuart et al., “A gene co-expression network for global discovery of conserved genetic modules,” Science, 302: 249 – 255, 2003. 34. G. S. Davidson et al., High Throughput Instruments, Methods, and Informatics for Systems Biology, Sandia National Laboratories, Albuquerque, NM, 2003. 35. Y. Koren, L. Carmel, and D. Harel, “Drawing huge graphs by algebraic multigrid optimization,” Multiscale Modeling and Simulation, 1(4): 645 – 673, 2003. 36. C. Walshaw, “A multilevel algorithm for force-directed graph drawing”, in Proceedings of the Eighth International Symposium on Graph Drawing (GD), Springer-Verlag, London, UK, 2000. 37. B. Hendrickson and R. Leland, “A Multi-Level Algorithm for Partitioning graphs”, in Supercomputing, ACM Press, New York, 1995. 38. S. Roweis and L. Saul, “Nonlinear dimensionality reduction by locally linear embedding,” Science, 290(5500): 2323– 2326, 2000. 39. J. B. Tenenbaum, V. D. Silva, and C. Langford, “A global geometric framework for nonlinear dimensionality reduction,” Science, 290(5500): 2319– 2323, 2000. 40. A. Ben-Hur, A. Elisseeff, and I. Guyon. “A stability based method for discovering structure in clustered data,” in R. Attman, A. K. Dunker, L. Hunter, T. E. Klein (Eds.), Pacific Symposium on Biocomputing, World Scientific Hackensack, NJ, 2002. 41. I. Guyon et al., “Gene selection for cancer classification using support vector machines,” Machine Learning, 46: 389 – 422, 2002. 42. J. Jaeger, R. Sengupta, and W. L. Ruzzo. “Improved gene selection for classification of microarrays,” in R. Attman, A. K. Dunker, L. Hunter, T. E. Klein (Eds.), Pacific Symposium on Biocomputing, World Scientific Hackensack, NJ, 2003. 43. B. Efron, “Bootstrap methods: Another look at the jackknife,” Ann. Stat., 7(37): 1–26, 1979. 44. P. Helman, R. Veroff, S. R. Atlas, C. Willman, “A Bayesian network classification methodology for gene expression data,” J. Comput. Biol., 11(4): 581 – 615, 2004. 45. I. T. Jolliffe, Principal Component Analysis, Springer-Verlag, New York, 2002. 46. M. Kirby, Geometric Data Analysis, Wiley, New York, 2001.