# Bayesian Functional Data Clustering for  Temporal  Microarray Data

Authors: Ping Ma; Wenxuan Zhong; Yang Feng; Jun S. Liu
Journal: International Journal of Plant Genomics
 (2008)

---

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


---

## Body

## 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 theith gene at time t be yit. To accommodate missing data that occasionally occurs
in microarray experiment, we denote ti=(t1,…,tni) and yi=(yi1,…,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(1)yi=μ(ti)+Zibi+ϵi,where μ(ti)=(μ(t1),…μ(tni))T is the
cluster's mean profile, bi∼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,σ2I) is the random
error independent of b and of each
other.By taking differentb vectors, we can
accommodate different nonrandom effects. For example, when bi=bi and Zi=1, the expression profile of the ith gene is
parallel to the mean profile μ (Figure 1). If bi=(bi1,bi2)T and Zi=(1,ti), 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 ℋ⊆{μ:M(μ)<∞} in which M(μ) is a square
seminorm, we can represent μ as(2)μ(t)=∑ν=1mdνϕν(t)+∑i=1qcjRM(sj,t),t∈[0,a],where {sj} is a set
consisting of all distinct {ti}, q is the number
of {sj}, and RM is the kernel
of ℋ. The choice of M(μ)=∫0a(d2μ/dt2)2dt yields the
cubic smoothing spline with(3)ϕ1(t)=1,ϕ2(t)=t,(4)RM(t1,t2)=∫0a(t1−u)+(t2−u)+du,where (⋅)+=max(⋅,0) [10].Writing (2) in a vector-matrix form, we
have(5)μ(ti)=Sidi+Rici,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(6)yi=Sidi+Rici+Zibi+ϵi.Denoting y=(y1T,…,ynT)T and S, R, Z, ϵ similarly, we
have the matrix representation(7)y=Sd+Rc+Zb+ϵ,where b=(b1T,…,bnT)T∼N(0,diag(B,…,B)).The prior distributions are specified as
follows:(8)d∼N(0,diag(δ1,…,δm)),c∼N(0,τ2I),σ2∼IG(ασ2,βσ2),τ2∼IG(ατ2,βτ2),B∼IW(ν0,S0−1),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:(9)[d∣b,c,σ2,δ,y]∼N(VdST(y−Rc−Zb)/σ2,Vd),[c∣d,b,σ2,τ2,y]∼N(VcRT(y−Sd−Zb)/σ2,Vc),[b∣d,c,σ2,B,y]∼N(VbZT(y−Sd−Rc)/σ2,Vb),[B∣b]∼IW(ν0+n,(S0+∑i=1nbibiT)−1),[τ2∣c]∼IG(ατ2+(q−m)/2,βτ2+cTc/2),[σ2∣d,b,c,y]∼IG(ασ2+n/2,βσ2+SSR),where Vd=(STS/σ2+diag(δ1−1,…,δm−1))−1, Vb=(ZTZ/σ2+diag(B−1,…,B−1))−1, Vc=(RTR/σ2+1/τ2I)−1, and SSR=(y−Sd−Rc−Zb)T(y−Sd−Rc−Zb).


## 2.2. The Mixture Model with Unknown  Number of Components

When more than one cluster is considered, we assume
that the expression of theith gene has a
Gaussian mixture distribution:(10)yi∼p1N(μ1,Σ1)+⋯+pKN(μK,ΣK),where μk and Σk=ZBkZT+σ2I are the mean
and covariance matrix for the kth component,
as given by (7); pk 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 variableJi for the ith gene so
that(11)yi∣Ji=j∼N(μj,Σj).When the number of Gaussian
components K is known, we
can get the joint posterior probability as(12)P(J,μ,Σ∣y)=π(μ,Σ)∏i=1npjiN(yi∣μji,Σji),where J=(j1,…,jn), μ=(μ1,…,μK), Σ=(Σ1,…,ΣK), and π(μ,Σ) is the joint
prior distribution.SinceK is unknown, we
used the following Bayesian information criterion (BIC):(13)BIC=−2logp(y∣MK,θ^K)+lKlogn,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 (p1,…,pK), the cluster proportions. Thus, given the cluster
indicator J, the posterior distribution of the p's is again a
Dirichlet distribution.Givenμ1,…,μK,B1,…,BK,σ2, we have the
conditional distribution of Ji:(14)p(Ji=j∣μ1,…,μK,B1,…,BK,σ2,y)=pjN(yi∣μj,ZBjZT+σ2I)∑k=1KpkN(yi∣μk,ZBkZT+σ2I).With an initial value ofJ, which gives rise to a partition of y:(y1J,…,yKJ), and the initial values of dk, bk, ck, Bk, where k=1,…,K, as well as σ2, we iterate the following iterative conditional
sampling steps:(i)
fori=1,…,n, draw a new ji from the
conditional distribution from (14) to replace the old one;(ii)
conditional onJ, sequentially
(a)
updatedk by a draw from [dk∣bk,ck,σ2,δ,ykJ], where k=1,…,K,(b)
updatebk from [bk∣dk,ck,σ2,Bk,ykJ], where k=1,…,K,(c)
updateck from [ck∣dk,bk,σ2,τk2,ykJ], where k=1,…,K,(d)
updateBk∼[Bk∣bk], and τk2∼[τk2∣ck], where k=1,…,K,(e)
updateσ2∼[σ2∣d,b,c,y],(f)
update(p1,…,pK)∼Di(n1+α1,…,nK+αK), where nj is the number
of genes in the jth 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 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.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 enrichmentP-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 enrichmentP-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 enrichmentP-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.


## 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 atwww.stat.uiuc.edu/~pingma/BayesianFDAClust.htm.



---

*Source: 231897-2008-04-17.xml*