Volume 10 Supplement 1

## Selected papers from the Seventh Asia-Pacific Bioinformatics Conference (APBC 2009)

- Research
- Open Access

# Principal component tests: applied to temporal gene expression data

- Wensheng Zhang
^{1}, - Hong-Bin Fang
^{2}and - Jiuzhou Song
^{1}Email author

**10 (Suppl 1)**:S26

https://doi.org/10.1186/1471-2105-10-S1-S26

© Zhang et al; licensee BioMed Central Ltd. 2009

**Published:**30 January 2009

## Abstract

### Background

Clustering analysis is a common statistical tool for knowledge discovery. It is mainly conducted when a project still is in the exploratory phase without any priori hypotheses. However, the statistical significance testing between the clusters can be meaningful in helping the researchers to assess if the classification results from implementing a clustering algorithm need to be improved, even after the cluster number has been determined by a well-established criterion. This is important when we want to identify highly-specific patterns through classification.

### Results

We proposed to use a principal component (PC) test, which is an implementation of an exact *F* statistic for the measures at multiple endpoints based on elliptical distribution theory, to assess the statistical significance between clusters. A challenge in the implementation is the choice of the number (q) of principal components to be considered, which can severely influence the statistical power of the method. We optimized the determination via validation according to a permutation test based on the clustering to be evaluated. The method was applied to a public dataset in classifying genes according to their temporal gene expression profiles.

### Conclusion

The results demonstrated that the PC testing were useful for determining the optimal number of clusters.

## Keywords

- Cluster Algorithm
- Cluster Number
- Agglomerative Hierarchical Cluster
- Functional Enrichment Analysis
- Silhouette Width

## Background

Data clustering is a common technique for statistical data analysis used in many fields [1], including machine learning, data mining, pattern recognition, and image analysis. Theoretically, clustering analysis identifies and classifies objects (or individuals) based on the similarity of the characteristics they possess. It seeks to minimize within-group variation and maximize between-group variation and results in a number of heterogeneous groups with homogeneous contents. The general categories of clustering methods include tree clustering (hierarchical clustering), block clustering, *k*-means clustering, and model-based clustering [1]. The evaluation of clustering analysis is a critical challenge in both theory and application.

The performance of clustering analysis can be assessed statistically in order to determine the appropriate clustering methods and cluster number [2]. Pseudo F statistic [3] is widely used for partitioning clustering algorithms, such as k-means, and has been included in the procedure FASTCLUST of SAS software [4]. BIC (Bayesian information criterion) is a well-established statistic based on standard statistical theory and fits model-based clustering procedures [5], which has been widely applied in bioinformatics [6–9]. Silhouette score [10] provides a measure of how well a data point was classified when it was assigned to a cluster according to both the tightness of the clusters and the separation between them. It has been used together with PAM (Partitioning Around Medoids) clustering algorithm [1]. Recently, the so-called Gap statistic was proposed [11], which can use the output of any clustering algorithms for the optimization of cluster number. Furthermore, clustering algorithms are commonly assessed from other angles, such as robustness, stability, consistency, and functional congruence of the members of the same cluster [2, 12–18].

On the other hand, while clustering analysis is mainly conducted when we are still in the exploratory phase of our research and do not have any prior hypotheses, the statistical significance testing between the clusters can be meaningful. The testing can help us to assess whether the classification results from running a clustering algorithm need to be improved, even after the cluster number has been determined by a well-established criterion. This is important in the clustering of genes on the basis of the temporal expression profiles. In order to extract specific knowledge about gene function from the expression profiles [19–21], researchers usually hope to have the number of clusters as large as possible but the contrasts between the clusters, each of which corresponds to a co-regulation pattern, should be statistically significant in general.

The significance testing between the clusters can be done by using Hotelling's *T*^{2}, the multivariate counterpart of Student's-*t* [22]. But when the number of measurement points is large and the size of samples is relatively small, the results from Hotelling's test are usually unstable [23]. Using the invariance of elliptical distribution theory, a type of exact *t* and *F* tests was proposed [23], which can be applied to high-dimension data with a small size of samples. The tests are based on the sum aggregates of original variables similar to O'Brien's method [24] but superior to the latter in maintaining the prescribed level of significance. Two direct implementations of the method are a one-fold principal component (PC) test corresponding to the exact *t* test and a multi-fold principal component test corresponding to the exact *F* test. The comparison of PC test and *T*^{2} clearly demonstrated the fact that the stabilizing effect of principal components and PC test made better use of the factor structure of the data of multiple end-points.

Microarray technology allows thousands of genes to be measured simultaneously on a single slide. Unsupervised learning on the basis of clustering analysis of microarray temporal gene expression data has been widely studied in order to discover classes of expression patterns and identify groups of genes that are regulated in a similar manner [7, 13, 19, 20]. However in literature the evaluation of clustering analysis was limited to the global assessment of clustering methods. In this paper, we proposed to use principal component tests based on the exact *F* test for multiple endpoint measures [23] to assess the significance of the contrasts between the gene clusters from different clustering algorithms and implemented it on a public data set. The testing can be conducted after the global evaluation for improving clustering analysis.

## Results

### Clustering and patterns

A major difference between CL1 and CL2 was that the two big clusters in the former were divided into two or more smaller classes in the latter. For example, the aggregate of the cluster1 and cluster2 in CL1 approximately corresponded to the aggregate of the cluster2, cluster7, cluster13, cluster14, cluster20 and cluster23 in CL2. In CL3, the major (62%) 483 gens were classified into the first cluster, which largely corresponded to the first biggest cluster in CL2. But the expression patterns of genes of these clusters were very different. Therefore, implementation of PAM with silhouette score criterion does not seem fit for the addressed dataset.

### Determination of q value

### Statistical evaluation of clustering

### Biological evaluation of clustering

Biologically functional enrichment analysis of the gene lists of eight clusters ^{a}

Cluster ID | Enriched PANTHER biological Processes | p-value |
---|---|---|

CL1-1 | lipid, fatty acid and steroid metabolism | |

(GO: 0006629) | 0.0003 | |

steroid metabolism (GO: 0008202) | 0.0006 | |

cholesterol metabolism (GO: 0008203) | 0.0253 | |

CL1-2 | blood clotting (GO: 0007596) | 0.0068 |

cell cycle (GO: 0007049) | 0.0079 | |

oncogenesis | 0.0101 | |

nucleoside, nucleotide and nucleic acid metabolism (GO: 0006139) | 0.0119 | |

embryogenesis (GO: 0009700) | 0.0435 | |

cell cycle control (GO: 0000074) | 0.0376 | |

mRNA transcription (GO: 0006366) | 0.0825 | |

cell proliferation and differentiation | 0.0843 | |

(GO: 0031054; GO: 0008283) cell structure and motility (GO: 0007010) | 0.0867 | |

immunity and defense (GO: 0006952) | 0.0953 | |

CL2-2 | cell cycle | 0.0325 |

CL2-7 | No | |

CL2-13 | No | |

CL2-14 | immunity and defense (GO: 0006952) | 0.0313 |

developmental processes (GO: 0007275) | 0.0432 | |

oncogenesis | 0.0926 | |

CL2-20 | nucleoside, nucleotide and nucleic acid metabolism (GO: 0006139) | 0.0796 |

CL2-23 | lipid, fatty acid and steroid metabolism (GO: 0006629) | 0.0084 |

steroid metabolism (GO: 0008202) | 0.0560 |

## Discussion

Clustering analysis is a widely used tool for knowledge discovery. Moreover, it is applied as a routine method in biology in the post-genomic era. The evaluation of clustering is a problem in its application. In this study, we compared the results of different clustering algorithms from a unique angle by testing the statistical significance of the contrasts between the clusters. In our knowledge, this paper is the first investigation of this kind. We used q-fold PC test which is an implementation of Lauter's exact *F* tests [23] for the measures of multiple endpoints. The method is superior to Hotelling's *T*^{2} [22] because of the stabilizing effects of the principal components, especially for the data with small sample size. This is important when we want to identify highly-specific patterns via clustering analysis.

The significance of the proposed clustering evaluation includes three aspects. Firstly, the results can tell us if the clustering is meaningful, at least from a statistical standpoint. A good clustering algorithm should meet a basic criterion, i.e., the clusters should be statistically distinguishable. In other words, all of the contrasts between the clusters should be statistically significant at a certain confidence level. Second, it can be helpful in the determination of cluster numbers. For example, in the analysis of temporal gene expression data mentioned above, both the BIC plot did not have the expected "U" shape. Thus, the determination based on a local minimum value may be equivocal and questionable. The results of the PC tests demonstrated that dividing the 483 genes (probes) into 18–20 clusters is appropriate. Finally, the method is extremely useful for the improvement of the results from a clustering analysis by demonstrating which clusters can be combined because of the lack of significant difference between them.

The number (q) of principal components to be considered is a challenge for the PC test. We optimized the determination via validation according to permutation test based on the clustering to be evaluated. In this way, the choice of q is determined by the data and clustering methods. It is superior to the choice based on cumulative energy content (CEC) because the latter needs an artificial threshold of the CEC percentage. More importantly, from the permutation test, we can assess the validity of the PC test itself in controlling type I error.

An alternative approach to the evaluation of clustering of genes based on the temporal expression profiling is biological validation. In this paper, we conducted biologically functional enrichment analysis of the gene lists of several clusters of interest. The results showed that the finer division of clusters from SSClust, a model-based clustering algorithm, can provide more specific relationships between clusters and biological functions.

It is worthy to note that the information from the biological validation is usually limited because the temporal gene expression profiles of the genes involved in a biological process can be very diverse, including, for instance, inverse co-regulation or co-regulation with a time lag or a combination of both [21, 27].

## Conclusion

The proposed PCA test method was applied to a public dataset in classifying genes according to their temporal gene expression profiles. The results demonstrated that the PC testing were useful for determining the optimal number of clusters. We also anticipate that the method could be used for pattern identification and similarity analysis.

## Methods

### Data

The initial data set, published by Iyer et al. [25], describes the transcription levels of genes detected by 517 gene probes, corresponding to 497 unique genes, during the first 24 h of the serum response in serum-starved human fibroblasts. By using an agglomerative hierarchical clustering method, the authors [25] detected 10 major gene expression profile clusters among the differentially expressed genes of the serum response. The ten classes contained 465 unique genes or 483 gene probes. Our work was focused on the data of these 483 gene probes with the log-transformed expression ratios as the variables. The gene symbols of 239 annotated genes were provided by Lagreid et al (2006).

### Principal component tests

The q-fold principal test (PC) used in this paper is implemented on the basis of a type of *t* or *F* statistic for high-dimension data.

Assume there are n individuals (genes) and from each one we have p observations at different time points. Assume p-dimensional distribution for x_{
i
}(i = 1, 2, n), i.e. x_{
i
}~ N(**μ**_{
i
}, **Σ**). Denote $X=({{x}^{\prime}}_{1},{{x}^{\prime}}_{2}\cdots {{x}^{\prime}}_{n}{)}^{\prime}$, a p × n matrix representing the gene expression. We have

X~ N_{p× n}(M, **Σ** ⊗ I_{
n
}),

where $M=({{\mu}^{\prime}}_{1},{{\mu}^{\prime}}_{2}\cdots {{\mu}^{\prime}}_{n}{)}^{\prime}$, **Σ** is a variance and covariance matrix.

**μ**

_{ 1 }=

**μ**

_{ 2 }= ...

**μ**

_{ n }, i.e.

^{(1)}and n

^{(2)}represent the numbers of genes in the two populations (clusters), respectively, vector k can be calculated with following equation,

**Λ**is the q × q diagonal matrix of q largest eigenvalues, then, Z= D'X has a matrix elliptical contoured distribution [28]. Based on the invariance of elliptically contoured distributions, if H

_{0}holds, the statistic

exactly follows F distribution with q and n-q-1 as the degrees of freedom [23],

where $G=Z\left({I}_{n}-\frac{1}{n}{1}_{n}{{1}^{\prime}}_{n}-k{k}^{\prime}\right)Z$. For a given n and p, the power of this statistic is dependent on the choice of q. When q = 1, the statistic (5) has t-distribution with degree of freedom n-2.

### Determination of q value

The number (q) of principal components to be considered is a challenge for the q-fold PC test. A solution is the choice based on cumulative energy content (CEC). However, the threshold of the CEC percentage has to be artificially determined. Here, we developed a permutation test based on the clustering to be evaluated. Let I_{c} be a vector containing the cluster IDs of the genes in the clustering. By shuffling, we get another vector ${\text{I}}_{\text{c}}^{\ast}$ which has all elements of I_{c} arranged in a random order. We, then, replace I_{c} with ${\text{I}}_{\text{c}}^{\ast}$ and carry our significance testing on the $\text{M}=\frac{1}{2}\text{k}\times (\text{k}-1)$ contrasts (k is the cluster numbers) between clusters using PC test with different q (q = 1, 2,...). For each q, we count the number (m) of the random contrasts with p-value smaller than the prescribed error of type-I at α, such as 0.05, and calculate the ratio R = m/M. Finally, we chose the minimum q which, R approximately equals to α. If the cluster number is small, the shuffling-testing procedure should be repeated several times.

### Clustering methods

The results from three clustering algorithms were evaluated in this paper. Following is a simple description of these methods.

#### Agglomerative hierarchical clustering

An agglomerative hierarchical clustering procedure produces a series of partitions of the data, P_{n}, P_{n-1},..., P_{1}, the first P_{n} consisting of n single object "clusters", the last P_{1}, consisting of a single group containing all n cases. At each stage the method joins together two clusters which are closest together (most similar) [19]. Differences between methods in this category arise because of the different ways of defining distance (or similarity) between clusters.

#### Model-based clustering with smoothing splines (SSClust)

where L is the likelihood for the mixture model, N is total gene number, k is the cluster number, and v_{i} is the numbers of free parameters for i^{th} cluster which is equivalent to the sum of the trace of the smoothing matrix [30]. A small BIC score indicates strong evidence for the corresponding clustering.

#### Partitioning Around Medoids (PAM)

*k*-means algorithm. It operates on the dissimilarity matrix of the given data set [1]. Compared with the ordinary

*k*-means, PAM is more robust, because it minimizes a sum of dissimilarities instead of a sum of squared Euclidean distances. PAM first computes

*k*representative objects, called medoids. A medoid can be defined as a characteristic a cluster, whose average dissimilarity to all the objects in the cluster is minimal. After finding the set of medoids, each object of the data set is assigned to the nearest medoid. That is, object

*i*is put into cluster

*v*

_{ i }, when medoid

*mv*

_{ i }is nearer than any other medoid

*m*

_{ w }. We used the

*pam*program in R package "

*cluster"*in Bioconductor, where the optimal number of clusters is selected on the silhouette plot. Silhouette score [10] is obtained by taking the mean of the average silhouette width for all clusters and silhouette width is defined as

where *a*(*i*) is the average distance of gene *i* to other genes in the same cluster, *b*(*i*) is the average distance of gene *i* to genes in its nearest neighboring cluster. Like BIC, a small silhouette score indicates evidence for the corresponding clustering.

### Functional enrichment analysis

A web-tool in PANTHER classification system [26] was used for the biologically functional enrichment analysis by comparing the lists of member genes contained in each cluster of interest with gene from H. Sapiens in NCBI. Only PANTHER biological processes, most of which can be exactly mapped to a Gene Ontology (GO) term [31], were investigated at detail. The p-values were firstly calculated on the basis of hyper-geometric distribution theory followed by correction for multiple testing using the Bonferroni method. Because the correction method is conservative, in the following text a biological process with adjusted p-value < 0.1 was considered as "significant".

## Declarations

### Acknowledgements

The authors are grateful to the three reviewers for their constructive comments. This study was supported by a USDA-NRI grant (NRI Proposal No.2008-00870).

This article has been published as part of *BMC Bioinformatics* Volume 10 Supplement 1, 2009: Proceedings of The Seventh Asia Pacific Bioinformatics Conference (APBC) 2009. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/10?issue=S1

## Authors’ Affiliations

## References

- Kaufman L, Rousseeuw P: Finding Groups in Data: An Introduction to Cluster Analysis. 1990, Wiley, New YorkView ArticleGoogle Scholar
- Milligan GW, Cooper MC: An examination of procedures for determining number of clusters in a data set. Psychometrika. 1985, 50: 159-179. 10.1007/BF02294245.View ArticleGoogle Scholar
- Calinski T, Harabasz J: A dendrite method for cluster analysis. Commun Statist. 1974, 3: 1-27. 10.1080/03610927408827101.View ArticleGoogle Scholar
- Institute S: SAS/STAT User's Guider. 2002Google Scholar
- Fraley C, Raftery AE: Model-Based Clustering, Discriminant Analysis, and Density Estimation. Journal of the American Statistical Association. 2002, 97: 611-631. 10.1198/016214502760047131.View ArticleGoogle Scholar
- Luan Y, Li H: Clustering of time-course gene expression data using a mixed-effects model with B-splines. Bioinformatics. 2003, 19 (4): 474-482. 10.1093/bioinformatics/btg014.View ArticlePubMedGoogle Scholar
- Ma P, Castillo-Davis CI, Zhong W, Liu JS: A data-driven clustering method for time course gene expression data. Nucleic Acids Res. 2006, 34 (4): 1261-1269. 10.1093/nar/gkl013.PubMed CentralView ArticlePubMedGoogle Scholar
- Yeung KY, Fraley C, Murua A, Raftery AE, Ruzzo WL: Model-based clustering and data transformations for gene expression data. Bioinformatics. 2001, 17 (10): 977-987. 10.1093/bioinformatics/17.10.977.View ArticlePubMedGoogle Scholar
- McLachlan GJ, Bean RW, Peel D: A mixture model-based approach to the clustering of microarray expression data. Bioinformatics. 2002, 18 (3): 413-422. 10.1093/bioinformatics/18.3.413.View ArticlePubMedGoogle Scholar
- Rousseeuw PJ: Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics. 1987, 20: 53-65. 10.1016/0377-0427(87)90125-7.View ArticleGoogle Scholar
- Tibshirani R, Walther G, Hastie T: Estimating the number of clusters in a dataset via the Gap statistic. Journal of the Royal Statistical Society B. 2001, 63: 411-423. 10.1111/1467-9868.00293.View ArticleGoogle Scholar
- Smolkin M, Ghosh D: Cluster stability scores for microarray data in cancer studies. BMC Bioinformatics. 2003, 4: 36-10.1186/1471-2105-4-36.PubMed CentralView ArticlePubMedGoogle Scholar
- Chen X, Jaradat SA, Banerjee N, Tanaka TS, Ko MSH, Zhang MQ: Evaluation and comparison of clustering algorithms in analyzing ES cell gene expression data. Statistica Sinica. 2002, 12: 241-262.Google Scholar
- Datta S, Datta S: Evaluation of clustering algorithms for gene expression data. BMC Bioinformatics. 2006, 7 (Suppl 4): S17-10.1186/1471-2105-7-S4-S17.PubMed CentralView ArticlePubMedGoogle Scholar
- Raatikainen KEE: Cluster analysis and workload classification. Performance Evaluation Review. 1993, 20 (4): 24-30. 10.1145/155775.155781.View ArticleGoogle Scholar
- Thalamuthu A, Mukhopadhyay I, Zheng X, Tseng GC: Evaluation and comparison of gene clustering methods in microarray analysis. Bioinformatics. 2006, 22 (19): 2405-2412. 10.1093/bioinformatics/btl406.View ArticlePubMedGoogle Scholar
- Rand WM: Objective criteria for the evaluation of clustering methods. Journal of the American Statistical Association. 1971, 66: 846-856. 10.2307/2284239.View ArticleGoogle Scholar
- Yeung KY, Haynor DR, Ruzzo WL: Validating clustering for gene expression data. Bioinformatics. 2001, 17 (4): 309-318. 10.1093/bioinformatics/17.4.309.View ArticlePubMedGoogle Scholar
- Eisen MB, Spellman PT, Brown PO, Botstein D: Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci USA. 1998, 95 (25): 14863-14868. 10.1073/pnas.95.25.14863.PubMed CentralView ArticlePubMedGoogle Scholar
- Wu LF, Hughes TR, Davierwala AP, Robinson MD, Stoughton R, Altschuler SJ: Large-scale prediction of Saccharomyces cerevisiae gene function using overlapping transcriptional clusters. Nat Genet. 2002, 31 (3): 255-265. 10.1038/ng906.View ArticlePubMedGoogle Scholar
- Lagreid A, Hvidsten TR, Midelfart H, Komorowski J, Sandvik AK: Predicting gene ontology biological process from temporal gene expression patterns. Genome Res. 2003, 13 (5): 965-979. 10.1101/gr.1144503.View ArticlePubMedGoogle Scholar
- Hotelling H: The generalization of Student's ratio. Ann Math Statist. 1931, 2: 360-378. 10.1214/aoms/1177732979.View ArticleGoogle Scholar
- Lauter J: Exact t and F tests for analyzing studies with multiple endpoints. Biometrics. 1995, 52: 964-970. 10.2307/2533057.View ArticleGoogle Scholar
- O'Brien PC: Procedures for comparing samples with multiple endpoints. Biometrics. 1985, 40: 1079-1087. 10.2307/2531158.View ArticleGoogle Scholar
- Iyer VR, Eisen MB, Ross DT, Schuler G, Moore T, Lee JC, Trent JM, Staudt LM, Hudson J, Boguski MS: The transcriptional program in the response of human fibroblasts to serum. Science. 1999, 283 (5398): 83-87. 10.1126/science.283.5398.83.View ArticlePubMedGoogle Scholar
- Thomas PD, Campbell MJ, Kejariwal A, Mi H, Karlak B, Daverman R, Diemer K, Muruganujan A, Narechania A: PANTHER: a library of protein families and subfamilies indexed by function. Genome Res. 2003, 13 (9): 2129-2141. 10.1101/gr.772403.PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang W, Song JZ: Term-tissue specific models for prediction of gene ontology biological processes using transcriptional profiles of aging in D. Melanogaster. 2007Google Scholar
- Fang K-T, Zhang J: Generalized multivariate analysis. 1990, Berlin, Heidelberg; Science Press Beijing and Springer-VerlagGoogle Scholar
- Rice JA, Wu CO: Nonparametric mixed effects models for unequally sampled noisy curves. Biometrics. 2001, 57 (1): 253-259. 10.1111/j.0006-341X.2001.00253.x.View ArticlePubMedGoogle Scholar
- Gu C: Smoothing Spline ANOVA Models. 2000, Springer-VerlagGoogle Scholar
- Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000, 25 (1): 25-29. 10.1038/75556.PubMed CentralView ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.