International Journal of Plant Genomics
Volume 2008 (2008), Article ID 231897, 4 pages
doi:10.1155/2008/231897
Research Article

Bayesian Functional Data Clustering for Temporal Microarray Data

Ping Ma,1 Wenxuan Zhong,1 Yang Feng,2 and Jun S. Liu2

1Department of Statistics, University of Illinois, Champaign, IL 61820, USA
2Department of Statistics, Harvard University, Cambridge, MA 02138, USA

Received 24 July 2007; Accepted 18 December 2007

Academic Editor: Nengjun Yi

Copyright © 2008 Ping Ma et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

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.

1. Introduction

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 [3–7], 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. Material and Methods

2.1. Mixed-Effect Representation of Gene Expression Profile

Let the expression value of the 𝑖 th gene at time 𝑡 be 𝑦 𝑖 𝑡 . To accommodate missing data that occasionally occurs in microarray experiment, we denote 𝐭 𝑖 = ( 𝑡 1 , … , 𝑡 𝑛 𝑖 ) and 𝐲 𝑖 = ( 𝑦 𝑖 1 , … , 𝑦 𝑖 𝑛 𝑖 ) 𝑇 , where 𝑛 𝑖 is the number of measurements of 𝑖 th gene. Our mixed-effect smoothing spline model [9] for genes in one cluster is 𝐲 𝑖 𝐭 = 𝝁 𝑖 + 𝑍 𝑖 𝐛 𝑖 + 𝝐 𝑖 , ( 1 ) where 𝝁 ( 𝐭 𝑖 ) = ( 𝜇 ( 𝑡 1 ) , … 𝜇 ( 𝑡 𝑛 𝑖 ) ) 𝑇 is the cluster's mean profile, 𝐛 𝑖 ∼ 𝑁 ( 0 , 𝐵 ) is the random effect to capture the intragene correlation, 𝑍 𝑖 is the known design matrix for the random effect, and 𝝐 𝑖 ∼ 𝑁 ( 0 , 𝜎 2 𝐼 ) is the random error independent of 𝐛 and of each other.

By taking different 𝐛 vectors, we can accommodate different nonrandom effects. For example, when 𝐛 𝑖 = 𝑏 𝑖 and 𝑍 𝑖 = 𝟏 , the expression profile of the 𝑖 th gene is parallel to the mean profile 𝝁 (Figure 1). If 𝐛 𝑖 = ( 𝑏 𝑖 1 , 𝑏 𝑖 2 ) 𝑇 and 𝑍 𝑖 = ( 𝟏 , 𝐭 𝑖 ) , the difference between the 𝑖 th gene profile and the mean profile is a linear function in time. More complicated structures such as periodicity can be modeled by letting the 𝑍 𝑖 be basis of a certain functional space.

231897.fig.001
Figure 1: A smoothing-spline mixed effect model for temporal gene expression.

By considering 𝜇 in a reproducing kernel Hilbert space ℋ ⊆ { 𝜇 ∶ 𝑀 ( 𝜇 ) < ∞ } in which 𝑀 ( 𝜇 ) is a square seminorm, we can represent 𝜇 as 𝜇 ( 𝑡 ) = 𝑚 𝜈 = 1 𝑑 𝜈 𝜙 𝜈 ( 𝑡 ) + 𝑞 𝑖 = 1 𝑐 𝑗 𝑅 𝑀 𝑠 𝑗 , 𝑡 , 𝑡 ∈ [ 0 , 𝑎 ] , ( 2 ) where { 𝑠 𝑗 } is a set consisting of all distinct { 𝑡 𝑖 } , 𝑞 is the number of { 𝑠 𝑗 } , and 𝑅 𝑀 is the kernel of ℋ . The choice of ∫ 𝑀 ( 𝜇 ) = 𝑎 0 ( 𝑑 2 𝜇 / 𝑑 𝑡 2 ) 2 𝑑 𝑡 yields the cubic smoothing spline with 𝜙 1 ( 𝑡 ) = 1 , 𝜙 2 ( 𝑡 ) = 𝑡 , ( 3 ) 𝑅 𝑀 𝑡 1 , 𝑡 2 = 𝑎 0 𝑡 1 − 𝑢 + 𝑡 2 − 𝑢 + 𝑑 𝑢 , ( 4 ) where ( ⋅ ) + = m a x ( ⋅ , 0 ) [10].

Writing (2) in a vector-matrix form, we have 𝝁 ( 𝐭 𝑖 ) = 𝑆 𝑖 𝐝 𝑖 + 𝑅 𝑖 𝐜 𝑖 , ( 5 ) where 𝑆 𝑖 is 𝑛 𝑖 × 𝑚 with the ( 𝑖 − , 𝜈 − ) th entry 𝜙 𝜈 ( 𝑡 𝑖 ) and 𝑅 is 𝑛 𝑖 × 𝑞 with the ( 𝑖 − , 𝑗 − ) th entry 𝑅 𝑀 ( 𝑡 𝑖 , 𝑠 𝑗 ) . Substituting (5) into (1), we have 𝐲 𝑖 = 𝑆 𝑖 𝐝 𝑖 + 𝑅 𝑖 𝐜 𝑖 + 𝑍 𝑖 𝐛 𝑖 + 𝝐 𝑖 . ( 6 ) Denoting 𝐲 = ( 𝐲 𝑇 1 , … , 𝐲 𝑇 𝑛 ) 𝑇 and 𝑆 , 𝑅 , 𝑍 , 𝝐 similarly, we have the matrix representation 𝐲 = 𝑆 𝐝 + 𝑅 𝐜 + 𝑍 𝐛 + 𝝐 , ( 7 ) where 𝐛 = ( 𝐛 𝑇 1 , … , 𝐛 𝑇 𝑛 ) 𝑇 ∼ 𝑁 ( 0 , d i a g ( 𝐵 , … , 𝐵 ) ) .

The prior distributions are specified as follows: 𝛿 𝐝 ∼ 𝑁 0 , d i a g 1 , … , 𝛿 𝑚 , 𝐜 ∼ 𝑁 0 , 𝜏 2 𝐼 , 𝜎 2 𝛼 ∼ I G 𝜎 2 , 𝛽 𝜎 2 , 𝜏 2 𝛼 ∼ I G 𝜏 2 , 𝛽 𝜏 2 , 𝜈 𝐵 ∼ I W 0 , 𝑆 0 − 1 , ( 8 ) 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: 𝐝 ∣ 𝐛 , 𝐜 , 𝜎 2 𝑉 , 𝜹 , 𝐲 ∼ 𝑁 𝑑 𝑆 𝑇 ( 𝐲 − 𝑅 𝐜 − 𝑍 𝐛 ) / 𝜎 2 , 𝑉 𝑑 , 𝐜 ∣ 𝐝 , 𝐛 , 𝜎 2 , 𝜏 2 𝑉 , 𝐲 ∼ 𝑁 𝑐 𝑅 𝑇 ( 𝐲 − 𝑆 𝐝 − 𝑍 𝐛 ) / 𝜎 2 , 𝑉 𝑐 , 𝐛 ∣ 𝐝 , 𝐜 , 𝜎 2 𝑉 , 𝐵 , 𝐲 ∼ 𝑁 𝑏 𝑍 𝑇 ( 𝐲 − 𝑆 𝐝 − 𝑅 𝐜 ) / 𝜎 2 , 𝑉 𝑏 , 𝜈 [ 𝐵 ∣ 𝐛 ] ∼ I W 0 𝑆 + 𝑛 , 0 + 𝑛 ∑ 𝑖 = 1 𝐛 𝑖 𝐛 𝑇 𝑖 − 1 , [ 𝜏 2 𝛼 ∣ 𝐜 ] ∼ I G 𝜏 2 + ( 𝑞 − 𝑚 ) / 2 , 𝛽 𝜏 2 + 𝐜 𝑇 , 𝐜 / 2 [ 𝜎 2 𝛼 ∣ 𝐝 , 𝐛 , 𝐜 , 𝐲 ] ∼ I G 𝜎 2 + 𝑛 / 2 , 𝛽 𝜎 2 , + S S R ( 9 ) where 𝑉 𝑑 = ( 𝑆 𝑇 𝑆 / 𝜎 2 + d i a g ( 𝛿 1 − 1 , … , 𝛿 𝑚 − 1 ) ) − 1 , 𝑉 𝑏 = ( 𝑍 𝑇 𝑍 / 𝜎 2 + d i a g ( 𝐵 − 1 , … , 𝐵 − 1 ) ) − 1 , 𝑉 𝑐 = ( 𝑅 𝑇 𝑅 / 𝜎 2 + 1 / 𝜏 2 𝐼 ) − 1 , and S S R = ( 𝐲 − 𝑆 𝐝 − 𝑅 𝐜 − 𝑍 𝐛 ) 𝑇 ( 𝐲 − 𝑆 𝐝 − 𝑅 𝐜 − 𝑍 𝐛 ) .

2.2. The Mixture Model with Unknown Number of Components

When more than one cluster is considered, we assume that the expression of the 𝑖 th gene has a Gaussian mixture distribution: 𝐲 𝑖 ∼ 𝑝 1 𝑁 ( 𝝁 1 , Σ 1 ) + ⋯ + 𝑝 𝐾 𝑁 ( 𝝁 𝐾 , Σ 𝐾 ) , ( 1 0 ) where 𝝁 𝑘 and Σ 𝑘 = 𝑍 𝐵 𝑘 𝑍 𝑇 + 𝜎 2 𝐼 are the mean and covariance matrix for the 𝑘 th component, as given by (7); 𝑝 𝑘 is the fraction of 𝑘 th component, and 𝐾 is the number of Gaussian components.

2.3. Class Labels and Cluster Numbers

To ease the computation, we introduce a “latent” membership labeling variable 𝐽 𝑖 for the 𝑖 th gene so that 𝐲 𝑖 ∣ 𝐽 𝑖 𝝁 = 𝑗 ∼ 𝑁 𝑗 , Σ 𝑗 . ( 1 1 ) When the number of Gaussian components 𝐾 is known, we can get the joint posterior probability as 𝑃 ( 𝐉 , 𝝁 , Σ ∣ 𝐲 ) = 𝜋 ( 𝝁 , Σ ) 𝑛 𝑖 = 1 𝑝 𝑗 𝑖 𝑁 ( 𝐲 𝑖 ∣ 𝝁 𝑗 𝑖 , Σ 𝑗 𝑖 ) , ( 1 2 ) where 𝐉 = ( 𝑗 1 , … , 𝑗 𝑛 ) , 𝝁 = ( 𝜇 1 , … , 𝜇 𝐾 ) , Σ = ( Σ 1 , … , Σ 𝐾 ) , and 𝜋 ( 𝝁 , Σ ) is the joint prior distribution.

Since 𝐾 is unknown, we used the following Bayesian information criterion (BIC): B I C = − 2 l o g 𝑝 𝐲 ∣ 𝑀 𝐾 , 𝜽 𝐾 + 𝑙 𝐾 l o g 𝑛 , ( 1 3 ) where 𝑀 𝐾 is the current model with parameters 𝜽 𝐾 , 𝜽 𝐾 is the estimate, and 𝑙 𝐾 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 𝐾 , and 𝐾 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 𝑚 𝑖 is ( 𝑚 𝑖 + 𝑐 ) / ( 𝑚 + 𝑐 ) , and for forming a new cluster of its own is 𝑐 / ( 𝑚 + 𝑐 ) , where 𝑚 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 , … , 𝛼 𝐾 ) for ( 𝑝 1 , … , 𝑝 𝐾 ) , the cluster proportions. Thus, given the cluster indicator 𝐉 , the posterior distribution of the 𝑝 's is again a Dirichlet distribution.

Given 𝝁 1 , … , 𝝁 𝐾 , 𝐵 1 , … , 𝐵 𝐾 , 𝜎 2 , we have the conditional distribution of 𝐽 𝑖 : 𝑝 𝐽 𝑖 = 𝑗 ∣ 𝝁 1 , … , 𝝁 𝐾 , 𝐵 1 , … , 𝐵 𝐾 , 𝜎 2 = 𝑝 , 𝐲 𝑗 𝐍 𝐲 𝑖 ∣ 𝝁 𝑗 , 𝑍 𝐵 𝑗 𝑍 𝑇 + 𝜎 2 𝐼 ∑ 𝐾 𝑘 = 1 𝑝 𝑘 𝐍 𝐲 𝑖 ∣ 𝝁 𝑘 , 𝑍 𝐵 𝑘 𝑍 𝑇 + 𝜎 2 𝐼 . ( 1 4 )

With an initial value of 𝐉 , which gives rise to a partition of 𝐲 ∶ ( 𝐲 𝐽 1 , … , 𝐲 𝐽 𝐾 ) , and the initial values of 𝐝 𝑘 , 𝐛 𝑘 , 𝐜 𝑘 , 𝐵 𝑘 , where 𝑘 = 1 , … , 𝐾 , as well as 𝜎 2 , we iterate the following iterative conditional sampling steps:(i) for 𝑖 = 1 , … , 𝑛 , draw a new 𝑗 𝑖 from the conditional distribution from (14) to replace the old one;(ii) conditional on 𝐉 , sequentially (a)update 𝐝 𝑘 by a draw from [ 𝐝 𝑘 ∣ 𝐛 𝑘 , 𝐜 𝑘 , 𝜎 2 , 𝜹 , 𝐲 𝐽 𝑘 ] , where 𝑘 = 1 , … , 𝐾 ,(b)update 𝐛 𝑘 from [ 𝐛 𝑘 ∣ 𝐝 𝑘 , 𝐜 𝑘 , 𝜎 2 , 𝐵 𝑘 , 𝐲 𝐽 𝑘 ] , where 𝑘 = 1 , … , 𝐾 ,(c) update 𝐜 𝑘 from [ 𝐜 𝑘 ∣ 𝐝 𝑘 , 𝐛 𝑘 , 𝜎 2 , 𝜏 2 𝑘 , 𝐲 𝐽 𝑘 ] , where 𝑘 = 1 , … , 𝐾 ,(d) update 𝐵 𝑘 ∼ [ 𝐵 𝑘 ∣ 𝐛 𝑘 ] , and 𝜏 2 𝑘 ∼ [ 𝜏 2 𝑘 ∣ 𝐜 𝑘 ] , where 𝑘 = 1 , … , 𝐾 ,(e) update 𝜎 2 ∼ [ 𝜎 2 ∣ 𝐝 , 𝐛 , 𝐜 , 𝐲 ] ,(f)update ( 𝑝 1 , … , 𝑝 𝐾 ) ∼ D i ( 𝑛 1 + 𝛼 1 , … , 𝑛 𝐾 + 𝛼 𝐾 ) , where 𝑛 𝑗 is the number of genes in the 𝑗 th cluster.

3. Results and Discussion

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 𝑃 -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.

fig2
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 𝑃 -value = 1 0 − 4 ) as well as cell rescue and defense are over-represented in this cluster (function enrichment 𝑃 -value = 1 0 − 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 𝑃 -value ≤ 1 0 − 1 4 ). Moreover, ribosome biogenesis are also over-represented (function enrichment 𝑃 -value ≤ 1 0 − 1 4 ). 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 𝑃 -value = 1 0 − 4 ), carbon utilization (functional enrichment 𝑃 -value = 1 0 − 2 ), and carbohydrate metabolism (function enrichment 𝑃 -value ≤ 1 0 − 1 0 ). 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.

4. Conclusions

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.

Acknowledgments

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

References

  1. J. O. Ramsay and B. W. Silverman, “Functional Data Analysis,” Springer, New York, NY, USA, 2005.
  2. J. O. Ramsay and B. W. Silverman, “Applied Functional Data Analysis: Methods and Case Studies,” Springer, New York, NY, USA, 2002.
  3. G. M. James and C. A. Sugar, “Clustering for sparsely sampled functional data,” Journal of the American Statistical Association, vol. 98, no. 462, 397 pages, 2003.
  4. Y. Luan and H. Li, “Clustering of time-course gene expression data using a mixed-effects model with B-splines,” Bioinformatics, vol. 19, no. 4, 474 pages, 2003.
  5. Y. Luan and H. Li, “Model-based methods for identifying periodically expressed genes based on time course microarray gene expression data,” Bioinformatics, vol. 20, no. 3, 332 pages, 2004.
  6. N. A. Heard, C. C. Holmes, and D. A. Stephens, “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, vol. 101, no. 473, 18 pages, 2006.
  7. P. Ma, C. I. Castillo-Davis, W. Zhong, and J. S. Liu, “A data-driven clustering method for time course gene expression data,” Nucleic Acids Research, vol. 34, no. 4, 1261 pages, 2006.
  8. G. Wahba, “Spline Models for Observational Data,” vol. 59 of CBMS-NSF Regional Conference Series in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 1990.
  9. C. Gu and P. Ma, “Optimal smoothing in nonparametric mixed-effect models,” Annals of Statistics, vol. 33, no. 3, 1357 pages, 2005.
  10. C. Gu, “Smoothing Spline ANOVA Models,” Springer, New York, NY, USA, 2002.
  11. L.-C. Lai, A. L. Kosorukoff, P. V. Burke, and K. E. Kwast, “Metabolic-state-dependent remodeling of the transcriptome in response to anoxia and subsequent reoxygenation in Saccharomyces cerevisiae,” Eukaryotic Cell, vol. 5, no. 9, 1468 pages, 2006.
  12. M. D. Robinson, J. Grigull, N. Mohammad, and T. R. Hughes, “FunSpec: a web-based cluster interpreter for yeast,” BMC Bioinformatics, vol. 3, 3 pages, 2002.