Logo of ijpgJournal's HomeManuscript SubmissionAims and ScopeAuthor GuidelinesEditorial BoardHome
Int J Plant Genomics. 2008; 2008: 231897.
Published online 2008 Apr 17. doi:  10.1155/2008/231897
PMCID: PMC2358942

Bayesian Functional Data Clustering for Temporal Microarray Data


We propose a Bayesian procedure to cluster temporal gene expression microarray profiles, based on a mixed-effect smoothing-spline model, and design a Gibbs sampler to sample from the desired posterior distribution. Our method can determine the cluster number automatically based on the Bayesian information criterion, and handle missing data easily. When applied to a microarray dataset on the budding yeast, our clustering algorithm provides biologically meaningful gene clusters according to a functional enrichment analysis.


Microarray technology enables the scientist to measure the mRNA expression levels of thousands of genes simultaneously. For a particular species of interest, one can make microarray measurements under many different conditions and for different types of cells (if it is a multicellular organism). Genes' expression profiles under these conditions often give the scientist some clues on biological roles of these genes. A group of genes with similar profiles are often “coregulated” or participants of the same biological functions.

When a series of microarray experiments are conducted sequentially during a biological process, we call the resulting dataset a “temporal” microarray dataset, which can provide insights on the underlying biology and help decipher the dynamic gene regulatory network. Clustering genes with similar temporal profiles is a crucial first step to reveal potential relationships among the genes.

Conventional clustering methods, such as the K-means and hierarchical clustering, do not take into consideration the correlation in the gene expression levels over time. Although it is possible to use a general multivariate Gaussian model to account for the correlation structure, such a model ignores the time order of the gene expressions. As evidenced in our example, the time factor is important in interpreting the results of gene expression clustering in temporal data. It is also possible to use an autoregression model to describe the gene expression time series, but such a model often requires stationarity, which is unlikely to hold in most temporal microarray data.

Recently, nonparametric analysis of data in the form of curves, that is, functional data, is subject to active research, see [1, 2] for a comprehensive treatment of functional data analysis; and curve-based functional clustering methods have emerged [37], but a rigorous assessment of the estimation precision is still lacking.

In this paper, we propose a Bayesian clustering method, which optimally combines the available information and provides a proper uncertainty measure for all estimated quantities. Our method is based on a mixture of mixed-effect smoothing splines models. For each cluster, we model its mean profile as a smoothing spline function and describe its individual gene's variation by a parametric random effect. Based on the theory of reproducing-kernel Hilbert spaces [8], we represent the mean expression curve as a linear combination of certain basis functions, which enables us to derive the full posterior distribution up to a normalizing constant. All the conditional distributions needed by a Gibbs sampler are also easy to compute and to sample from. Our method automatically takes care of the missing data and infers the number of clusters in the data. Using the method, we analyzed a microarray dataset of budding yeast, we found that the majority of the clusters we had obtained are enriched for known and expected biological functions.

Our method is not restricted to temporal microarray data, and can be applied to all curve clustering problems, especially for sparsely and irregularly sampled temporal data.


2.1. Mixed-effect representation of gene expression profile

Let the expression value of the ith gene at time t be yit. To accommodate missing data that occasionally occurs in microarray experiment, we denote t i = (t 1,…, t ni) and y i = (y i1,…, yini)T, where ni is the number of measurements of ith gene. Our mixed-effect smoothing spline model [9] for genes in one cluster is


where μ(t i) = (μ(t 1),…μ(tni))T is the cluster's mean profile, b i ~ N(0, B) is the random effect to capture the intragene correlation, Zi is the known design matrix for the random effect, and ϵi ~ N(0, σ 2 I) is the random error independent of b and of each other.

By taking different b vectors, we can accommodate different nonrandom effects. For example, when b i = bi and Zi = 1, the expression profile of the ith gene is parallel to the mean profile μ (Figure 1). If b i = (b i1, b i2)T and Zi = (1, t i), the difference between the ith gene profile and the mean profile is a linear function in time. More complicated structures such as periodicity can be modeled by letting the Zi be basis of a certain functional space.

Figure 1
A smoothing-spline mixed effect model for temporal gene expression.

By considering μ in a reproducing kernel Hilbert space H [subset, dbl equals] {μ : M(μ) < ∞} in which M(μ) is a square seminorm, we can represent μ as


where {sj} is a set consisting of all distinct {ti}, q is the number of {sj}, and RM is the kernel of H. The choice of M(μ) = ∫a 0 (d 2 μ/dt 2)2 dt yields the cubic smoothing spline with



where (·)+ = max(·, 0) [10].

Writing (2) in a vector-matrix form, we have


where Si is ni × m with the (i, ν)th entry ϕ ν(ti) and R is ni × q with the (i, j)th entry RM(ti, sj). Substituting (5) into (1), we have


Denoting y = (y T 1,…, y T n)T and S, R, Z, ϵ similarly, we have the matrix representation


where b = (b T 1,…, b T n)T ~ N(0, diag(B,…, B)).

The prior distributions are specified as follows:


where IG and IW are inverse-Gamma and inverse-Wishart distributions, respectively.

These priors lead to the following full conditional posteriors, which are used in our Gibbs sampler:


where Vd = (ST S/ σ 2 + diag(δ −1 1,…, δ −1 m))−1, Vb = (ZT Z/ σ 2 + diag(B −1,…, B −1))−1, Vc = (RT R/ σ 2 + 1/τ 2 I)−1, and SSR = (yS dR cZ b)T (yS dR cZ b).

2.2. The mixture model with unknown number of components

When more than one cluster is considered, we assume that the expression of the ith gene has a Gaussian mixture distribution:


where μk and Σk = ZBkZT + σ 2 I are the mean and covariance matrix for the kth component, as given by (7); p k is the fraction of kth component, and K is the number of Gaussian components.

2.3. Class labels and cluster numbers

To ease the computation, we introduce a “latent” membership labeling variable ji for the ith gene so that


When the number of Gaussian components K is known, we can get the joint posterior probability as


where J = (j 1,…, j n), μ = (μ 1,…, μ K), Σ = (Σ1,…, ΣK), and π(μ, Σ) is the joint prior distribution.

Since K is unknown, we used the following Bayesian information criterion (BIC):


where MK is the current model with parameters θK, θ^K is the estimate, and lK is the total number of parameters in our model. A small BIC score indicates the adequacy of the corresponding model. An alternative to our current approach (i.e., each clustering configuration is equally likely given the number of clusters K, and K is determined by BIC) is to use a Polya Urn prior (also called the “Chinese restaurant” process), which postulates that when a new member comes in, its a priori probability for joining an existing cluster of size mi is (mi + c)/(m + c), and for forming a new cluster of its own is c/(m + c), where m is the total number of existing members. This prior, however, favors unbalanced cluster configurations (e.g., very large and very small clusters) and may not be appropriate in our applications.

2.3.1. Gibbs Sampling from the Posterior

To complete our Bayesian analysis, we employ the Dirichlet prior Di (α 1,…, α K) for (p 1,…, p K), the cluster proportions. Thus, given the cluster indicator J, the posterior distribution of the p's is again a Dirichlet distribution.

Given μ 1,…, μ K, B 1,…, B K, σ 2, we have the conditional distribution of ji:


With an initial value of J, which gives rise to a partition of y : (y J 1,…, y J K), and the initial values of d k, b k, c k,Bk, where k = 1,…, K, as well as σ 2, we iterate the following iterative conditional sampling steps:

  1. for i = 1,…, n, draw a new ji from the conditional distribution from (14) to replace the old one;
  2. conditional on J, sequentially
    1. update d k by a draw from [d k | b k, c k, σ 2, δ, y J k], where k = 1,…, K,
    2. update b k from [b k | d k, c k, σ 2, Bk, y J k], where k = 1,…, K,
    3. update c k from [c k | d k, b k, σ 2, τ 2 k, y J k], where k = 1,…, K,
    4. update Bk ~ [Bk | b k], and τ 2 k ~ [τ 2 k | c k], where k = 1,…, K,
    5. update σ 2 ~ [σ 2 | d, b, c, y],
    6. update (p 1,…, p K) ~ Di(n 1 + α 1,…, nK + α K), where nj is the number of genes in the jth cluster.


To study oxygen-responsive gene network, Lai et al. [11] used cDNA microarray to monitor the gene expression changes of wild-type budding yeast (Saccharomyces cerevisiae) under aerobic condition in galactose medium. Under aerobic condition, the oxygen concentration was lowered gradually until oxygen was exhausted during a period of ten minutes. Microarray experiments were conducted at 14 time points under aerobic condition. A reference sample was obtained from a pooled RNA collected from all time points for hybridization.

For the analysis, Lai et al. [11] normalized gene expression after time 0 to gene expression of time 0 to set the common starting point. They identified 2388 genes whose expression is differentially expressed at one or more time points. Using our method, 2388 genes was clustered to 31 clusters, which yields the smallest BIC. FunSpec [12] was used for gene annotation and biological function enrichment analysis, where the Bonferroni-corrected functional enrichment P-values based on hypergeometric distributions are reported. We found 23 clusters out of 31 clusters discovered have biological functions over-represented. Among them, estimated mean gene expression profiles of three clusters are given in Figure 2.

Figure 2
Estimated mean expression curves for cluster A, B, and C (from top to bottom) discovered in the yeast aerobic expression data.

In cluster A, which consists of 40 genes, the estimated mean expression goes up progressively as oxygen level goes down, which suggests that the genes in this cluster were transiently upregulated in response to aerobisis. Accordingly, genes involved in stress response (function enrichment P-value = 10−4) as well as cell rescue and defense are over-represented in this cluster (function enrichment P-value = 10−4). Furthermore, genes involved in molecular functions of oxidoreductase and coproporphyrinogen oxidase are also presented, which explains the upregulation of the gene expression levels.

We have 92 genes in cluster B, where the estimated mean gene expression drops down at the beginning rapidly and then goes up gradually. In this cluster, 34 genes are involved in protein synthesis (function enrichment P-value ≤ 10−14). Moreover, ribosome biogenesis are also over-represented (function enrichment P-value ≤ 10−14). These processes were affected by oxygen level initially, but were quickly adjusted to high expression levels to maintain living of yeast.

Contrast to cluster B, cluster C (68 genes) consists of genes involved in galactose fermentation (function enrichment P-value = 10−4), carbon utilization (functional enrichment P-value = 10−2), and carbohydrate metabolism (function enrichment P-value ≤ 10−10). The initial upregulation of gene expression under aerobic condition can be partly explained by the fact that the cell increases the energy uptaking through the carbon utilization as oxygen level goes down; but as the oxygen level continues to drop down, these processes are replaced by the more energy-efficient processes, which drives the expression levels of genes to be downregulated.


Conventional clustering methods do not take into consideration the correlation in the gene expression levels over time. Multivariate Gaussian models and time series analysis cannot model the time factor and correlation properly. These limitations can be readily overcome by the full Bayesian approach developed here. Although certain prior distributions and the related hyperparameters need to be input by the user, we found the clustering results rather robust to variations in such inputs. Moreover, our Bayesian clustering algorithm serves as a platform to incorporate more biological knowledge. Open source R code is available at www.stat.uiuc.edu/~pingma/BayesianFDAClust.htm.


The authors thank C. I. Castillo-Davis for his help in designing Figure 1 and for the many constructive suggestions and discussion in the early stage of this research. The authors thank Ji Young Kim for designing the software website. The authors also thank Kurt Kwast for providing the yeast microarray data.


1. Ramsay JO, Silverman BW. New York, NY, USA: Springer; 2005. Functional Data Analysis.
2. Ramsay JO, Silverman BW. New York, NY, USA: Springer; 2002. Applied Functional Data Analysis: Methods and Case Studies.
3. James GM, Sugar CA. Clustering for sparsely sampled functional data. Journal of the American Statistical Association. 2003;98(462):397–408.
4. 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. [PubMed]
5. Luan Y, Li H. Model-based methods for identifying periodically expressed genes based on time course microarray gene expression data. Bioinformatics. 2004;20(3):332–339. [PubMed]
6. 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. Journal of the American Statistical Association. 2006;101(473):18–29.
7. Ma P, Castillo-Davis CI, Zhong W, Liu JS. A data-driven clustering method for time course gene expression data. Nucleic Acids Research. 2006;34(4):1261–1269. [PMC free article] [PubMed]
8. Wahba G. Vol. 59. Philadelphia, Pa, USA: SIAM; 1990. Spline Models for Observational Data. (CBMS-NSF Regional Conference Series in Applied Mathematics).
9. Gu C, Ma P. Optimal smoothing in nonparametric mixed-effect models. Annals of Statistics. 2005;33(3):1357–1379.
10. Gu C. New York, NY, USA: Springer; 2002. Smoothing Spline ANOVA Models.
11. Lai L-C, Kosorukoff AL, Burke PV, Kwast KE. Metabolic-state-dependent remodeling of the transcriptome in response to anoxia and subsequent reoxygenation in Saccharomyces cerevisiae. Eukaryotic Cell. 2006;5(9):1468–1489. [PMC free article] [PubMed]
12. Robinson MD, Grigull J, Mohammad N, Hughes TR. FunSpec: a web-based cluster interpreter for yeast. BMC Bioinformatics. 2002;3:3–35. [PMC free article] [PubMed]

Articles from International Journal of Plant Genomics are provided here courtesy of Hindawi Publishing Corporation
PubReader format: click here to try


Related citations in PubMed

See reviews...See all...

Cited by other articles in PMC

See all...


Recent Activity

Your browsing activity is empty.

Activity recording is turned off.

Turn recording back on

See more...