## An empirical study on Principal Component Analysis f

or clustering gene expression data Ka Y ee Y eung W alter L. Ruzzo T echnical Report UW

-CSE-2000-11-03 No v ember, 2000

Department of Computer Science & Engineering

University of W

ashington Seattle, W

A 98195

# An empirical study on Principal Component Analysis for clustering gene expression data Ka Y ee Y eung, W alter L. Ruzzo Dept of Computer Science and Engineering, Univ ersity of W ashington kayee, ruzzo

### cs.w ashington.edu No v 1, 2000
Abstract There is a great need to dev elop analytical methodology to analyze and to exploit the information contained in gene expression data. Because of the large number of genes and the complexity of biological netw orks, clustering is a useful exploratory technique for analysis of gene expression data. Other classical techniques, such as principal component analysis (PCA), ha v e also been applied to analyze gene expression data. Using dif ferent data analysis techniques and dif ferent clustering algorithms to analyze the same data set can lead to v ery dif ferent conclusions. Our goal is to study the ef

fectiv eness of principal components (PC’ s) in capturing cluster structure. In other w ords, we empirically compared the quality of clusters obtained from the original data set to the quality of clusters obtained from clustering the PC’ s using both real gene expression data sets and synthetic data sets. Our empirical study showed that clustering with the PC’ s instead of the original v ariables does not necessarily impro v e cluster quality
. In particular, the first few PC’ s (which contain most of the v ariation in the data) do not necessarily capture most of the cluster structure. W e also showed that clustering with PC’ s has dif ferent impact on dif

ferent algorithms and dif

ferent similarity metrics.

11



## 1 Introduction DNA microarrays of

fer the first great hope to study variations of many genes simultaneously [Lander, 1999 ]. Large amounts of gene expression data ha ve been generated by researchers. There is a great need to develop analytical methodology to analyze and to exploit the information contained in gene expression data [Lander, 1999 ]. Clustering analysis attempts to divide objects into groups such that objects within the same group are more similar to each other than objects in other groups. Because of the large number of genes and the complexity of biological netw orks, clustering is a useful exploratory technique for analysis of gene expression data. Since genes with related functions tend to ha ve similar expression patterns, possible roles for genes with unknown functions can be suggested based on the known functions of some other genes that are placed in the same cluster [Chu et al. , 1998 ]. Clustering is sometimes used as a preprocessing step in inferring regulatory netw orks. For example, [Chen et al. , 1999 ] used clustering to identify ORF’ s that ha ve similar expression patterns to reduce the size of the regulatory netw ork to be inferred. Many clustering algorithms ha ve been proposed for gene expression data. F or example, [Eisen et al. , 1998 ] applied the a

verage-link hierarchical clustering algorithm to identify groups of co-regulated yeast genes. [Ben-Dor and Y

akhini, 1999 ] reported success with their CAST algorithm. [T amayo et al. , 1999 ] used self-or

ganizing maps to identify clusters in the yeast cell cycle and human hematopoietic dif fer

entiation data sets. Other techniques, such as principal component analysis (PCA), ha ve also been proposed to analyze gene expression data. Principal component analysis (PCA) ([Dunteman, 1989 ], [Everitt and Dunn, 1992 ], [Jollif fe, 1986 ]) is a classical technique to reduce the dimensionality of the data set by transforming to a new set of variables to summarize the features of the data set. In particular, [Raychaudhuri et al. , 2000 ] applied PCA to the sporulation data set 1

. Using dif ferent data analysis techniques and dif ferent clustering algorithms to analyze the same data set can lead to very dif ferent conclusions. F or example, [Chu et al. , 1998 ] identified seven clusters in the sporulation data set using the Cluster softw are [Eisen et al. , 1998 ], but [Raychaudhuri et al. , 2000 ] claimed that there are no clusters present in the same data set when the data points are viewed in the space of the first tw o principal components (PC’ s). In this paper, we empirically in vestigate the ef

fectiveness of PCA as a preprocessing step in cluster analysis using both real gene expression data sets with external clustering criteria and synthetic data sets. 

## 2 Principal Component Analysis (PCA) 2.1 An Example of PCA The central idea of principal component analysis (PCA) is to reduce the dimensionality of the data set while retaining as much as possible the variation in the data set. Principal components (PC’ s) are linear transformations of the original set of variables. PC’ s are uncorrelated and ordered so that the first few PC’ s contain most of the v ariations in the original data set [Jollif fe, 1986 ]. The first PC has the geometric interpretation that it is a new coordinate axis that maximizes the v ariation of the projections of the data points on the new coordinate axis. Figure 1 shows a scatterplot of some fictitious data points in tw o dimensions ( and  ). The points show an elliptical shape, and the first PC is in the direction of the principal axis of this ellipse (marked   in Figure 1). The second 1 Sporulation is the process in which diploid cells undergo meiosis to produce haploid cells in reproduction of yeast. The sporulation data set [Chu et al. , 1998 ] shows the temporal expression patterns of 97% of yeast genes over seven successive time points in the sporulation of yeast.2 2 PRINCIP AL COMPONENT AN

AL YSIS (PCA)
PC is orthogonal to the first PC and is marked   in Figure 1. If the data points are projected onto the first PC, most of the variation of the tw o dimensional data points w ould be captured in one dimension.

                         !! "" ## $$ % & ' ( ) * + , x

x PC PC

1

2 1 2

Figure 1: An example illustrating PCA 

### 2.2 Definitions of PCA experimental conditions (variables) p 1genes or clones (observations) 1 j

X

n i i,j

vector x j class

Figure 2: Notations for gene expression data Let - be a gene expression data set with . genes and / experimental conditions. In this report, our goal is to cluster the genes in the data set, hence the experimental conditions are the v ariables. Let 0 132 be a column vector of expression levels of all the . genes under experimental condition 4 (see Figure 2). A PC is a linear transformation of the experimental conditions. Let 0 57698;:=< >@?  BA3CED >

0 1F2 be the G th PC. In particular, the first PC, 0 5BH , can be written as :

< >@?  BA 

D >

0 1F2 . Let I be the cov ariance matrix of the data, and 0 A

6 be a column vector of all the AJCED > ’ s, i.e ., 0 ALK 6 8NM A3CED  PO

A3CED  EOPQPQPQRO ASCED <UT
. T o derive the first PC, we ha

ve to find 0 A

H that maximizes VXWZY M[:\< >]?  BA 

D >

0 1F2

T

8 0 A

K H

I^0 A H , subject to the constraint 0 A K H

0 A

H

8`_
. It can be shown that 0 A

H is the eigen

vector corresponding to the lar gest eigen

value, a

 , of I , and VXWZY M 0 5

H

T

8 a

 [Jollif fe, 1986 ]. In general, the G th PC, 0 5F6b8N:c< >]?  BASCRD >

0 1F2 , can be derived by maximizing VFWUY M :\< >@? 

ASCED >

0 1F2

T , such that 0 A

K 6

0 A

6d8e_ and 0 A

K 6

0 ALf 8hg , where i\jkG
. It can be shown that 0 A

6 is an eigen

vector of I corresponding to its G th largest eigen

value a

C , and VXWZY M 0 5Z6 T

8 a

C [Jollif fe, 1986 ].2.3 Choosing the number of PC’ s 3
In the case of gene expression data, the population co variance matrix I is not known. The sample co v ariance matrix can be used instead. Let D > be the gene expression level of gene i under experimental condition 4
. The sample co variance between conditions 4 and G , M 4

O

G

T , can be calculated as   

:   ?



M 

D > 
	 T

M

 D C


	 T , where 8  

:   ?





D > . From the derivation of PC’ s, the G th PC can be interpreted as the direction that maximizes the variation of the projections of the data points such that it is orthogonal to the first G

 _ PC’ s, and the G th PC has the G th largest variance among all PC’ s. Since most of the variation of high dimensional data points can be captured in reduced dimensions defined by the first few PC’ s, PCA is often used in visualization of high dimensional data points. 

### 2.3 Choosing the number of PC’ s Since the v ariance of PC’ s are ordered, usually the first  (  9/ , where / is the number of experiments in the data set) PC’ s are used in data analysis. The next question is how we should choose  , the number of first PC’ s to be retained, to adequately represent the data set. There are some common rules of thumb to choose the number of components to retain in PCA. Most of the rules are informal and ad-hoc. The first common rule of thumb is to choose  to be the smallest integer such that a chosen percentage of total variation is exceeded. In [Raychaudhuri et al. , 2000 ], the first tw o components which represent over 90% of the total variation in the sporulation data were chosen. Another common approach uses a scr

ee gr aph , in which the G th eigen value is plotted against the component number, G
. The number of components  is chosen to be the point at which the line in the scree graph is “steep” to the left but “not steep” to the right. The main problem with these approaches is that they are very subjective. There are some more formal approaches in the literature, but in practice, they tend not to w ork as well as the ad-hoc approach [Jollif fe, 1986 ]. 

### 2.4 Co variance versus correlation matrices In the PCA literature, some authors prefer to define PC’ s using the corr elation matrix instead of the co variance matrix. The correlation between a pair of variables is equivalent to the covariance divided by the product of the standard deviations of the tw o variables. Extracting the PC’ s as the eigen vectors of the correlation matrix is equivalent to computing the PC’ s from the original variables after each has been standardized to ha

ve unit variance. PCA based on covariance matrices has the potential drawback that the PC’ s are highly sensitive to the unit of measurement. If there are lar ge dif ferences between the variances of the v ariables, then the first few PC’ s computed with the covariance matrix are dominated by the variables with large v ariances. On the other hand, defining PC’ s with the correlation matrix has the drawback that the data is arbitrarily re-scaled to ha ve unit variance. The general rule of thumb is to define PC’ s using the correlation matrix if the variables are of dif ferent types [Jollif

fe, 1986 ]. 

### 2.5 Application of PCA in cluster analysis In the clustering literature, PCA is sometimes applied to reduce the dimension of the data set prior to clustering. The first few (say  ,  / ) PC’ s are usually used (for example, [Jollif fe et al. , 1980 ]). Most clustering algorithms require a measure of pairwise similarity or dissimilarity between obser

vations as input. There are tw o popular pairwise similarity metrics in clustering gene expression data: Euclidean distance (for example, [W en et al. , 1998 ]) and correlation coef ficient (for example, [Eisen et al. , 1998 ]). The pairwise Euclidean distance between tw o objects is unchanged after the PCA step if all / PC’ s are used. When Euclidean distance is used as the similarity metric, using the first4 2 PRINCIP AL COMPONENT AN

AL YSIS (PCA)

 PC’ s simply pro vides an approximation to the similarity metric [Jollif fe, 1986 ]. When correlation coef ficient is used as the similarity metric, the pairwise correlation coef ficient between tw o objects is not the same after the PCA step even if all / PC’ s are used. There is no simple relationship between the correlation coef

ficients of the same pair of objects with and without PCA. In general, the extra computation to find the PC’ s f ar outweighs any reduction in running time for using fewer PC’ s to compute the Euclidean distance [Jollif fe, 1986 ]. So, the hope for using PCA prior to cluster analysis is that PC’ s may “extract” the cluster structure in the data set. Figure 3 is a fictitious situation in which the PCA preprocessing step before cluster analysis may help. The first PC is in the direction of inter-cluster separation (the blue dotted line) in Figure 3. Projection of the data points on the first PC clearly highlights the separation between the tw o clusters in the data. However, PCA does not help in all situations. For example, in Figure 4, the first PC is in the direction of  . Projection of the data points onto the first PC does not preserve the separation between the tw o clusters in the data.

 

          

         

   !! " ## $$ %% && ' ( )) ** + ,

-------------------- .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... .................... 

### x

x 1 2

Figure 3: An example illustrating PCA helps

in cluster analysis. / 0 1 2

3 4 5 6 7 8 9 : ;; << = > ? @ A B CC DD
E F GG HH

I J KK L M N OO P QQ R S T U V W X

Y Z [[ \\ ]] ^^ __ ``

aa b cc d ee ff 

### x

x 1 2

Figure 4: An example illustrating PCA does not help in cluster analysis. In addition to the fictitious examples above illustrating the possible pros and cons of PCA on cluster analysis, [Chang, 1983 ] showed theoretically that the first few PC’ s may not contain cluster information under certain assumptions. Assuming that the data is a mixture of tw o multivariate normal distributions with dif ferent means but with an identical within-cluster covariance matrix, [Chang, 1983 ] derived a relationship between the distance of the tw o subpopulations and any subset of PC’ s, showing that the set of PC’ s with the largest eigen v alues does not necessarily contain more cluster structure information (the distance between the tw o subpopulations is used as a measure of discriminatory power for cluster structures). He also generated an artificial example in which there are tw o classes, and if the data points are visualized in tw o dimensions, the tw o classes are only well-separated in the subspace of the first and last PC’ s. In [Raychaudhuri et al. , 2000 ], PCA w as applied to the sporulation data set [Chu et al. , 1998 ]. The data points were visualized in the subspace of the first tw o PC’ s, and they showed a unimodal distribution. [Raychaudhuri et al. , 2000 ] concluded that the sporulation data may not contain any clusters based on visualization. W ith [Chang, 1983 ]’ s theoretic results and the possibility of the situation in Figure 4 in mind, it is clear that clustering with the PC’ s instead of the original variables does not ha ve universal success. However, the theoretical results in [Chang, 1983 ] are true only under an unrealistic assumption for gene expression data ( i.e ., there are tw o classes and each of the classes is generated according to the multivariate normal distribution with a common covariance matrix). Therefore, there is a need to in

vestigate the ef

fectiveness of PCA as a preprocessing step to cluster analysis on gene5

expression data before any conclusions are drawn. Our report is an attempt for such an empirical study
. 

## 3 Ov erview of Our Methodology Our goal is to empirically in vestigate the ef

fectiveness of clustering gene expression data using PC’ s instead of the original v ariables. Our methodology is to run a clustering algorithm on a given data set, and then apply the same algorithm to the PC’ s of the same data set. Then, the clustering results with and without PCA are compared against an external criterion. The details of the experiments will be discussed in the following sections. 

### 3.1 Data sets W e used tw o gene expression data sets which ha ve external criteria, and four sets of synthetic data to ev aluate the ef

fectiveness of PCA. In this report, we use the w ord class to refer to a group in our external criteria that is used to assess clustering results. The w ord cluster is used to refer to clusters obtained by a clustering algorithm. W e assume both classes and clusters are disjoint. 3.1.1 Gene expression data sets The ovary data: A subset of the o v ary data set ([Schummer et al. , 1999 ], [Schummer, 2000 ]) is used. The o v ary data set is generated by hybridizing randomly selected cDNA ’ s to membrane arrays. The subset of the o v ary data set we used contains 235 clones (clones are portions of genes) and 24 samples, 7 of which are derived from normal tissues, 4 from blood samples, and the remaining 13 from o varian cancers in v arious stages of malignancy. The tissue samples are the experimental conditions. The 235 clones were sequenced, and they correspond to 4 dif ferent genes. The numbers of clones corresponding to each of the four genes are 58, 88, 57, and 32 respectively
. W e expect clustering algorithms to separate the four dif

ferent genes. Hence, the four genes form the four classes of external criteria for this data set. Dif ferent clones may ha

ve dif ferent hybridization intensities. Therefore, the data for each clone is normalized across the 24 experiments to ha ve mean 0 and variance 1 2

. The yeast cell cycle data: The second gene expression data set we used is the yeast cell cycle data set [Cho et al. , 1998 ] which shows the fluctuation of expression levels of approximately 6000 genes over tw o cell cycles (17 time points). By visual inspection of the raw data, [Cho et al. , 1998 ] identified 420 genes which peak at dif

ferent time points and categorized them into fi ve phases of cell cycle. Out of the 420 genes they classified, 380 genes were classified into only one phase (some genes peak at more than one phase in the cell cycle). Since the 380 genes were identified by visual inspection of gene expression data according to the peak times of genes, we expect clustering results to correspond to the fi ve phases to a certain degree. Hence, we used the 380 genes that belong to only one class (phase) as our external criteria. The 17 time points are the experimental conditions. The data is normalized to ha ve mean 0 and v ariance 1 across each cell cycle as suggested in [T amayo et al. , 1999 ]. 2 If the correlation matrix is used instead of the covariance matrix in PCA, the tissue samples (experiments) are the variables and are standardized across all the clones.6 3 O VER VIEW OF OUR METHODOLOGY 3.1.2 Synthetic data sets Since the array technology is still in its inf ancy [Lander, 1999 ], the “real” data may be noisy, and clustering algorithms may not be able to extract all the classes contained in the data. There may also be information in real data that is not known to biologists. Furthermore, synthetic data sets pro vide us with inexpensive replicates of the data to increase the reliability of our empirical study
. Therefore, we w ould like to complement our empirical study of the ef fectiveness of PCA with synthetic data, for which the classes are known. T o the best of our knowledge, modeling gene expression data sets is an ongoing ef fort by many researchers, and there is no well-established model to represent gene expression data. The following four sets of synthetic data represent our preliminary ef forts on synthetic gene expression data generation. W e do not claim that any of the four sets of synthetic data capture most or all of the characteristics of gene expression data. Each of the synthetic data has strengths and weaknesses. By using all four sets of synthetic data to ev aluate the ef fectiveness of PCA on clustering, we hope to achieve a thorough comparison study capturing many dif ferent aspects of gene expression data. The first tw o synthetic data sets represent attempts to generate replicates of the ovary data set by randomizing dif

ferent aspects of the original data. The last tw o synthetic data sets are generated by modeling expression data with a mathematical formula. In each of the four synthetic data sets, ten replicates are generated. Ideally, more replicates w ould be more desirable. However, the algorithms ha ve very long running time 3

. In each replicate, 235 observ ations and 24 variables are randomly generated. Mixture of normal distrib utions on the ovary data: V isual inspection of the o vary data suggests that the data is not too f

ar from normal. The expression levels for dif ferent clones of the same gene are not identical due to the f

act that the clones represent dif ferent portions of the cDNA. Figure 5 shows the distribution of the expression levels in a normal tissue in a dif ferent class (gene) from the o vary data. W e found that the distributions of the normal tissue samples are typically closer to normal distributions than those of tumor samples, for example, Figure 6. Even though some of the tumor tissues from some classes (genes) do not closely follow the normal distribution, we generate the data using a mixture of multivariate normal distributions in this synthetic data set. The sample co variance matrix and the mean vector of each of the four classes (genes) in the o vary data are computed. The size of each class in the synthetic data is the same as the ovary data. Each class in the synthetic data is generated according to a multivariate normal distribution with the sample co variance matrix and the mean vector of the corresponding class in the ovary data. This synthetic data set preserves the covariance between the tissue samples in each gene. It also preserves the mean vectors of each class. The weakness of this synthetic data set is that the assumption of the underlying multivariate normal distribution for each class may not be true for real genes. Randomly permuted ovary data: No underlying distribution of the gene expression data is assumed in this synthetic data set. The size of each class in this synthetic data set is again the same as the ovary data. The random data for an artificial gene in class (where 8 _

OPQPQPQRO  ) under experimental condition 4 (where 4 8 _

OPQPQPQRO  ) can be generated by randomly sampling (with replacement) the expression levels under condition 4 in the same class of the ovary data. This data set does not assume any underlying distribution. However, any possible correlation between tissue samples (for example, the normal tissue samples may be correlated) is not preserved due to the independent random sampling of the expression levels from each experimental condition. 3 It takes approximately 3 to 4 hours to run the modified greedy algorithm (see Section 3.5) with one clustering algorithm on one replicate on a Pentium 500.3.1 Data sets 7

-1.5 -1.0 -0.5 0.0 0.5 1.00246810

expression levels of a normal tissue in class 3frequency

Figure 5: Histogram of the distribution of the expression levels in a normal tissue for a gene (class) in the ovary data

-1.5 -1.0 -0.5 0.0 0.5 1.0024681012

expression levels of a tissue sample in class 1frequency

Figure 6: Histogram of the distribution of the expression levels in a tumor tissue for a gene (class) in the ovary data8 3 O VER VIEW OF OUR METHODOLOGY Hence, the resulting sample co variance matrix of this randomly permuted data set w ould be close to diagonal. However, inspection of the ovary data shows that the sample covariance matrices are not too f ar from diagonal. Therefore, this set of randomly permuted data may be reasonable replicates of the original ovary data set. Cyclic data with different class sizes: This synthetic data set models ten cyclic classes. The obser

vations can be interpreted as genes and the experimental conditions can be interpreted as time points. In this data set, the cyclic beha vior of genes (observations) is modeled by the function. Classes are modeled as genes that ha ve similar peak times over the time series data. Dif ferent classes ha

ve dif ferent phase shifts and ha ve dif ferent sizes. Let  D > be the simulated expression level of gene i and condition 4 in this simulated data set with ten classes. Let  D >

8  >  a

> M A



   3M i

O 4

T

T , where SM i

O 4

T

8 M

 

> 

  

C  

T [Zhao, 2000 ]. A

 represents the a

verage expression level of gene i , which is chosen according to the normal distribution with mean 0 and standard deviation 2.  is the amplitude control for gene i , which is chosen according to a normal distribution with mean 3 and standard deviation 0.5. SM i

O 4

T models cyclic time series data. In this synthetic data set, each cycle is assumed to span 8 time points. G is the class number, which is chosen according to Zipf ’ s Law [Zipf, 1949 ], which allows us to model classes with dif ferent sizes. Since dif ferent classes ha

ve dif ferent values of G , dif ferent classes are represented by dif ferent phase shifts of the function. a

> is the amplitude control of condition 4 , chosen according to the normal distribution with mean 3 and standard deviation 0.5. 

> represents an additive experimental error, chosen according to the standard normal distribution. Each observ ation (row) is normalized to ha ve mean 0 and variance 1 before PCA or any clustering algorithm is applied. This synthetic model suf fers from the drawback that the form of the data may not be realistic for gene expression data, and the ad-hoc choice of the parameters for the distributions of A

 ,  , a

> , and  > . However, there is evidence that the form of the model is reasonable for the yeast cycle data [Zhao, 2000 ]. Cyclic data with spiky classes: This synthetic data has the same form as the cyclic data with dif ferent class sizes. Again, there are ten synthetic classes, and the data are generated from the equation for  D > . However, the class number ( G in SM i

O 4

T ) is generated according to the uniform distribution, and not Zipf ’ s law . Hence, the class sizes are approximately the same in this data. Some genes show “spik y” beha vior, i.e ., their expression levels are changed sharply over a short period of time. W e use the term spik y classes to refer to classes in which genes show spik y beha vior
. Spik

y classes, are modeled by raising the function in 3M i

O 4

T to higher powers. Thus, dif ferent classes are modeled by phase shifts or dif ferent “spikiness” ( i.e ., dif ferent powers of the function in SM i

O 4

T ). This synthetic data set suf fers the same drawback as the cyclic data with dif ferent class sizes: the model may not be realistic, and the ad-hoc choice of the parameters for the distributions of A

 ,  ,

a

> , and  > . However, this synthetic data set hopes to capture a more complicated form of the real expression data by modeling classes with dif ferent shapes ( i.e ., dif

ferent spikiness) in addition to dif ferent peak times. 

### 3.2 PCA on the data sets W e use the co variance matrix to define PC’ s for both real and synthetic data. There are tw o main reasons for that. First, the variables in our case are the array experiments and hence are of the same type. In particular, for the ov ary data described in Section 3.1, all the experiments are scaled to ha ve the same median prior to PCA. Second, we computed the PC’ s using both the correlation and the co variance matrices for the ov ary and yeast cell cycle data and there is no major dif ference between3.3 Agreement between tw o partitions 9
the tw o sets of PC’ s. The first 14 PC’ s defined using the covariance matrix of the ovary data account for over 90% of the total v ariation of the data. On the yeast cell cycle data, the first 8 PC’ s defined using the covariance matrix account for o ver 90% of the total variation. The scree graphs for these tw o data sets are shown in Figures 21 and 22 in Appendix A. 

### 3.3 Agreement between two partitions In order to compare clustering results against external criteria, a measure of agreement is needed. Since we assume that both the external criteria and clustering results are disjoint, measures of agreement between tw o partitions can be used. In the statistics literature, many measures of agreement were proposed and ev aluated (for example, [Rand, 1971 ], [Milligan et al. , 1983 ],[Hubert and Arabie, 1985 ], [Milligan and Cooper, 1986 ] and many others). Given a set of . objects 8 

OPQPQPQEO 

  , suppose  8  

OPQPQPQ O

	  and 8 V



OPQPQPQRO V  represent tw o dif

ferent partitions of the objects in such that   ?







8 8 >@?  V

> and 

  

  8  8 V

>



V

>

 for _

 i  8 i    and _

 4  8 4   
. Suppose that  is our external criterion and is a clustering result. Let W be the number of pairs of objects that are placed in the same class in  and in the same cluster in ,  be the number of pairs of objects in the same class in  but not in the same cluster in , be the number of pairs of objects in the same cluster in but not in the same class in  , and  be the number of pairs of objects in dif ferent classes and dif

ferent clusters in both partitions. The quantities W and  can be interpreted as agreements, and  and as disagreements. The Rand index [Rand, 1971 ] is simply  	!"	# . The Rand index lies between 0 and 1. When the tw o partitions agree perfectly, the Rand index is 1. The problem with the Rand index is that the expected value of the Rand index of tw o random partitions does not take a constant value (say zero). The adjusted Rand index proposed by [Hubert and Arabie, 1985 ] assumes the generalized hyper geometric distribution as the model of randomness, i.e ., the  and partitions are picked at random such that the number of objects in the classes and clusters are fixed. Let .

 > be the number of objects that are in both class 

 and cluster V

> . Let .

 %$ and .

$ > be the number of objects in class 

 and cluster V

> respectively
. The notations are illustrated in T able 1. Class or Cluster V

 V

 QPQPQ V Sums   .



 .

  QPQPQ .

 .



$   .



 .



 QPQPQ .

 .



$

. . . .

. . . . . .

. . . . . 	 .



 .



 QPQPQ .

 .



$ Sums .

$  .

$  QPQPQ .

$ .

$&$ Z8 . T able 1: Notation for the contingency table for comparing tw o partitions. The general form of an index with a constant expected value is  

(' 

'

<

')#+*,'-  

.' / 

 /102/  

(' 

'

<

')#+*,'-  

.' , which is bounded abo ve by 1, and takes the value 0 when the index equals its expected v alue. Under the generalized hyper geometric model, it can be shown [Hubert and Arabie, 1985 ] that: 354 6)7  D > 98 .

 >  ;:=< > 8 4 6-7  8

.

 %$  ?: 7 > @8 .

$ >  A:=< >CB 8

.  D: (1)10 3 O VER VIEW OF OUR METHODOLOGY
The expression W

  can be simplified to a linear transformation of :  D >

8

.

 >  ;:

. W ith simple algebra, the adjusted Rand index [Hubert and Arabie, 1985 ] can be simplified to: :  D >

8

.

 > 

:  : 

8

.

 %$  ?: :

>

8

.

$ > 

:  B

8

.  D:

  : 

8

.

 %$ 

:

 :

>

8

.

$ > 

:   : 

8

.

 %$ 

:

:

>

8

.

$ > 

:  B

8

. 

: (2) Example 1 in Appendix B illustrates how the adjusted Rand index is computed. The Rand index for comparing the tw o partitions in Example 1 is 0.711, while the adjusted Rand index is 0.313. The Rand index is much higher than the adjusted Rand index, which is typical. Since the Rand index lies between 0 and 1, the expected v alue of the Rand index (although not a constant value) must be greater than or equal to 0. On the other hand, the expected value of the adjusted Rand index has v alue zero and the maximum v alue of the adjusted Rand index is also 1. Hence, there is a wider range of v alues that the adjusted Rand index can take on, thus increasing the sensitivity of the index. In [Milligan and Cooper, 1986 ], many dif ferent indices were evaluated for measuring agreement between tw o partitions in hierarchical clustering analysis across dif ferent hierarchy levels ( i.e ., with dif ferent numbers of clusters), and they recommended the adjusted Rand index as the index of choice. In this report, we adopt the adjusted Rand index as our measure of agreement between the external criteria and clustering results. 

### 3.4 Clustering algorithms and similarity metrics W e implemented three clustering algorithms: the Cluster Af finity Sear

ch T echnique (CAST) [Ben-Dor and Y akhini, 1999 ], the hierarchical aver

ag e-link algorithm, and the k-means algorithm (with a verage-link initialization) [Jain and Dubes, 1988 ]. CAST : W

e implemented the pseudo-code for CAST given in [Ben-Dor and Y akhini, 1999 ] with tw o additional heuristics that ha ve been added to B I O C L U S T , the implementation of CAST by its authors. The CAST algorithm takes as input the pairwise similarities of objects and a parameter  which is a real number between 0 and 1. The parameter  is a similarity threshold to decide whether an object is added to or removed from a cluster
. Thus, varying the parameter  changes the number of clusters formed. Please refer to Appendix C for more details of the algorithm. Hierarchical average-link: Agglomerative hierarchical algorithms build clusters bottom up. Initially, each object is in its own cluster
. In each step, the tw o clusters with the greatest cluster similarity are merged. This process is repeated until the desired number, G , of clusters is produced. In a verage-link, the cluster similarity criterion is the a verage pairwise similarity between objects in the tw o clusters. Refer to [Jain and Dubes, 1988 ] and [Anderber g, 1973 ] for detailed discussions on hierarchical algorithms. The a

verage-link clustering algorithm is used by [Eisen et al. , 1998 ] to analyze gene expression data. K-means: The number of clusters, G , is an input to the k-means clustering algorithm. Clusters are described by centr

oids , which are cluster centers, in the algorithm. In our implementation of k-means [Jain and Dubes, 1988 ], the initial centroids consist of the clustering results from a verage-link. Each object is assigned to the centroid (and hence cluster) with the closest Euclidean distance. New centroids of the G clusters are computed after all objects are assigned. The steps of assigning objects to centroids and computing new centroids are repeated until no objects are moved between clusters. Similarity metrics: There are tw o popular similarity metrics used in the gene expression analysis community: Euclidean distance (for example, [W en et al. , 1998 ]) and correlation coef ficient (for3.5 Our approach 11

example, [Eisen et al. , 1998 ]). In our experiments, we evaluated the ef fectiveness of PCA on cluster

ing analysis with both Euclidean distance and correlation coef ficient, namely, CAST with correlation coef ficient4 , a verage-link with both correlation and distance, and k-means with both correlation and distance. If Euclidean distance is used as the similarity metric, the minimum number of components in sets of PC’ s ( 

 ) considered is 2. If correlation is used, the minimum number of components ( 

 ) considered is 3 because there are at most 2 clusters if 2 components are used (when there are 2 components, the correlation coef ficient is either 1 or -1, see Appendix D for details). 

### 3.5 Our approach Given a data set with an external criterion, our evaluation methodology consists of the following steps: 1. A clustering algorithm is applied to the given data set, and the adjusted Rand index with the external criterion is computed. 2. PCA is applied to the same data set. The same clustering algorithm is applied to the data after the PCA preprocessing step using the first  PC’ s (where  8  

OPQPQPQ O / , and 

 is the number of components we start with, which is either 2 or 3 as explained in Section 3.4, and / is the number of experimental conditions in the data). The adjusted Rand index is computed for each of the clustering results using the first  PC’ s. 3. The same clustering algorithm is applied to the data after the PCA preprocessing step using dif ferent sets of PC’ s. The objective in this step is to find a set of PC’ s that gives a “high” adjusted Rand index. One w ay to determine the set of PC’ s that gives the maximum adjusted Rand index is by exhaustive search. However, exhaustive search is very computationally intensive: for  components, there are 8 / 

: possible sets of PC’ s that we ha ve to cluster
. The ovary data has 24 experimental conditions and when  8N_  , the exhaustive search gives us 8

 _



: (approximately 2.7 million) possible sets of PC’ s to cluster
. Since our objective is to show that the highest adjusted Rand index is not necessarily achieved by the first  PC’ s, it suf

fices to show that there exists a set of PC’ s that achieves higher adjusted Rand index than the first PC’ s. A simple strate gy we implemented is the gr eedy approach. Let   be the minimum number of components that we start with. In the greedy approach, we start with the exhaustive search for the minimum number of components,   . Denote the optimum set of components as / . For each  (where  8eM 



 _

T

OPQPQPQEO / ), one additional component that is not already in /



 is added to the set of components, the data with all the . genes under this set of components is clustered, and the adjusted Rand index is computed. The additional component that achieves the maximum adjusted Rand index is added to /



 to form / . In the greedy approach, we implicitly assume that a set of PC’ s that achieves a high adjusted Rand index for  components is a good candidate for achieving a high adjusted Rand index for   _ components (for  8  

OPQPQPQ O

M /

 _

T ). Since the assumption for the greedy approach may not be satisfied, we implemented a modified gr eedy approach. The modified greedy approach requires another parameter, G , which is an inte ger indicating the number of “best” solutions to keep in each search step. Denote the optimum G sets of components as / 8  /

OPQPQPQ O C /

 , where  8 



OPQPQPQ O / . The modified greedy approach 4 When Euclidean distance is used in CAST , the algorithm usually does not con verge in practice.12 4 RESUL TS AND DISCUSSION
also starts with an exhaustive search for the minimum number of components,   . However, G sets of components which achieve the top G adjusted Rand indices were stored. For each  (where  8 M 





_

T

OPQPQPQEO / ) and each of the  / (where i

8 _

OPQPQPQEO G ), one additional component that is not already in  /



 is added to the set of components, the subset of data with the extended set of components is clustered, and the adjusted Rand index is computed. The top G sets of  components that achieves the highest adjusted Rand indices are stored in 

/ . The modified greedy approach allows the search to ha ve more choices in searching for a set of components that gives a high adjusted Rand index. Note that when G

8 _ , the modified greedy approach is identical to the simple greedy approach, and when G 8 8

/ 

:
, the modified greedy approach is reduced to exhaustive search. So the choice for G is a tradeof f between running time and quality of solution. In our experiments, G is set to be 3. 

## 4 Results and Discussion W e ran our experiments on tw o gene expression data sets and four synthetic data sets. In this section, the results of the experiments will be presented. Before the detailed results are presented for each set of experiments, here are our overall conclusions from our empirical study: W e found that the PCA preprocessing step does not necessarily impro ve cluster quality, i.e ., the adjusted Rand indices of the clustering results on the data after PCA are not necessarily higher than the adjusted Rand indices of the clustering results on the original data on both real and synthetic data. W e also showed that in most cases, the first  components (where  8  

OPQPQPQRO / ) do not necessarily give the highest adjusted Rand index, i.e ., there exists another set of  components that achieves a higher adjusted Rand index than the first  components. There are no clear trends regarding the choice of the optimal number of PC’ s over all the data sets and over all the clustering algorithms and over the dif ferent similarity metrics. There is no ob vious relationship between cluster quality ( i.e ., adjusted Rand index) and the number or set of PC’ s used. In most cases, the modified greedy approach achieves higher adjusted Rand indices than the simple greedy approach. In the following sections, the detailed experimental results on each data set is presented. F or some of the results, graphs plotting the adjusted Rand index against the number of components are shown. Usually the adjusted Rand index without PCA, the adjusted Rand index of the first  components, and the adjusted Rand indices using the greedy and modified greedy approaches are shown in each graph. Note that there is only one value for the adjusted Rand index computed with the original v ariables (without PCA), while the adjusted Rand indices computed using PC’ s vary with the number of components. 

### 4.1 Gene expression data 4.1.1 The ovary data Figure 7 shows the result of our experiments on the ovary data using CAST [Ben-Dor and Y akhini, 1999 ] as the clustering algorithm and correlation coef ficient as the similarity metric. The adjusted Rand in-4.1 Gene e xpression data 13

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3 	

	


	            

F i r s t  !

" # g r e e d y $%

& ' ( ' ) &+*

, ) ) & - Figure 7: Adjusted Rand index against the number of components using CAST and correlation on the ovary data. dices using the first  components (where  8 /. OPQPQPQEO  ) are mostly lower than that without PCA. However, the adjusted Rand indices using the greedy and modified greedy approaches for 4 to 22 components are higher than that without PCA. This shows that clustering with the first  PC’ s instead of the original v ariables may not help to extract the clusters in the data set, and that there exist sets of PC’ s (other than the first few which contain most of the variation in the data) that achieve higher adjusted Rand indices than without PCA. Moreo ver, the adjusted Rand indices computed using the greedy and modified greedy approaches are not very dif ferent for this data set using the CAST algorithm and correlation.

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3 02123547698:<;>=9:735?2:207690A@CB D E F G H I J E K L M E

N OQP

R S g r e e d y f i r s t T%O

U V W V X UZY

[ X X U \ Figure 8: Adjusted Rand index against the number of components using k-means and correlation on the ovary data. Figures 8 and 9 show the adjusted Rand indices using the k-means algorithm on the ovary data with the correlation and Euclidean distance as similarity metrics respectively
. Figure 8 shows that the adjusted Rand indices using the first  components tends to increase from below the index without PCA to abo ve that without PCA as the number of components increases. However, the results using the same algorithm but Euclidean distance as the similarity metric show a very dif ferent picture (Figure 9): the adjusted Rand indices are high for first 2 and 3 PC’ s and then drop drastically to below that without PCA. Manual inspection of the clustering result of the first 4 PC’ s using k-means and Euclidean distance shows that tw o classes are combined in the same cluster while the clustering result of the first 3 PC’ s separates the 4 classes, showing that the drastic drop in the adjusted Rand index reflects de gradation of cluster quality with additional PC’ s. When the data points are visualized in the space of the first and second PC’ s, the four classes are reasonably well-separated in the Euclidean space. However, when the data points are visualized in the space of the second and fourth PC’ s, tw o classes overlap. The degradation of cluster quality with additional PC’ s is probably because classes are not14 4 RESUL TS AND DISCUSSION

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5

2 4 6 8 1 0 1 2 1 4 1 6 1 8 2 0 2 2 2 4 	

	
             

f i rs t   

! " g r e e d y #$

% & ' & ( %*)

+ ( ( % , Figure 9: Adjusted Rand index against the number of components using k-means and Euclidean distance on the ov ary data.

very well-separated in the Euclidean space of the higher PC’ s, and hence, it is more dif ficult for the clustering algorithms to extract the classes (see Appendix E). Figures 8 and 9 also show that dif ferent similarity metrics may ha ve very dif

ferent ef fect on the use of PCA as a preprocessing step to cluster analysis. The adjusted Rand indices using the modified approach in Figure 8 show an irregular pattern. In some instances, the adjusted Rand index computed using the modified greedy approach is even lower than that using the first few components and that using the greedy approach. This shows that our heuristic assumption for the greedy approach is not alw ays valid, i.e ., a set of PC’ s that achieve a high adjusted Rand index for  components may not be a good candidate for achieving a high adjusted Rand index for   _ components (for  8  

OPQPQPQ O

M /

 _

T ). Nevertheless, the greedy and modified greedy approaches show that there exists other sets of PC’ s that achieve higher adjusted Rand indices than the first few PC’ s most of the time. The results using the hierarchical a verage-link algorithm with correlation coef ficient and Euclidean distance as similarity metrics show a similar pattern to the results using k-means (graphs not shown here). Note that the adjusted Rand index without PCA using CAST and correlation (0.664) is much higher than that using k-means (0.563) and a verage-link (0.572) using the same similarity metric. Manual inspection of the clustering results without PCA shows that only CAST clusters mostly contain clones from each class, while k-means and a verage-link clustering results combine tw o classes into one cluster
. This confirms that higher adjusted Rand indices reflect higher cluster quality with respect to the external criteria. W

ith the first  components, CAST with correlation has a similar range of adjusted Rand indices to the other algorithms (approximately between 0.55 to 0.68). The rule of thumb of choosing the first 14 PC’ s to co ver 90% of the total variation in the data w ould ha ve a detrimental ef

fect on cluster quality if CAST with correlation, k-means with distance, or a verage-link with distance is the algorithm being used. When correlation is used (Figures 7 and 8), the adjusted Rand index using all 24 PC’ s is not the same as that using the original variables. On the other hand, when Euclidean distance is used (Figure 9), the adjusted Rand index using all 24 PC’ s is the same as that with the original variables. This is because the Euclidean distance between a pair of genes using all the PC’ s is the same as that using the original variables. But correlation coef ficient is not preserved after PCA (Section 2.5).4.2 Synthetic data 15

0
. 3 0
. 3 5 0
. 4 0
. 4 5 0
. 5 0
. 5 5

3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 	

	

	             

 

 g r e e d y f i r s t !"

# $ % $ & #('

) & & # * Figure 10: Adjusted Rand index against the number of components using CAST and correlation on the yeast cell cycle data. 4.1.2 The yeast cell cycle data Figure 10 shows the result on the yeast cell cycle data using CAST [Ben-Dor and Y akhini, 1999 ] as the clustering algorithm and correlation coef ficient as the similarity metric. The adjusted Rand indices using the first 3 to 7 components are lower than that without PCA, while the adjusted Rand indices with the first 8 to 17 components are comparable to that without PCA.

0
. 4 0
. 4 2 0
. 4 4 0
. 4 6 0
. 4 8 0
. 5 0
. 5 2 0
. 5 4

2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 	

	
	             

 !

" # g r e e d y f i r s t $%

& ' ( ' ) &+*

,) ) & - Figure 11: Adjusted Rand index against the number of components using k-means and Euclidean distance on the yeast cell cycle data. Figure 11 shows the result on the yeast cell cycle data using k-means and Euclidean distance. The adjusted Rand indices without PCA are relatively high compared to those using PC’ s. Figure 11 on the yeast cell cycle data shows a very dif ferent picture than Figure 9 on the ovary data. This shows that the ef fectiveness of clustering with PC’ s depends on the data set being used. The results on the yeast cell cycle data sets using k-means with correlation and a verage-link (with both correlation and Euclidean distance) are not shown here. 

### 4.2 Synthetic data 4.2.1 Mixture of normal distrib utions on the ovary data Figure 12 shows the results of our experiments on the synthetic mixture of normal distributions on the ovary data using CAST [Ben-Dor and Y akhini, 1999 ] as the clustering algorithm and correlation coef ficient as the similarity metric. The lines in Figure 12 represent the a verage adjusted Rand indices over the 10 replicates of the synthetic data, and the error bars represent one standard deviation from the mean for the modified greedy approach and for using the first  PC’ s. The error bars show that16 4 RESUL TS AND DISCUSSION

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5 0
. 8 0
. 8 5

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3 	

	

	
            

N O P R S g r e e d y F i r s t T%O

U V W V X UZY

[ X X U \ Figure 12: Adjusted Rand index against the number of components using CAST and correlation on the mixture of normal distributions on the ovary data. the standard deviations using the modified greedy approach tend to be lower than that using the first  components. A careful study also shows that the modified greedy approach has lower standard deviations than the greedy approach (data not shown here). This shows that the modified greedy approach is more robust than the greedy approach in identifying a set of components with a high adjusted Rand index. The error bars for the case without PCA are not shown for clarity of the figure. The standard deviation for the case without PCA is 0.064 for this set of synthetic data, which w ould overlap with those using the first components and the modified greedy approach. A manual study of the experimental results from each of the 10 replicates (details not shown here) shows that 8 out of the 10 replicates show very similar patterns to the a verage pattern in Figure 12, i.e ., most of the cluster results with the first  components ha

ve lower adjusted Rand indices than that without PCA, and the results using the greedy and modified greedy approach are slightly higher than that without PCA. In the following results, only the a verage patterns will be shown. Figure 12 shows a similar trend to real data in Figure 7, but the synthetic data has higher adjusted Rand indices for the clustering results without PCA and with the greedy and modified greedy approaches. As in the case with real data (Figure 7), the adjusted Rand indices with the first PC’ s lie below that without PCA, and those with the greedy and modified greedy approach for 4 to 22 components are abo ve that without PCA.

0
. 3 0
. 3 5 0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5 0
. 8

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3  "!$#&%('*),+"-.)/!$01)/2&%(2&3 4 5 6 7 8 9 : ; 6 < = > 6

N OQP

R S g r e e d y F i r s t T%O

U V W V X UZY

[ X X U \ Figure 13: A

verage adjusted Rand index against the number of components using k-means and correlation on the mixture of normal distributions of the ovary data. The a verage adjusted Rand indices using the k-means algorithm with the correlation and Euclidean distance as similarity metrics are shown in Figure 13 and Figure 14 respectively
. In Figure 13, the adjusted Rand indices using the first  components gradually increase as the number of components increases, while in Figure 14, the adjusted Rand indices using the first  indices are mostly below that4.2 Synthetic data 17

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7

2 4 6 8 1 0 1 2 1 4 1 6 1 8 2 0 2 2 2 4 2,"!$#&%('*),+"-.)/!$01)/2&%(2&3 4 5 6 7 8 9 : ; 6 < = > 6

N OQP

R S g r e e d y f i r s t T%O

U V W V X UZY

[ X X U \ Figure 14: A

verage adjusted Rand index against the number of components using k-means and Euclidean distance on the mixture of normal distributions of the ovary data. without PCA. The results using a

verage-link and correlation (not shown here) are similar to that of k-means and correlation. The a

verage adjusted Rand indices using a verage-link and Euclidean distance (not shown here) tend to decrease as the number of components is increased. 4.2.2 Randomly permuted ovary data

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5 0
. 8 0
. 8 5 0
. 9

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3 07123547698:A; =9:235?7:207690A@CB         

N O P R S g r e e d y f i r s t T%O

U V W V X U!Y

[ X X U \ Figure 15: A

verage adjusted Rand index against the number of components using CAST and correlation on the randomly permuted o vary data. Figures 15, 16, and 17 show the a verage adjusted Rand indices using CAST and correlation, k-means with correlation, and k-means with Euclidean distance on the randomly permuted o vary data respectively
. The general trend is very similar to the results on the mixture of normal distributions in Section 4.2.1. The a

verage adjusted Rand indices computed from CAST clusters in Figure 15 using the first  PC’ s lie below that using the original variables (similar to Figure 12). When the k-means algorithm is used with correlation as the similarity metric in Figure 16, the a verage adjusted Rand indices using the first  PC’ s tend to increase as the number of components increases (similar to Figure 13). When the modified greedy or the greedy approach is used, the a verage adjusted Rand indices for all clustering results (except CAST with 3 components) with all of the algorithms are abo ve that without PCA.18 4 RESUL TS AND DISCUSSION

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5 0
. 8 0
. 8 5 0
. 9

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3 02123547698:<;>=9:735?2:207690A@CB         

N OQP

R S g r e e d y f i r s t T%O

U V W V X UZY

[ X X U \ Figure 16: A

verage adjusted Rand index against the number of components using k-means and correlation on the randomly permuted ov ary data.

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5 0
. 8 0
. 8 5

2 4 6 8 1 0 1 2 1 4 1 6 1 8 2 0 2 2 2 4    &      

N OQP

R S g r e e d y f i r s t T+O

U V W V X U!Y

[ X X U \ Figure 17: A

verage adjusted Rand index against the number of components using k-means and Euclidean distance on the randomly permuted ovary data. 4.2.3 Cyclic data with different cluster sizes The results using this synthetic data set are very dif ferent than those using the mixture of normal distributions and the randomly permuted ovary data. The a verage adjusted Rand index without PCA (0.941) is very high. The high adjusted Rand index indicates that the data set is very clean and the clustering algorithms reco ver most of the classes. Figure 18 shows the a

verage adjusted Rand indices using CAST and correlation. Manual inspection shows that in 7 out of the 10 replicates, the adjusted Rand index without PCA is perfect ( i.e ., 1). The a verage adjusted Rand indices using the first  components are much lower than that without PCA in Figure 18. Note that there is sharp decline in the a verage adjusted Rand index when all 24 components are used in the greedy approach in Figure 18. This is no accident. A close inspection shows that 9 out of the 10 replicates show a drastic decline in adjusted Rand index when all the components are used with the greedy approach. In all of the 9 replicates, the additional component which contributes to the sharp decline in the adjusted Rand index is the third PC. The results of the modified greedy approach are not shown since the greedy approach achieves almost perfect adjusted Rand indices. Figure 19 shows the a

verage adjusted Rand indices with the k-means algorithm and Euclidean distance as the similarity metric. Again, the adjusted Rand index without PCA is very high. But in this case, using the first  components achieve higher or comparable adjusted Rand indices to that without PCA.4.2 Synthetic data 19

0
. 4 0
. 5 0
. 6 0
. 7 0
. 8 0
. 9 1 1
. 1

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3 
	A
	


	         

 g r e e d y f i rs t Figure 18: A

verage adjusted Rand index against the number of components using CAST and correlation on the cyclic data with dif ferent cluster sizes.

0
. 9 9 0
. 9 9 2 0
. 9 9 4 0
. 9 9 6 0
. 9 9 8 1 1
. 0 0 2

2 4 6 8 1 0 1 2 1 4 1 6 1 8 2 0 2 2 2 4 2,"!$#&%('*)"+,-.)/!$01)/2&%(213 4            

N OQP

R S g r e e d y f i r s t Figure 19: A

verage adjusted Rand index against the number of components using k-means and Euclidean distance on the cyclic data with dif ferent cluster sizes. 4.2.4 Cyclic data with spiky clusters The general trend using this set of synthetic data sets is very dif ferent than the other synthetic data sets (even the data set in Section 4.2.3): using PCA (with the first components or the greedy or modified greedy approach) helps to achieve higher adjusted Rand indices.

0
. 4 0
. 4 5 0
. 5 0
. 5 5 0
. 6 0
. 6 5 0
. 7 0
. 7 5 0
. 8

3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3 ! #"%$'&()*&+!,#&-#"%#. / 0 1 2 3 4 5 6 1 7 8 9 1

N O P R S g r e e d y f i r s t T%O

U V W V X U!Y

[ X X U \
Figure 20: A

verage adjusted Rand index against the number of components using CAST and correlation on the cyclic data with spik y clusters. Figure 20 shows an example of results applying CAST and correlation to this set of synthetic data. Unlike the results with other real or synthetic data sets, the adjusted Rand indices of CAST clusters with the first  PC’ s lie abo ve that without PCA in Figure 20. The results using other algorithms also show a similar trend.20 5 CONCLUSIONS AND FUTURE W ORK 

## 5 Conclusions and Future W ork Our experiments on tw o real gene expression data sets and four sets of synthetic data show that clustering with the PC’ s instead of the original variables does not necessarily impro ve cluster quality
. Our empirical study shows that the traditional wisdom that the first few PC’ s that contain most of the v ariation in the data may help to extract cluster structure is generally not true. W e also show that there usually exists some sets of  PC’ s that achieve higher quality of clustering results than the first  PC’ s. Our empirical results show that clustering with PC’ s has dif ferent impact on dif

ferent algorithms and dif ferent similarity metrics. When CAST is used with correlation as the similarity metric, clustering with the first  components usually gives a lower adjusted Rand index than clustering with the original v ariables (this is true in both of the real gene expression data sets and in 3 out of the 4 synthetic data sets). On the other hand, when k-means is used with correlation as the similarity metric, using all of the PC’ s in cluster analysis instead of the original variables gives higher or similar adjusted Rand indices on all of our real and synthetic data sets. When Euclidean distance is used as the similarity metric, clustering (either with k-means or a verage-link) using the first couple of PC’ s usually achieves higher or comparable adjusted Rand indices to without PCA, but the adjusted Rand indices drop sharply with more PC’ s. Since the Euclidean distance computed with the first  PC’ s is just an approximation to the Euclidean distance computed with all the experiments, the first couple of PC’ s probably contain most of the cluster information while the last PC’ s are mostly noise. There is no clear indication from our results what should be the number of PC’ s to use in the case of Euclidean distance. Using the number of first PC’ s chosen by the rule of thumb to cover 90% of the total variation in the data is too many in the case of Euclidean distance on the ovary data and yeast cell cycle data. Based on our empirical results, we recommend against using the first few PC’ s if CAST and correlation is used to cluster a gene expression data set. On the other hand, we recommend using all of the PC’ s if k-means and correlation is used instead. However, the increased adjusted Rand indices using the “appropriate” PC’ s with k-means and a verage-link are comparable to that of CAST using the original v ariables in many of our results. Therefore, choosing a good clustering algorithm is as important as choosing the “appropriate” PC’ s. There does not seem to be any general relationship between cluster quality ( i.e ., adjusted Rand index with an external standard) and the number of components used based on the results on both real and synthetic data sets. The choice of the first few components is usually not optimal (except when Euclidean distance is used), and usually may even achieve lower adjusted Rand indices than without PCA. There usually exists another set of PC’ s that achieves higher adjusted Rand indices than clustering with the original variables or with the first  PC’ s. However, there does not seem to be any general trend for the the set of components chosen by the greedy or modified greedy approach that has a high adjusted Rand index. Usually, there are no external criteria a vailable for real gene expression data, so it w ould be very useful if a rule to choose PC’ s for cluster analysis is a vailable. A careful manual inspection of our empirical results shows that the first tw o PC’ s are usually chosen in the exhaustive search step for the set of   components that give the highest adjusted Rand indices. In f act, when CAST is used with correlation as the similarity metric, the 3 components found in the exhaustive search step always include the first tw o PC’ s on all of our real and synthetic data sets. The first tw o PC’ s are usually returned by the exhaustive search step when k-means with correlation, or k-means with Euclidean distance, or a verage-link with correlation is used on all of our data sets except the synthetic data set with spik y clusters. On our real gene expression data sets, the first 2 PC’ s are alw ays returned by the exhaustive search step except when k-means with Euclidean distance is applied21

to the o v ary data. The patterns of results using the synthetic mixture of normal distributions and the randomly permuted o v ary data are very similar to those on the actual ovary data. This implies that these tw o synthetic data sets probably ha

ve similar complexities as the real gene expression data set. On the other hand, the patterns of results on both cyclic data sets are very dif ferent than those on the yeast cell cycle and those on the o v ary data, implying that the cyclic data sets may not be as satisf actory models for synthesizing gene expression data as the mixture of normal and the randomly permuted models. In particular, the cyclic data set with dif

ferent cluster sizes achieves very high (close to 1) adjusted Rand indices even without PCA, which is very dif ferent than what we observe on real data. There are a few possible directions of both empirical and theoretical future w ork. Empirically, it w ould be interesting to generate more replicates ( i.e ., more than 10) for each set of synthetic data to see if the standard deviations from the a verage pattern w ould go down. Furthermore, it w ould be interesting to check if increasing the parameter G (the number of best sets of components to keep in each search step) w ould significantly impro ve the adjusted Rand index computed. Comparing the adjusted Rand indices using a set of random PC’ s to those computed with the greedy approach w ould also be interesting. Our observ ation abo ve suggests that the first tw o PC’ s are usually chosen in the exhaustive search step. Generating a set of random PC’ s that alw ays includes the first tw o PC’ s, and then applying clustering algorithms and computing the adjusted Rand indices may also lead to interesting insights. In terms of future theoretical w ork, it is interesting to develop time-ef ficient approximation algorithms (other than the greedy and modified greedy approaches) to compute a set of PC’ s that achieves a high adjusted Rand index. Developing other models for generating synthetic gene expression data is definitely of interest. In addition, formally testing the normality of the ovary data w ould be useful to support the synthetic data model of the mixture of multivariate normal distributions. Our empirical study shows that the ef fectiveness of PCA on cluster analysis depends on the particular data set, the clustering algorithm and the similarity metric used. F or most real gene expression data sets, an external criterion to assess clustering results is not a vailable. It w ould be very valuable to develop a methodology that does not require an external criterion to evaluate the ef fectiveness of PCA as a preprocessing step. In our previous w ork [Y eung et al. , 2000 ], we proposed a methodology that estimates the predictive power of clustering algorithms. W e believe that our methodology in [Y eung et al. , 2000 ] can be modified to in vestigate the ef

fectiveness of clustering with the PC’ s instead of the original variables. It w ould be interesting to compare the evaluation results using our methodology with the results in this paper using external criteria. T o conclude, we believe that our empirical study is one step forw ard to in vestigate the ef

fectiveness of clustering with the PC’ s instead of the original variables. 

## Ackno wledgement I w ould like to thank Michel Schummer from the Department of Molecular Biotechnology at University of W

ashington for his unpublished ovary data set and his feedback. I w ould also like to thank Chris Frayley from the Statistics Department at University of W ashington for sharing with me the Splus code to generate a mixture of Gaussian variables. I w ould also like to thank Lue Ping Zhao at the Fred Hutchinson Cancer Research Center for his ideas of the cyclic models for gene expression data. I w ould also like to thank Amir Ben-Dor at Agilent Laboratories for sharing the additional heuristics implemented in B I O C L U S T with me. In addition, I w ould like to thank my friends Mathieu Blanchette, Jeremy Buhler and Don P atterson for their comments on an early draft of this paper. Finally, I w ould like to thank my advisor Larry Ruzzo, and my committee members Mark Campbell, Pedro Domingos, Phil Green, Da

vid Haynor and Martin T ompa for their feedback and suggestions during the examina-22 B EXAMPLE ILLUSTRA TING THE ADJUSTED RAND INDEX
tion and for the writeup. 

## Appendix A Scree graphs f

or the o v ary and y east cell cycle data

0 0
. 0 5 0
. 1 0
. 1 5 0
. 2 0
. 2 5 0
. 3 0
. 3 5

1 3 5 7 9 1 1 1 3 1 5 1 7 1 9 2 1 2 3  
		eigenvalue

Figure 21: Scree graph for the ovary data

0 0
. 0 5 0
. 1 0
. 1 5 0
. 2 0
. 2 5 0
. 3 0
. 3 5 0
. 4 0
. 4 5

1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 1 3 1 4 1 5 1 6 1 7 	

	eigenvalue

Figure 22: Scree graph for the yeast cell cycle data Figure 21 shows the scree graph for the ovary data. There is sharp change of steepness in Figure 21 at 3 components, and another gentle change at 6 components. Similarly, Figure 22 shows the scree graph for the yeast cell cycle data. Again, there is a sharp change of steepness at 3 components, and another gentle change at 5 components. These tw o examples illustrate that the scree graph approach to decide the number of first PC’ s to be used is very ad-hoc and subjective. 

## B Example illustrating the adjusted Rand index The following example illustrates how the adjusted Rand index (discussed in Section 3.3) is computed. Example 1 is a contingency table in the same form as in T able 1. W is defined as the number of pairs of objects in the same class in  and same cluster in , hence W can be written as :  D >

8

.

 > 

:

. In Example 1, W

8 8

 

:  8

 

: 8 
.  is defined as the number of pairs of objects in the same class in  but not in the same cluster in
. In terms of the notation in T able 1,  can be written as : 

8

.

 %$  ?:  :  D >

8

.

 >  ;: In Example 1, 

8 8

 

:  8

 

:  8

 

:  

8 
. Similarly, is defined as the number of pairs of objects in the same cluster in but not in the same23 Class or Cluster V

 V

 V Sums   1 1 0 2   1 2 1 4  0 0 4 4 Sums 2 3 5 .

8 _Rg Example 1

class in  , so can be written as :

>

8

.

$ > 

:



:  D >

8

.

 > 

: 8 8

 

:

 8

. 

:

 8

 

:

 

8 
.  is defined as the number of pairs of objects that are not in the same class in  and not in the same cluster in
. Since W

 

  

8 8

.  D:
, 

8 8

_Rg  :  

 

 

8 


. The Rand index for comparing the tw o partitions in Example 1 is

    8 g

Q 

_

_ , while the adjusted Rand index is

 

 



  



  



 

 





8 g

Q . X_ . (see Section 3.3 for the definitions of the Rand and adjusted Rand indices). The Rand index is much higher than the adjusted Rand index. 

## C Details of the CAST algorithm The Cluster Af

finity Search T

echnique (CAST) is an algorithm proposed by [Ben-Dor and Y akhini, 1999 ] to cluster gene expression data. The input to the algorithm includes the pairwise similarities of the genes, and a cutof

f parameter  (which is a real number between 0 and 1). The clusters are constructed one at a time. The current cluster under construction is called   < '  . The af finity of a gene  , W

M 

T , is defined to be the sum of similarity v alues between  and all the genes in  <

'  . A gene  is said to ha ve high af finity if W

M 

T   

 <

' 


. Otherwise,  is said to ha

ve low af finity
. Note that the af

finity of a gene depends on the genes that are already in  <

'  . The algorithm alternates between adding high af finity genes to  <

'  , and removing low af

finity genes from  <

'  .  <

'  is closed when no more genes can be added to or removed from it. Once a cluster is closed, it is not considered any more by the algorithm. The algorithm iterates until all the genes ha ve been assigned to clusters and the current   < '  is closed. When a new cluster  <

'  is started, the initial af finity of all genes are 0 since  <

'  is empty. One additional heuristic that the authors [Ben-Dor and Y akhini, 1999 ] implemented in their softw are B I O C L U S T is to choose a gene with the maximum number of neighbors to start a new cluster
. Another heuristic is that after the CAST algorithm con verges, there is an additional iterative step, in which all clusters are considered at the same time, and genes are moved to the cluster with the highest a verage similarity
. 

## D Correlation coefficient when there are 2 components When there are 2 components, the correlation coef ficient is either 1 or -1. Suppose there are tw o genes   and 

 with tw o components. Let  D > (where i

O 4 8 _

O  ) be the expression level of gene i under component 4
. The correlation coef

ficient between 

 and 

 can be simplified to: M 

D 

 

D 

T M



D 

 

D 

T  M 

D 

 

D 

T

 M



D 

 

D 

T

 (3)24 E VISU ALIZA TION OF THE CLUSTERING RESUL T WITH K-MEANS AND EUCLIDEAN DIST ANCE ON THE O

V AR Y D A T A
Since the denominator in Equation 3 represents the product of the norms of genes 

 and 

 , the denominator must be positive. From Equation 3, the correlation coef ficient between genes 

 and 

 is 1 if M 

D 

 

D 

T M



D 

 

D 

T g , the correlation coef

ficient is -1 if M 

D 

 

D 

T M



D 

 

D 

T

j g
. If 

D 

8 

D  or 

D 

8 

D  , the correlation coef

ficient is undefined. Since there are only tw o possible values that the correlation coef ficient can take when there are tw o components, there are at most tw o clusters. 

## E Visualization of the clustering result with k-means and Euclidean distance on the o v ary data The results on the ov ary data set using k-means and Euclidean distance in Figure 9 show that the adjusted Rand indices are high for first 2 and 3 PC’ s and then drop drastically to below that without PCA. When the four classes of the ovary data are viewed in the space of the first tw o PC’ s (Figure 23), the four classes are reasonably well-separated in the Euclidean space. In f act, when the clustering result using the first 3 PC’ s is viewed in the space of the first tw o PC’ s, the four clusters correspond mostly to the four classes in Figure 24. However, tw o classes (class 1 and class 4) overlap in Figure 25, in which the classes are viewed in the space of the second and fourth PC’ s. In f act, manual inspection shows that the tw o o verlapping classes were combined into one cluster by k-means with Euclidean distance when four PC’ s are used. The fourth PC probably is mostly noise, which makes it more dif ficult for k-means to extract the four classes.

- 3 - 2 - 1 0 1 2 3 4

- 6 - 4 - 2 0 2 4

   

    Figure 23: V

isualization of the four classes from the ovary data in the space of the first tw o PC’ s.

- 3 - 2 - 1 0 1 2 3 4

- 6 - 4 - 2 0 2 4

P C 1PC2

                                    Figure 24: V

isualization of the four clusters using the first 3 PC’ s and k-means with Euclidean distance from the ov ary data in the space of the first tw o PC’ s.REFERENCES 25

- 3 - 2 - 1 0 1 2 3 4 5

- 4 - 2 0 2 4

   

   Figure 25: V

isualization of the four classes from the ovary data in the space of the second and fourth PC’ s. 

## References [Anderber g, 1973] Anderber

g, M. R. (1973) Cluster analysis for applications
. Academic Press. [Ben-Dor and Y

akhini, 1999] Ben-Dor, A. and Y akhini, Z. (1999) Clustering gene expression patterns. In RECOMB99: Pr

oceedings of the Thir

d Annual International Confer ence on Computational Molecular Biology , L

yon, France. [Chang, 1983] Chang, W
. C. (1983) On using principal components before separating a mixture of tw o multivariate normal distributions. Applied Statistics , 32 , 267–275. [Chen et al. , 1999] Chen, T., Filko v, V
. and Skiena, S. S. (1999) Identifying gene re gulatory netw orks from experimental data. In RECOMB99: Pr oceedings of the Thir

d Annual International Confer

ence on Computational Molecular Biology , L yon, France. [Cho et al. , 1998] Cho, R. J., Campbell, M. J., W inzeler, E. A., Steinmetz, L., Conw ay, A., W odicka, L., W olfsber g, T
. G., Gabrielian, A. E., Landsman, D., Lockhart, D. J. and Da vis, R. W
. (1998) A genome-wide transcriptional analysis of the mitotic cell cycle. Molecular Cell , 2 , 65–73. [Chu et al. , 1998] Chu, S., DeRisi, J., Eisen, M., Mulholland, J., Botstein, D., Brown, P
. O. and Her

skowitz, I. (1998) The transcriptional program of sporulation in budding yeast. Science , 282 , 699–705. [Dunteman, 1989] Dunteman, G. H. (1989) Princiapl Components Analysis
. Sage Publications. [Eisen et al. , 1998] Eisen, M. B., Spellman, P
. T., Brown, P

. O. and Botstein, D. (1998) Cluster analysis and display of genome-wide expression patterns. Pr oceedings of the National Academy of Science USA , 95 , 14863–14868. [Everitt and Dunn, 1992] Everitt, B. S. and Dunn, G. (1992) Applied Multivariate Data Analysis
. Oxford University Press. [Hubert and Arabie, 1985] Hubert, L. and Arabie, P
. (1985) Comparing partitions. J ournal of Classification , 193–218. [Jain and Dubes, 1988] Jain, A. K. and Dubes, R. C. (1988) Algorithms for Clustering Data
. Prentice Hall, Englew ood Clif

fs, NJ.26 REFERENCES

[Jollif fe, 1986] Jollif

fe, I. T
. (1986) Principal Component Analysis
. New Y ork : Springer-V

erlag. [Jollif fe et al. , 1980] Jollif

fe, I. T., Jones, B. and t. Morgan, B. J. (1980) cluster analysis of the elderly at home: a case study
. Data analysis and Informatics , 745–757. [Lander, 1999] Lander, E. S. (1999) Array of hope. Natur e Genetics , 21 , 3–4. [Milligan and Cooper, 1986] Milligan, G. W
. and Cooper, M. C. (1986) A study of the comparability of external criteria for hierarchical cluster analysis. Multivariate Behavior al Resear ch , 21 , 441– 458. [Milligan et al. , 1983] Milligan, G. W ., Soon, S. C. and Sokol, L. M. (1983) The ef fect of cluster size, dimensionality, and the number of clusters on reco very of true cluster structure. IEEE T ransactions on P attern Analysis and Machine Intellig ence , 5 , 40–47. [Rand, 1971] Rand, W

. M. (1971) Objective criteria for the evaluation of clustering methods. J ournal of the American Statistical Association , 66 , 846–850. [Raychaudhuri et al. , 2000] Raychaudhuri, S., Stuart, J. M. and Altman, R. B. (2000) Principal components analysis to summarize microarray experiments: application to sporulation time series. In P acific Symposium on Biocomputing , vol. 5. [Schummer, 2000] Schummer, M. (2000) Manuscript in preparation. [Schummer et al. , 1999] Schummer, M., Ng, W
. V ., Bumgarner, R. E., Nelson, P
. S., Schummer, B., Bednarski, D. W

., Hassell, L., Baldwin, R. L., Karlan, B. Y
. and Hood, L. (1999) Compartive hybridization of an array of 21500 ovarian cdnas for the discovery of genes overexpressed in ovarian carcinomas. An International J ournal on Genes and Genomes , 238 , 375–385. [T amayo et al. , 1999] T

amayo, P ., Slonim, D., Mesirov, J., Zhu, Q., Kitareew an, S., Dmitro vsk y, E., Lander, E. S. and Golub, T
. R. (1999) Interpreting patterns of gene expression with selfor ganizing maps: methods and application to hematopoietic dif ferentiation. Pr

oceedings of the National Academy of Science USA , 96 , 2907–2912. [W en et al. , 1998] W

en, X., Fuhrman, S., Michaels, G. S., Carr, D. B., Smith, S., Barker, J. L. and Somogyi, R. (1998) Large-scale temporal gene expression mapping of central nervous system development. Pr

oceedings of the National Academy of Science USA , 95 , 334–339. [Y eung et al. , 2000] Y

eung, K. Y

., Haynor, D. R. and Ruzzo, W
. L. (2000) V

alidating clustering for gene expression data. T

ech. Rep. UW

-CSE-00-01-01, Dept. of Computer Science and Engineer

ing, University of W

ashington. [Zhao, 2000] Zhao, L. P

. (2000) Personal communications. [Zipf, 1949] Zipf, G. K. (1949) Human behavior and the principle of least ef fort
. Addison-W

esley.