Conditional clustering of temporal expression profiles
© Wang et al; licensee BioMed Central Ltd. 2008
Received: 01 November 2007
Accepted: 11 March 2008
Published: 11 March 2008
Many microarray experiments produce temporal profiles in different biological conditions but common cluster techniques are not able to analyze the data conditional on the biological conditions.
This article presents a novel technique to cluster data from time course microarray experiments performed across several experimental conditions. Our algorithm uses polynomial models to describe the gene expression patterns over time, a full Bayesian approach with proper conjugate priors to make the algorithm invariant to linear transformations, and an iterative procedure to identify genes that have a common temporal expression profile across two or more experimental conditions, and genes that have a unique temporal profile in a specific condition.
We use simulated data to evaluate the effectiveness of this new algorithm in finding the correct number of clusters and in identifying genes with common and unique profiles. We also use the algorithm to characterize the response of human T cells to stimulations of antigen-receptor signaling gene expression temporal profiles measured in six different biological conditions and we identify common and unique genes. These studies suggest that the methodology proposed here is useful in identifying and distinguishing uniquely stimulated genes from commonly stimulated genes in response to variable stimuli. Software for using this clustering method is available from the project home page.
Cluster analysis of gene expression data is commonly used to group gene expression measurements, cross-sectionally or longitudinally, into categories of genes that have similar patterns of expression. Various clustering methods have been proposed to analyze microarray data [1–4], and hierarchical and k-means clustering have been applied to the discovery and characterization of the regulatory mechanisms of several processes and organisms [5–10]. Time course microarray experiments allow investigators to look at the behaviors of genes over time and hence introduce another dimension of observation . In the past few years, several methods for clustering temporal gene expression data have been proposed that use autoregressive models (CAGED) , hidden Markov models , mixture models with variations of the EM algorithm  or B splines , and Bayesian hierarchical models with an agglomerative search using smoothing splines . However, it has been suggested that autoregressive models  and some other clustering procedures [13, 17] are better suited to cluster long temporal gene expression data, possibly measured at regularly spaced time points . We introduced an extension to CAGED that overcomes this limitation and is more suitable to cluster short gene expression profiles in . This method uses polynomial models to describe the expression profiles and uses proper prior distributions for the model parameters to make the result of clustering invariant under linear transformations of time .
Furthermore, the algorithm finds the optimal number of clusters during the clustering process using a Bayesian decision-theoretic approach.
The common objective of all these different methods is to cluster gene expression temporal profiles observed in one specific biological condition, but in many experiments the temporal expression profiles are observed under different biological conditions: for example Diehn and coauthors  measured the gene expression profile of T cells under various activating stimulations to identify those genes that react uniquely to specific activating stimulations. To account for different biological conditions, Storey et al.  introduced an F statistic for differential analysis of time course expression data that produced a principled way to find the genes that are expressed differently in two experimental conditions. Here, the objective is different: we are not only interested in discovering the genes that have different dynamics in different experimental conditions. We wish to be able to simultaneously discover these genes and also group them together with other genes that have the same temporal expression profiles.
Other authors have already proposed solutions to this problem. Heard et al.  have proposed to extend their model-based clustering of temporal expression profiles in multiple conditions by modeling the concatenated time series. This approach has the advantage of increasing the robustness of clustering by simultaneously using the information from the different experiments. However, this gain of robustness comes at the price that only concatenated time series are clustered, so those genes that have a similar profile in one experimental condition but different profiles in the other conditions will be lost because they are allocated to different clusters. Ng et al  proposed a mixture model with random effect that uses the EM algorithm to cluster gene expression data from either time course experiment or experiments with replicates. The flexible nature of linear mixed model gives a unified approach to cluster data collected from various designed experiments and, in principle, this method is applicable to clustering temporal data measured in different conditions. However, similarly to the method in , it only clusters the concatenated time series from multiple conditions.
Here we propose an extension of our polynomial-based method to cluster short expression profiles measured in different conditions. This method that we call conditional clustering stratifies the data according to the experimental conditions and performs separate cluster analysis within the strata, then attempts to merge the resulting clusters if the merging could improve a Bayesian metrics. Because clustering within each stratum does not use all the information for those genes that have common patterns across more than one experimental condition, we further use an iterative procedure to find those genes that have unique expression profile under a specific condition and genes that have common expression profiles under two, or even more conditions. An application of this novel procedure to real data produces biologically meaningful gene sets.
Results and discussion
Cluster analysis within conditions
In matrix notation we can write
x g = Fβ g + ε g (1)
where , F is the n × (p + 1) design matrix in which the i th row is , β g = (μ g , βg 1, ..., β gp ) T is the vector of regression coefficients, and is the vector of errors. We assume that the errors are normally distributed, with and , for any t i , and p is the polynomial order. To complete the specification of the Bayesian model, we assume the standard conjugate normal-gamma prior density for the parameters β g and τ g so that the marginal distribution of τ g and the distribution of β g , conditional on τ g , are
τ g ~ Gamma(α1, α2)
β g |τ g ~ N(β0, (τ g R0)-1)
where E(τ g ) = α1α2, E(1/τ g ) = 1/((α1 -1)α2), and E(β g ) = β0 = (β00, β01, ..., β0p) T , and R0 is the (p + 1) × (p + 1) identity matrix. The prior hyper-parameters α1, α2, β0 are identical for all genes. We use the data of the non-expressed genes to specify the prior hyper-parameters, as suggested in [19, 25]. To compare different partitions of the genes, we compute the posterior probability of different clustering models so that, given the observed gene expression profiles, the best clustering model is the one with maximum posterior probability. This method was originally suggested in  and works as follows. Let M c denote the model with c clusters of gene expression data, where each cluster C k groups the set of expression profiles generated by the same polynomial models with coefficients β k and variance 1/τ k , k = 1, ..., c. Each cluster contains m k genes that are jointly modeled as
x k = F k β k + ε k
Note that we now label the vectors x g assigned to the same cluster C k with the double script k g , in which k denotes the cluster membership, and g = 1, ..., m k is the index for the genes in this cluster. Here m = ∑ k m k , ε k is the vector of uncorrelated errors with E(ε k ) = 0 and V(ε k ) = 1/τ k .
When all clustering models are a priori equally likely, the posterior probability p(M c |x) is proportional to the marginal likelihood f(x|M c ) that becomes our probabilistic scoring metric.
To make the computation feasible, the same agglomerative, finite-horizon heuristic search strategy introduced by Ramoni et al in  is used. This heuristic search orders the merging of profiles by their distance, so that the closest profiles are tested for merging first. The procedure first calculates the m(m - 1) pair-wise distances for the m genes and then attempts to merge the two closest genes into one cluster. If this merging increases the likelihood it is accepted, the two genes are assigned to the same clusters and their profile is replaced by the average profile that is a point-by-point average of the expression of the two genes. The same procedure is repeated for the new set of m - 1 profiles. If this merging is not accepted, the heuristic tries to merge the next second closest genes to see if their merging increases the likelihood or not. This procedure continues until an acceptable merging is found, otherwise it stops. This strategy automatically determines the best number of clusters when it does not find a pair of profiles to be merged into the same cluster and therefore stops. Note that the decision to merge profiles is based on the posterior probability and the distance between profiles is only used to speed up computations. In the implementation we have two possible choices of distances: the Euclidean distance and the negative of the correlation coefficient.
To find the best polynomial order, we repeat the conditional cluster analysis for different values p, and compute the Bayesian Information Criterion (BIC)  of the cluster set for each p and the value that optimizes the BIC is chosen. Let q be the number of parameters in the model, then q = c(p + 2) when the model does not have intercept and q = c(p + 3) when the model has intercept. The BIC of the clustering model M c is
BIC = -2 log f(x|M c ) + q log(n) (5)
where f(x|M c ) is the marginal likelihood of the model specified in Equation 4. We use the same p for clustering in each condition and the final merging of all conditions to ensure consistency of the overall clustering model.
We evaluated this clustering algorithm in three simulated data, and we demonstrated that this algorithm performs well in identifying the correct number of clusters as well as the correct cluster assignments . These results are consistent with other evaluations, where we showed that the combination of a distance driven search with a Bayesian scoring metrics leads to correctly identifying clusters in an efficient way [26, 27]. The results were compared to a competing program STEM  which is also designed for clustering short temporal gene expression data.
Cluster analysis across all experimental conditions
Cluster the gene expression profiles within each experimental condition.
Use all the clusters generated in step 1 as input of a new cluster analysis across all conditions using the same algorithm. During this step we try the merging of these initial clusters using the marginal likelihood in Equation 4 as scoring metric until no merging can improve the likelihood.
Identify the GCP and GUP in the clusters derived in step 2. If there are no GCP, so there are no clusters merging the profiles of the same gene in two or more conditions, go to step 4. If there are GCP go to step 5.
There are no GCP. The resulting clusters are the final clusters for GUP.
Remove those GCP from the analysis, but keep their cluster assignments. Re-cluster the remaining data containing GUP using our clustering algorithm.
Identify the GCP and GUP in the clusters produced at step 5. If there are no GCP, go to step 7. If there are GCP go to step 5.
Take the clusters containing the GCP removed during the iterations, try to merge these clusters of GCP to see if the merging could improve the marginal likelihood. The resulting clusters are the final clusters for GCP.
The results from this procedure give us two sets of clusters: clusters of GCP and clusters of GUP. The GCP clusters contain expression profiles for genes that have the same behavior in at least two experimental conditions. The GUP clusters have expression profiles of genes that have unique expression patterns in a particular experimental condition. We should notice that there are multiple conditions here and to choose the optimal polynomial order p we proceed as follows. First, for a fixed p, we perform the clustering in each condition, then we try to merge the clusters from all the conditions if that improves the likelihood, as mentioned above, and this is the initial results. The optimal p is selected by comparing BIC of the various initial results using p = 1, ..., n - 1.
We conducted three simulation studies to evaluate the performance of the proposed Bayesian conditional clustering algorithm. The first two simulations examine the effects of sample size and variability on the accuracy of the cluster produced as initial results. The third simulation examines the effectiveness of the whole clustering algorithm in finding GCP and GUP. In all three simulations we generated normalized patterns.
We clustered the simulated data using the conditional clustering algorithm and used a statistic proposed by Rand  to evaluate the similarity between the results and the true cluster assignment. The rationale of this statistic is that, given two sets of clusters, the pairs of objects that are either assigned to the same cluster or split across clusters in both sets show similar cluster assignments, so the statistic is simply the proportion of pairs of objects that are assigned consistently in the two sets. We borrow the example used in  to describe this statistic more in details. Suppose we wish to cluster six objects, say a, b, c, d, e, f, and consider two ways of grouping them: group 1 consists of the two clusters (a, b, c), (d, e, f), and group 2 consists of the three clusters (a, b), (c, d, e), and (f). The six objects can be grouped into 15 possible pairs and we can label the elements of each of the 15 pairs as either "assigned to the same cluster" in both groups, "assigned to different clusters" in both groups, or mixed in the other cases, based on the cluster assignments in the two groups. For example, the pair ab is assigned to the same cluster in both groups, the pair ac is mixed because it is assigned to the same cluster in the first group but a and c are assigned to different clusters in the second group, and the pair ad is split across two clusters in both groups. In the two groups above the two pairs ab, de are assigned to the same clusters, and the seven pairs ad, ae, af, bd, be, bf, cf are split across clusters in both groups so that Rand statistic is 9/15 = 0.6.
Cross classification table of the iterative clustering results for the simulated 700 genes.
Our Results Assignment
We applied our new clustering algorithm to the gene expression data from . The experiment used cDNA microarrays to study the genomic expressions of human T cells in an experimental control condition, and in response to stimulations of CD3, CD28, their co-stimulation CD3/CD28, lectin phytohemagglutinin (PHA), and a combination of the calcium ionophore ionomycin and the phorbol ester phorbol 12-myristate 13-acetate (PMA/lo). In each of these conditions, the expression profile of human T cells was observed at 0, 1, 2, 6, 12, 24, 48 hours. We used the expression profiles analyzed in , but removed the profiles with missing data so that we analyzed the expression levels of 2,362 genes.
Biological categories identified by EASE. NS = not significant
Category analysis with EASE
Total genes in clusters
Common genes in clusters
Specific non-overlapping genes for each stimulation
TF and Immune-Hs
Gene Ontology categories
GO CATEGORY: TFS AND IMMUNE
Categories of overlapping genes in ALL stimulations: DNA synthesis, cell cycle
Categories of genes in ALL stimulations: DNA binding, TFs and immune, TF binding, DNA synthesis, cell cycle
CEBPG; CHUK; CXCR4; HIF1A; HLA-E; ICAM1; IL8; IRAK1; IRF2; NFATC1; NFKB2; POU2F2; RELB; STAT1; STAT5B; TBX21; TCF7; VDR; VIPR1; WT1
Categories of CD3/CD28 specific genes (overlapping genes were subtracted): TF and immune, DNA binding, TF binding, chromatin, DNA synthesis, DNA binding, RNA processing, cell cycle
CEBPG; CHUK; HIF1A; ICAM1; IL8; IRAK1; NFATC1; NFKB2; POU2F2; RELB; STAT1; STAT5B; TBX21; VDR; WT1
Temporal microarray experiments across several experimental conditions are frequently conducted, and it is very important to recognize the specific features built in this design. In this paper we introduce a clustering technique that incorporates the experimental conditions in the analysis. More importantly, we develop an iterative procedure to distinguish, from a set of resulting clusters, clusters of common genes and unique genes. Our simulation study shows that this clustering algorithm can correctly identify the number of clusters, and find the common and unique genes with low error rate. The application of this technique to the analysis of data from  gives us a set of very interesting clusters biologically, and a new perspective to look at this data.
There are some limitations in this work that can stimulate further research. Other methods, for example bi-clustering, can be applied to cluster data measured in two conditions. Bi-clustering has the advantage of simultaneously clustering rows and columns, and it has been applied to co-clustering of genes and experimental conditions to find differentially expressed genes . However, the current versions of bi-clustering do not take into account the temporal nature of the data. Nevertheless, it would be interesting to compare the approaches. Another limitation of this work is that we considered the situation in which the sampling frequency is the same across experimental conditions. It is not straightforward to generalize our approach to experiments comparing expression profiles measured with different sampling frequency. One possibility is to make the sampling frequencies homogeneous by adding the time points that are missing and use, for example, imputations to fill in the missing observations or Markov Chain Monte Carlo methods to estimate the Bayesian score.
Availability and requirements
Project name: Bayesian conditional clustering of temporal expression profiles.
Project home page: http://people.bu.edu/sebas/condclust/index.htm.
Operating system: Microsoft Windows XP and Vista.
Programming language: R
This research was supported by NIH grant NHGRI R01 HG003354-01A2. The authors thank Al Ozonoff, Gheorghe Doros and Mike LaValley for their suggestions to improve and evaluate the algorithm of conditional clustering.
- Alter O, Brown PO, Botstein D: Singular value decomposition for genome-wide expression data processing and modeling. Proc Natl Acad Sci USA 2000, 97: 10101–10106. 10.1073/pnas.97.18.10101PubMed CentralView ArticlePubMedGoogle Scholar
- Butte AJ, Tamayo P, Slonim D, Golub TR, Kohane IS: Discovering functional relationships between rna expression and chemotherapeutic susceptibility using relevance networks. Proc Natl Acad Sci USA 2000, 97: 12182–12186. 10.1073/pnas.220392197PubMed CentralView ArticlePubMedGoogle Scholar
- Eisen M, Spellman P, Brown P, Botstein D: Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci USA 1998, 95: 14863–14868. 10.1073/pnas.95.25.14863PubMed CentralView ArticlePubMedGoogle Scholar
- Tamayo P, Slonim D, Mesirov J, Zhu Q, Kitareewan S, Dmitrovsky E, Lander ES, Golub TR: Interpreting patterns of gene expression with self-organizing maps: Methods and application to hematopoietic differentiation. Proc Natl Acad Sci USA 1999, 96: 2907–2912. 10.1073/pnas.96.6.2907PubMed CentralView ArticlePubMedGoogle Scholar
- Golub R, Slonim DK, Tamayo P, Huard C, Gaasenbeek M, Mesirov JP, Loh ML, Coller H, Downing JR, Caligiuri MA, Bloomfield CD, Lander ES: Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring. Science 1999, 286: 531–537. 10.1126/science.286.5439.531View ArticlePubMedGoogle Scholar
- Iyer VR, Eisen MB, Ross DT, Schuler G, Moore T, Lee JCF, Trent JM, Staudt LM, Hudson J Jr, Boguski MB, Lashkari D, Shalon D, Botstein D, Brown PO: The transcriptional program in the response of human fibroblasts to serum. Science 1999, 283(5398):83–87. 10.1126/science.283.5398.83View ArticlePubMedGoogle Scholar
- Lapointe J, Li C, Higgins JP, Van de Rijn M, Bair E, Montgomery K, Ferrari M, Egevad L, Rayford W, Bergerheim U, Ekman P, DeMarzo AM, Tibshirani R, Botstein D, Brown PO, Brooks JD, Pollack JR: Gene expression profiling identifies clinically relevant subtypes of prostate cancer. Proc Natl Acad Sci 2004, 101(3):811–816. 10.1073/pnas.0304146101PubMed CentralView ArticlePubMedGoogle Scholar
- Sorlie T, Tibshirani R, Parker J, Hastie T, Marron JS, Nobel A, Deng S, Johnsen H, Pesich R, Geisler S, Demeter J, Perou CM, Lonning PE, Brown PO, Borresen-Dale A, Botstein D: Repeated observation of breast tumor subtypes in independent gene expression data sets. Proc Natl Acad Sci 2003, 100(14):8418–8423. 10.1073/pnas.0932692100PubMed CentralView ArticlePubMedGoogle Scholar
- Spielman RS, Bastone LA, Burdick JT, Morley M, Ewens WJ, Cheung VG: Common genetic variants account for differences in gene expression among ethnic groups. 2003, 39(2):226–231.Google Scholar
- Wen X, Fuhrman S, Michaels GS, Carr DB, Smith S, Barker JL, Somogyi R: Large-scale temporal gene expression mapping of central nervous system development. Proc Natl Acad Sci USA 1998, 95: 334–339. 10.1073/pnas.95.1.334PubMed CentralView ArticlePubMedGoogle Scholar
- Bar-Joseph Z: Analyzing time series gene expression data. Bioinformatics 2004, 20(16):2493–503. 10.1093/bioinformatics/bth283View ArticlePubMedGoogle Scholar
- Ramoni M, Sebastiani P, Kohane IS: Cluster analysis of gene expression dynamics. Proc Natl Acad Sci USA 2002, 99(14):9121–9126. 10.1073/pnas.132656399PubMed CentralView ArticlePubMedGoogle Scholar
- Schliep A, Schonhuth A, Steinhoff C: Using hidden markov models to analyze gene expression time course data. Bioinformatics 2003, 19: i264-i272. 10.1093/bioinformatics/btg1036View ArticleGoogle Scholar
- Ma P, Castillo-Davis CI, Zhong W, Liu JS: A data-driven clusteirng method for time course gene expression data. Nuc Acids Res 2006, 34(4):1261–1269. 10.1093/nar/gkl013View 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/btg014View ArticlePubMedGoogle Scholar
- Heard NA, Holmes CC, Stephens DA: A quantitative study of gene regulation involved in the immune response of anopheline mosquitoes: An application of bayesian hierarchical clustering of curves. J Amer Statist Assoc 2006, 101(473):18–29. 10.1198/016214505000000187View ArticleGoogle Scholar
- Bar-Joseph Z, Gerber G, Jaakkola TS, Gifford DK, Simon I: Continuous representations of time series gene expression data. J Comput Biol 2003, 10(3–4):341–356. 10.1089/10665270360688057View ArticlePubMedGoogle Scholar
- Ernst J, Nau GJ, Bar-Joseph Z: Clustering short time series gene expression data. Bioinformatics 2005, 21(Suppl 1):i159-i168. 10.1093/bioinformatics/bti1022View ArticlePubMedGoogle Scholar
- Wang L, Ramoni MF, Sebastiani P: Clustering short gene expression profiles. Lecture Notes in Computer Science: Research in Computational Molecular Biology: 10th Annual International Conference, RECOMB 2006, Venice, Italy, April 2–5. Proceedings 2006, 60–68. [http://www.springerlink.com/content/t2447h8566270k01/]Google Scholar
- Kass RE, Raftery A: Bayes factors. J Amer Statist Assoc 1995, 90: 773–795. 10.2307/2291091View ArticleGoogle Scholar
- Diehn M, Alizadeh AA, Rando OJ, Liu CL, Stankunas K, Botstein D, Crabtree GR, Brown PO: Genomic expression programs and the integration of the CD28 costimulatory signal in t cell activation. Proc Natl Acad Sci USA 2002, 99(18):11796–11801. 10.1073/pnas.092284399PubMed CentralView ArticlePubMedGoogle Scholar
- Storey JD, Xiao W, Leek JT, Tompkins RG, Davis RW: Significance analysis of time course microarray experiments. Proc Natl Acad Sci 2005, 102(36):12837–12842. 10.1073/pnas.0504609102PubMed CentralView ArticlePubMedGoogle Scholar
- Heard NA, Holmes CC, Stephens DA, Hand DJ, Dimopoulos G: Bayesian coclustering of anopheles gene expression time series: Study of immune defense response to multiple experimental challenges. Proc Natl Acad Sci USA 2005, 102(47):16939–16944. 10.1073/pnas.0408393102PubMed CentralView ArticlePubMedGoogle Scholar
- Ng SK, McLachlan GJ, Wang K, Ben-Tovim Jones L, Ng S-W: A mixture model with random-effects components for clustering correlated gene-expression profiles. Bioinformatics 2006, 22(14):1745–1752. 10.1093/bioinformatics/btl165View ArticlePubMedGoogle Scholar
- Sebastiani P, Xie H, Ramoni MF: Bayesian analysis of comparative microarray experiments by model averaging. Bayes Anal 2006, 707–732. [http://ba.stat.cmu.edu/journal/2006/vol01/issue04/sebastiani.pdf]Google Scholar
- Ramoni M, Sebastiani P, Cohen PR: Bayesian clustering by dynamics. Mach Learn 2002, 47(1):91–121. 10.1023/A:1013635829250View ArticleGoogle Scholar
- Sebastiani P, Ramoni M, Kohane IS: Bayesian model-based clustering of gene expression dynamics. In The Analysis of Microarray Data: Methods and Software. Springer, New York, NY; 2003:409–427.Google Scholar
- Rand WM: Objective criteria for the evaluation of clustering methods. J Amer Statist Assoc 1971, 66(336):846–850. 10.2307/2284239View ArticleGoogle Scholar
- Hosack DA, Dennis G Jr, Sherman BT, Clifford Lane H, Lempicki RA: Identifying biological themes within lists of genes with EASE. Genome Biol 2003, 4(10):R70. 10.1186/gb-2003-4-6-p4PubMed CentralView ArticlePubMedGoogle Scholar
- Kluger Y, Basri R, Chang JT, Gerstein M: Spectral bioclustering of microarray data: Coclustering genes and conditions. Genome Res 2003, 13(4):703–716. 10.1101/gr.648603PubMed CentralView ArticlePubMedGoogle Scholar
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.