# Bayesian Mixture Model Analysis for Detecting Differentially Expressed Genes

Authors: Zhenyu Jia; Shizhong Xu
Journal: International Journal of Plant Genomics
 (2008)

---

## Abstract

Control-treatment design is widely used in microarray gene
expression experiments.
The purpose of such a design is to detect genes that express
differentially between the control and the treatment. Many
statistical procedures have been developed to detect
differentially expressed genes, but all have pros and cons and
room is still open for improvement. In this study, we propose a
Bayesian mixture model approach to classifying genes into one of
three clusters, corresponding to clusters of downregulated,
neutral, and upregulated genes, respectively. The Bayesian method
is implemented via the Markov chain Monte Carlo (MCMC) algorithm.
The cluster means of down- and upregulated genes are sampled from
truncated normal distributions whereas the cluster mean of the
neutral genes is set to zero. Using simulated data as well as data
from a real microarray experiment, we demonstrate that the new
method outperforms all methods commonly used in differential
expression analysis.


---

## Body

## 1. Introduction

Current
microarray technology allows us to measure the expression of thousands of genes
simultaneously in a small chip. Many advanced statistical methods have been
developed to analyze data generated from this technology. For example, a
Bayesian model has been proposed to normalize gene expression data and select
genes that are differentially expressed [1]. Jörnsten et al. [2] developed a
Bayesian method to impute missing values prior to any statistical analysis.
Efron et al. [3], Broët et al.
[4], Edwards et al. [5],  Do et al. [6]
used a Bayesian approach to identify differentially  expressed  genes. Methods of
[7–9] are examples for Bayesian clustering
analysis.In this study, we focus on detecting differentially
expressed genes in two-sampled designs, in which only two conditions (control
and treatment) are examined and each condition is replicated several times. The
data are assumed to have been properly imputed and normalized as needed. Much
effort has been made on this study to detect genes whose expression levels
respond to the treatment.Numerous methods have been suggested based on theP-values
of a test statistic. The P-values are obtained from separate tests and
reported for all genes [10, 11]. Cui and Churchill [12] reviewed the test statistics
for differential expression for microarray experiments. The gene-specific
summary statistic provides a basis for the rank ordering of genes and the
creation of a short list of genes declared as differentially expressed.
However, genes are studied separately in the sense that calculating the test statistic
for one gene does not use information from other genes. In limma (linear model
for microarray data), Smyth [13] cleverly borrowed information from the
ensemble of genes to make inference for individual gene based on the moderate
t-statistic. Some other researchers also took advantages of shared information
by examining data jointly. Efron et al. [3]  proposed a mixture model
methodology implemented via an empirical Bayes approach. For each gene, the
expression levels were mapped to a single z-score, which is assumed to arise
from two distributions (affected and unaffected). Pan [14] used a mixture model
to cluster genes based on gene specific t-statistic. They are examples where
the idea of clustering was first used in differential expression analysis
(clustering methods are usually applied to microarray data generated from
multiple conditions). Similarly, the methods of [4–6, 15] all used a mapped
quantity, like the z-score or the t-statistic, to classify genes into different
groups. These methods are different from the proposed Bayesian mixture model in
that genes were clustered based on a summary statistic for these methods while
our Bayesian method clusters genes based on the pattern of expression. It is
more desirable to use the expression profile than to use a summary score for
cluster analysis because information is bound to be lost when mapping the
expression profile into a single score. Recently, Newton et al. [16]
developed a
new method in which gene expression are directly used for classification.
Newton et al. [16] used a hierarchical model that consists of two levels. The
first level of the model describes the conditional distributions of the
observed measurement of gene expression given the means of the control (μ1i) and the
treatment (μ2i), where i indexes the
gene. The second level describes the distribution of μ1i and μ2i jointly by a
mixture of three multivariate distributions with constrained parameters of μ1>μ2, μ1=μ2, and μ1<μ2. Genes become linked by virtue of having μ1i and μ2i drawn from a
common distribution. Parameter estimation was accomplished via the
Expectation-Maximization (EM) algorithm [17].For the Bayesian clustering method presented in this
study, the observed gene expression levels are described by a regression model
as done by Gusnanto et al. [15]. For each gene, the irrelevant intercept is removed
by a special normalization scheme. The slope of the regression represents the
difference of expression under the two conditions. The Bayesian method is
implemented via the Markov chain Monte Carlo (MCMC) algorithm. The regression
coefficient of each gene is assumed to be sampled from a mixture of three
normal distributions with constrained parameters. The three distributions are N(βk,νk) for k=1,2,3, where β1<0, β2=0, and β3>0 are the
constrained parameters. The proposed new method actually turns the problem of a
complicated multivariate mixture distribution of Newton et al. [16] into that of
a univariate mixture distribution.The new Bayesian method proposed in this study (Method
I) is compared to five different methods that are commonly used in differential
expression analysis. These four methods, called Methods II, III, IV, V, and VI
respectively, are described below. Method II is the regularizedt-test
(SAM, Tusher et al.
[10]), in which a t-like score is calculated for each gene.
Genes with scores greater than an adjusted threshold are said to have altered
expressions under the treatment. Permutations are used to calculate the false
discovery rate (FDR, Benjamini and Liu [18]). Method III is the model-based cluster
analysis of Pan [14] where the variable used for clustering is the t-test
statistic. Genes that are assigned into a cluster with a mean t value
significantly different from zero are declared to be differentially expressed.
Method IV is the hierarchical mixture model method of Newton et al. [16], where
the expected expression values of the ith gene under
the control (μ1i) and the
treatment (μ2i) are assumed
to arise from a mixture of three distributions. Clustering is actually
conducted based on the latent parameters (μ1i and μ2i). The cluster
of genes with μ1i≠μ2i are declared as
differentially expression genes. A short list of significant genes is created
based on FDR at 0.05. Method V is the linear model and empirical Bayes method
(limma) of Smyth [13]. In limma, a hierarchical linear model is used to
describe the gene expression levels. A moderated t-statistic, with extra
information borrowed from other genes, is calculated for each gene. Adjusted P-values
are calculated for genes based on the moderated t-statistics to achieve FDR
control. Method VI is the mixture mixed model analysis of Gusnanto et al. [15], where
a similar linear model is used to describe gene expressions. However, the
variable used for clustering is the mean difference between the treatment and the
control. The mathematical differences between the proposed method and the four
methods compared will be evident after the new method is introduced. Advantages
of the new method over the four methods will be demonstrated in Section 3 and
further addressed in the Section 4.


## 2. Methods

## 2.1. Linear Model

LetYij(i=1,…,mandj=1,…,N) be the log transformed expression level for gene i from chip j, where m is the total
number of genes and N is the total
number of chips. Let Xj be a
dichotomous variable indicating the status of the chip with Xj=0 for the control
and Xj=1 for the
treatment. We use the following regression model to describe Yij: (1)Yij=αi+Xjγi+εij,  where αi is the
intercept of the regression model, γi is the
regression coefficient and εij is the residual
error. Note that the regression coefficient represents the difference in gene
expression between the control and the treatment (μ2−μ1). A large γi implies that
gene i is
differentially expressed. The intercept αi represents the
mean expression level of gene i in the control (μ1). The intercept αi is irrelevant
to the differential expression analysis and thus should be removed from the
model (see next section for detail). The regression coefficient γi is assumed to
be sampled from one of three normal distributions, corresponding to the
downregulated genes, the neutral genes, and the upregulated genes,
respectively. The residual error is assumed to be N(0,σ2) distributed. In
matrix notation, model (1) can be expressed as (2)Yi=Wαi+Xγi+εi,  where Yi={Yij}j=1N is a column
vector for the expressions of gene i across all
samples, X={Xj}j=1N is the
incidence matrix (a column vector) for the linear model, εi={εij}j=1N is a vector for
the residual errors, and W is an N×1 unity vector,
that is, a vector of 1 for all elements.


## 2.2. Normalization

As stated earlier,αi is irrelevant
to the differential expression analysis. We propose to remove the effect of α by using a
linear contrast. This process, which is called normalization, is similar to the
way of removing the fixed effects in the restricted maximum likelihood (REML)
method for variance component analysis [19]. For each gene, we
subtract the observed expression level by the average level across all chips,
that is, (3)Yi*=Yi−WY¯i, where Y¯i.=N−1∑j=1NYij is the mean
expression for gene i across all the
chips. The mean expression can also be expressed in matrix notation as (4)Y¯i.=(WTW)−1WTYi.. Therefore, the normalized expression is described as a
linear combination of the original expressions (5)Yi*=[I−W(WTW)−1WT]Yi=LYi,  where (6)L=I−W(WTW)−1WT.  We can prove
that LW=0. Therefore, (7)Yi*=LYi=LXγi+Lεi. Let X*=LX and εi*=Lεi, we have (8)Yi*=X*γi+εi*.  Now the
intercept has been removed by using linear contrasts LYi in place of the
original Yi. The residual variance-covariance matrix of the
linear contrasts is Var(εi*)=LLTσ2. Note that the rank of matrix LLT is N−r(W)=N−1, indicating that the inverse of matrix LLT does not exist.
The model-based cluster analysis, however, requires the inverse. Therefore, we
delete the last row of matrix L to form a new (N−1)×N matrix, called L*, to build the linear contrasts. This treatment is
equivalent to deleting the last element of Yi*, that is, the dimension of Yi* is (N−1)×1. The dimensions of matrix X* and εi* also change
into (N−1)×1 accordingly.
Let R=L*(L*)T be an (N−1)×(N−1) matrix of full
rank. We now have Var(εi*)=Rσ2, which will be used in the model-based cluster
analysis.To simplify the notation, let us defineyi=Yi*, x=X* and ϵi=εi*, all with a dimension of (N−1)×1. The new model with the intercept removed becomes (9)yi=xγi+ϵi.  Note that the
special way of normalization described above only applies to linear effects.
Different methods should be used for adjusting nonlinear effects. In subsequent
analysis, we assume that normalization via linear contrasts has been conducted
and thus model (9) will be the basis for parameter estimation and statistical
inference.


## 2.3. Mixture Model and the Bayesian Setup

Conditional onγi, the probability density of yi is (10)p(yi∣γi,σ2)=ϕn(yi;xγi,Rσ2),  where n=N−1, ϕn(yi;xγi,Rσ2) is the
n-dimensional normal density of variable yi with mean xγi and
variance-covariance matrix Rσ2. The gene-specific effect γi is assumed to
be sampled from one of three normal distributions (clusters), N(βk,νk) for k=1,2,3. We define variable zi as the cluster
assignment of gene i, where zi=k if gene i comes from
cluster k. For notational simplicity, we define η(zi,k) for k=1,2,3 as a redundant
variable to zi, where η(zi,k)=1 if zi=k and η(zi,k)=0 otherwise. Let η(zi)={η(zi,k)}k=13 for η(zi,k)=0,1 and ∑k=13η(zi,k)=1. The new variable η(zi) will be used in
place of zi in subsequent
derivation. Let β={βk}k=13 and ν={νk}k=13. The density of the mixture distribution of γi is (11)p(γi∣η(zi),β,ν)=∑k=13η(zi,k)ϕ1(γi;βk,νk).  Note that
variable η(zi,k) is unknown and
it is treated as missing value. In fact, inferring the probability distribution
of η(zi,k) is the main
purpose of the proposed mixture model analysis. Define π={πk}k=13 as the mixing
proportions of the three components for πk≥0 and ∑k=13πk=1. The distribution of η(zi) is multinomial
with one observation, (12)p(η(zi)∣π)=∏k=13πkη(zi,k).We have introduced the probability densities of the
data, the regression coefficients, and the missing values for a single gene. We
now combine the densities of individual genes to form the joint density of all
genes. The probability density of the datay={yi}i=1m is (13)p(y∣γ,σ2)=∏i=1mp(yi∣γi,σ2).  The density of γ={γi}i=1m is (14)p(γ∣η,β,ν)=∏i=1mp(γi∣η(zi),β,ν).  The density of η={η(zi)}i=1m is (15)p(η∣π)=∏k=13πkmk,  where mk=∑i=1mη(zi,k) for ∑k=13mk=m. The joint density of (y,γ,η) is (16)p(y,γ,η∣θ)=p(y∣γ,σ2)p(γ∣η,β,ν)p(η∣π),  where θ=(σ2,β,π,ν) is the
parameter vector.The next step of the analysis is to find a suitable
prior forθ. The following vague prior is chosen for each of the
variance components, σ2∼Inv−χ2(0,0) and νk∼Inv−χ2(0,0) for k=1,2,3. They are called the scaled inverse chi-square
distributions with zero degrees of freedom and zero-scale parameter. The actual
forms of the priors are p(σ2)=1/σ2 and p(νk)=1/νk for k=1,2,3. The scaled inverse chi-square distribution is
conjugate, and thus the posterior distribution of the variance is also scaled
inverse chi-square. The three clusters of genes are distinguished by the
overall patterns of responses to the treatment. The first cluster consists of
all the downregulated genes, the second cluster represents all the neutral genes,
and the third cluster contains all the upregulated genes. The characteristics
of the three clusters can be represented by enforcing the following constraints
on the means of the three clusters: β1<0, β2=0, and β3>0. Therefore, the prior distributions of the means of
the three clusters are p(βk)=1\  (βk∈Ωk). Here, we adopt a special notation to indicate that p(βk)=1 if (βk∈Ωk) and p(βk)=0 if (βk∉Ωk), where Ω1≡(β1<0), Ω2≡(β2=0), and Ω3≡(β3>0). The prior for π is the
multivariate generalization of the Beta distribution p(π∣δ)=D(π∣δ,δ,δ), called the Dirichlet distribution, where δ=1 is used in this
study. The joint prior for all the parameters is p(θ)=(1/σ2)∏k=13(1/νk). The posterior distribution of all the unknowns p(γ,η,θ∣y) is proportional
to the joint distribution of all variables, (17)p(γ,η,θ∣y)∝p(y,γ,η,θ)=p(y∣γ,σ2)p(γ∣η)p(η∣π)p(θ).  Statistical
inference of this distribution is the theme of the Bayesian analysis. There is
no explicit form for the joint distribution (17). Therefore, we draw observations of the unknowns
from the conditional distributions. Fortunately, with the current Bayesian
setup, the conditional distribution of any single variable, given all other
variables, has a known distribution. Therefore, the MCMC process can be
proceeded exclusively using the Gibbs sampler.


## 2.4. Markov Chain Monte Carlo

The detailed
steps of the Markov chain Monte Carlo (MCMC) process are described as follows.Step 1.
Sett=0 and initialize
all variables (γ(t),η(t),β(t), σ2(t),ν(t),π(t)).Step 2.
Simulateγi(t+1) from γi∼N(γ¯i,sγi2) conditional on zi=k, where (18)sγi2=[xTR−1x+σ2/νk]−1σ2γ¯i=[xTR−1x+σ2/νk]−1[xTR−1y+(σ2/νk)βk],  where all variables in the right-hand side take the
most current values. In this step, the most current values of the variables in
the right-hand side are the values at iteration t. This statement also holds for subsequent steps,
except that the most current values of the variables in the right-hand side are
values at iteration t or t+1, depending on whether the variable has been updated
or not in the current sweep.Step 3.
Simulateη(zi)(t+1) from (19)η(zi)∼Multinomial(1,πi1,πi2,πi3), a multinomial distribution with one observation, where πik is
(20)πik=πkp(γi∣η(zi,k),βk,νk)∑k′=13πk′p(γi∣η(zi,k′),βk′,νk′).Step 4.
Simulateβk(t+1) from a
truncated normal distribution (21)βk∼N(β¯k,sβk2)1(βk∈Ωk). The mean and the variance of the normal distribution
are (22)β¯k=1∑i=1mη(zi,k)∑i=1mη(zi,k)γi,sβk2=νk∑i=1mη(zi,k),  respectively. Note that β2=0 is enforced and
no sampling for β2 is required.
The general inverse transformation method [20] is used to sample β1 and β3 from the
corresponding truncated normal distributions.Step 5.
Simulateσ2(t+1) from (23)σ2∼Inv−χ2(mn,∑(yi−xγi)TR−1(yi−xγi)), where m×n is the degree of freedom, and the term with  summation is the scale parameter of the
inverse chi-square distribution.Step 6.
Simulateνk(t+1) from (24)νk∼Inv−χ2(mk,∑η(zi,k)(γi−βk)2), where mk=∑i=1mη(zi,k).Step 7.
Simulateπ(t+1) from (25)π∼Dirichlet(m1+1,m2+1,m3+1).So far, every variable has been updated. We then
incrementt by 1(t=t+1) and repeat
Steps 3–7 until the chain gets sufficiently long to allow reliable post-MCMC
inference about the parameters of interest. Genei is assigned to group k if η¯(zi,k)=max⁡{η¯(zi,k′)}k′=13, where η¯(zi,k) is the
posterior mean of η(zi,k) calculated from
the posterior sample.Schemes for sampling variables from the aforementioned
distributions are discussed by Gelman et al. [21].


## 3. Applications

## 3.1. Analysis of Simulated Data

We simulated
the expression ofm=1000 genes from six
different groups on N=16 microarray
chips. Chips 1–8 represent the control, and chips 9–16 represent the treatment.
The true parameter values of the six groups of genes are shown in Table 1. If we ignore the intercepts, the six groups of genes
really come from three functional clusters: Cluster 1 represents the (m1=20) downregulated
genes (Group 3); Cluster 2 contains the (m2=960) neutral genes
(Groups 2 and 4); and Cluster 3 consists of the (m3=20) upregulated
genes (Groups 1, 5, and 6). The true parameters of the original six groups and
the parameters for the combined three clusters are given in Table 1. Note that there are actually three components (Group
1, 5, and 6) for the combined Cluster 3. The weighted β of three
simulated components is treated as the true β for this
cluster, that is, 1.2×0.25+0.8×0.25+0.9×0.5=0.95. The expression patterns are shown in Figure 1 for
the original six groups and Figure 2 for the combined three clusters.Table 1
Parameters used in the simulation experiment and their
estimates from the Bayesian mixture model analysis.
Parameterαkβkπkνkσ2TrueGroup 13.51.20.0050.010.03Group 2400.3450.04Group 34.5−10.020.02Group 4000.6150.04Group 500.80.0050.01Group 600.90.010.01True
(combined)Cluster 1—−10.0200.020.03Cluster 2—00.9600.04Cluster 3—0.950.0200.03EstimatedCluster 1—−0.9840.0230.030.03Cluster 2—00.9370.04Cluster 3—0.7020.0400.11Original expression patterns for the six
simulated groups of genes. In each plot, the first half represents the observed
data from chips 1–8 and the second half represents the observed data from
chips 9–16.Expression patterns for the three combined
clusters after normalization. See Figure1 for the legends.The data were analyzed using the Bayesian mixture model
analysis reported in this study (Method I). We set the number of iterations
equal to 12000 in the MCMC process. The results generated from the first 2000
iterations were discarded for the burn-in period. For the remaining 10000
iterations, we saved one observation in every 20 iterations to reduce series
correlation between consecutive observations. Therefore, the posterior sample
contained 500 observations. The posterior mean of each variable is considered
as the Bayesian estimate. Hypothesis tests were only performed on the cluster
meansβ rather than on
individual genes. Genes assigned into a significant cluster were simultaneously
declared to be differentially expressed. From Table 1, we can see that the
estimated parameters agree well with the true parameters. Genes are correctly
assigned into three clusters except that two genes in Cluster 2 (the neutral
cluster) are incorrectly assigned into the other two clusters (down- and
upregulated clusters) (see Table 2). Upon enforcing constraints on the cluster
means, we are able to fix the number of clusters to three, rather than finding
the optimal number of clusters using the Bayes factor or its BIC
[22] approximation.Table 2
Numbers of genes assigned into each of the three
clusters for six different methods of differential expression analysis. (The sum of each
column within a method represents the true number of genes simulated from that
cluster and the sum of each row represents the number of genes assigned into
that cluster.)
MethodEstimateTrueSumType I errorCluster 1Cluster 2Cluster 3ICluster 12010210.002Cluster 209580958Cluster 3012021IICluster 12020502250.459Cluster 205190519Cluster 3023620256IIICluster 100000.033Cluster 209280928Cluster 320322072IVCluster 120370570.073Cluster 208900890Cluster 30332053VCluster 12015301730.325Cluster 206480648Cluster 3015920179VICluster 12050250.011Cluster 209490949Cluster 3062026Sum (method)20960201000An interesting question is what would happen if the
cluster means are not constrained? To answer this question, we also analyzed
the simulated data with unconstrainedβ. We varied the number of clusters from 3 to 6 and
found out that 4 was the optimal number of clusters (maximum value of Bayes
factor). One cluster mean was close to zero. Another cluster had a negative
mean. The third and fourth clusters all had positive means but the two means
were almost identical. Therefore, Clusters 3 and 4 may be well considered as a
single cluster. When treated as three clusters, the results were similar to
those reported earlier when the cluster means were restricted (data now shown).For comparison, the data were also analyzed using the
five methods (Methods II–VI) described in Section1. The numbers of genes
classified into each of the three clusters are given in Table 2 for all the six
methods, including Method I (the new method). Cutoff value 0.05 was used for
gene detection in Method II and Method V to achieve FDR control. All methods
have successfully found the 40 truly differentially expressed genes. However,
Methods II, III, IV, V, and VI detected more false-positive genes than Method
I. The model-based cluster analysis of Pan (Method III) required BIC to decide
the optimal number of clusters. For simplicity of comparison, we only used the
result of Method III under three clusters. We noticed that Method III always
put the down- and upregulated genes into a cluster with a large variance, but
placed the neutral genes into the other two cluster with small variances. The
two clusters with small variances can be combined as a single cluster. Because
all the methods have successfully detected the 40 nonneutral genes, none of
them have suffered the Type II error. However, different methods tended to have
different numbers of genes from the neutral cluster misplaced into the
functional clusters. This means that different methods have different Type I
errors.LetN(k,k′) be the number
of genes from cluster k assigned into
cluster k′ for k,k′=1,2,3. Then N(k)=∑k′=13N(k,k′) is the number
of genes in cluster k. Recall that Clusters 1 and 3 contain differentially
expressed genes and Cluster 2 is reserved for the neutral genes. The empirical
Type I error is then defined as (26)Type I error=N(2,1)+N(2,3)N(2). The empirical Type II error is defined as (27)Type II error=N(1,2)+N(3,2)N(1)+N(3). The empirical power is defined as (28)Power=1−Type II error. The empirical Type I errors for all the five methods
are listed in Table 2. Method I (the new method) have the smallest Type I
error. Note that this method of validation is valid only for simulation study
where the true gene assignment is known.


## 3.2. Analysis of Mice Data

We analyzed a
cDNA microarray dataset published by Dudoit et al. [11]. The data are publicly
available on the website (http://www.stat.berkeley.edu/users/terry). The data
consist of expression measurements of 6342 genes from a study of lipid
metabolism in mice [23]. This experiment was set up to find genes
that are differentially expressed in the livers of mice with very low HDL
cholesterol levels (treatment) in contrast to a group of normal inbred mice
(control). The treatment group consists of eight scavenger receptor BI (SR-BI)
transgenic mice and the control group consists of eight normal inbred mice.We first analyzed the data with the Bayesian mixture
model method proposed in this study (Method I). We used exactly the same setup
of the MCMC as that used in the simulation experiment. The estimated parameters
are given in Table3. Only β3 (cluster mean
of up-regulated genes) is significantly different from zero. The mixing
proportions of the three clusters show that majority of the genes (82.46%) are
neutral. Note that the estimated mixing proportions (π^1=0.1724, π^2=0.8246, π^3=0.0030) are some
nuisance parameters, which do not necessarily match the actual numbers of genes
assigned to the three clusters (m1=0, m2=6332, m3=10). Report
should be made based on the actual numbers of genes assigned to the clusters.Table 3
Parameters estimated for the mice data using the
Bayesian mixture model analysis.
βkπkνkσ2Cluster 1−0.00650.17240.06160.0904Cluster
200.82460.0006Cluster
30.64180.00302.0266All six methods that were used to analyze the
simulated data were applied here to analyze the mice data. The numbers of genes
assigned to the three clusters are given in Table4 for all the five methods.
Clearly, Methods I and VI were similar to each other and Methods II, III, IV,
and V placed much more genes into the nonneutral clusters, a phenomenon which
has been observed early in the simulation experiment.Table 4
Numbers of genes assigned to the three clusters for
the mice data for six different methods.
MethodCluster 1Cluster 2Cluster 3(down)(neutral)(up)I0633210II406182120III6662760IV13630029V12629337VI0632913In the original study [23], individual P-values
were calculated based on the t-test statistics for all the genes and then the
adjusted P-values were calculated using the method described in
Westfall and Young [24]. Five genes (represented by eight array elements) were
detected by Callow et al. [23] as differentially expressed between the
treatment and the control (see Table 5). Note that several array elements may
represent the same gene. We found that gene EST AI746730 is not included in the
dataset and thus was not analyzed with any of the methods (marked as NA in
Table 5).Table 5
Some differentially expressed genes for the mice data
detected by six different methods.
Gene or elementMethodIIIIIIIVVaSR-BI (ID: 3)√√√√SR-BI (ID: 783)√√√√√√SR-BI (ID: 1581)√√√√√√Glutathione-S-transferase√√√√√√β-globin√√√√Cytochrome P450 2B10EST AI746730NANANANANANAbGrowth factor receptor bound pr2√√√√√√Creatine Kinase, muscle√√√√√Ubiquitin-conjugating enzyme E2H√√√√√√Capping protein alpha 2√√√√√√Diff. Ass. Prot. 13 kD√√√√√√Blank (ID: 2374)√√√√√√Blank (ID: 2375)√Blank (ID: 2377)√√Blank (ID: 2379)√√√Blank (ID: 2384)√a: genes
reported in Callow et al. [23]; b: a subset of genes not reported in Callow et al.  [23].The Bayesian mixture model analysis detected 10 genes
(or spots), four of which were reported in Callow et al. [23]. The expressions
of the 10 genes are plotted in Figure 3 (the first two rows). Significant
differences between the control and treatment groups can be easily seen. Note
that element 2374 (the first plot of the second row of Figure 3) was a “blank”
spot on the chip. However, it was detected by all methods due to the  remarkable
difference between the control and the treatment. The third row of the plots
represents five randomly selected genes which were not detected by any of the
methods. Genes of this kind have the same expression levels under the control
and treatment. The last row of the plots (in Figure 3) shows the expression
patterns of five genes that were not detected by the Bayesian mixture model
analysis (Method I), but detected by the other five methods (Methods II, III,
IV, V, and VI). When examining these five genes, we found that the differences
between the treatment and the control are not obvious. The statistical
significance may be caused entirely by the outliers. Therefore, the five
methods (II, III, IV, V, and VI) are perhaps too sensitive to outliers. For
example, the first plot (5203) of the last row in Figure 3 is the spot that
represents β-globin in the
dataset. We cannot tell much difference between the control and the treatment.
Genes with numbers 2375, 2377, 2379, and 2384 are all “blank” spots. The
observed differences between the control and the treatment are purely caused by
chance, yet these blank spots were detected by Methods II, III, IV, V, or VI.
Finally, Method III [14] and Method V [13] missed gene SR-BI at
array element number 3, which is known to have altered the expression between
the control and the SR-BI transgenic mice.Expression patterns of some genes from the mice data. The first two rows (representing ten genes) are genes detected by
all methods. The third row (five genes) are genes detected by none of the
methods. The last row (five genes) are those detected by Methods II–VI but
not by Methods I.When we examined the plots of the ten genes detected
by the Bayesian mixture model analysis, we found that seven of them have
increased the expression level by the treatment while three of them have
decreased the expression level by the treatment. Interestingly, our method
assigned the three genes with negative regression coefficients into the cluster
with positive regression coefficients. If we look at the estimated parameters
(Table3) again, we realized that both Clusters 1 and 2 may be considered as
the neutral cluster (mean close to zero and variance very small). The third
cluster contains both the up-and downregulated genes because it has a much
larger variance than the other two clusters. The same phenomenon also occurred
for Method VI [15], but the issue was not discussed by
Gusnanto et al. [15]. A possible explanation is that there are too few genes with
negative regression coefficients, which made them hard to form a single cluster
by themselves. However, why were these three genes assigned into Cluster 3 (β>0) instead of
Cluster 2 (β=0), which is
closer to the three genes in terms of the cluster mean? The reason may be that
Cluster 3 has a larger variance than Cluster 2.A separate simulation experiment has been conducted to
verify the notion thatν plays a more
important role than β in calculating
the posterior probability of cluster assignment πik (data not
shown).


## 4. Discussion

Similar to Method IV [16], our method is
based on a hierarchical model in which the parameters of interest (gene
effects) are further described by a mixture distribution. Clustering is made
based on the parameter rather than on a summary statistic such as the t-like
statistical score. This allows us to incorporate the error structure of the
gene expression profile into the linear model (see (10)), and thus capture the most information from the
data. This explains why our method is different from (or even better than) the
other four methods (Methods II, III, V, and VI). But why is our method better
than (or different from) Method IV? This may be contributed by the two differences
between Methods I and IV. One difference is that the incorporated normalization
scheme in our model allows us to deal with only the effect of differential
expression (regression coefficient) whereas Method IV deals with effects of
gene expressions in both the control and the treatment. In other words, we are
dealing with a single variable (regression coefficient, the only parameter of
interest) but Method IV deals with two variables (intercept and regression
coefficients). The dimension reduction from two to one and the simplified
Gaussian mixture distribution of our model may largely contribute to the higher
efficiency of our method. The second difference between our method and Method
IV is that we used the probabilities of cluster assignment to classify genes
into three clusters, and genes assigned into the neutral cluster are excluded
from the list of differentially expressed genes, but Method IV goes one step
further by employing an FDR criterion to select a list of differentially
expressed genes. It is obvious that the FDR generated list of genes depends
largely on the subjective cutting rule set by the investigator. We consider
that the extra step of FDR analysis after gene clustering is not only redundant
but also leads to subjective decision for differential expression analysis.Method II [10] sorts genes based on
the P-value for a regularized t-statistic calculated from each gene.
Information from other genes plays no role in calculating the P-value
for the current gene of interest. This, combined with the subjectiveness of the
cutting rule for gene selection, may largely explains the difference between
Method II and the proposed new method.Method III [14] uses a t-statistic as the data point
for cluster analysis. The original gene expression profile is converted into
the t score. Some information may be lost during the conversion because the
method fails to incorporate the error for calculating the t score. In addition,
the t score is the differential gene expression divided by the standard
deviation of the difference. However, the calculated standard deviation is
subject to large error when the number of replicates per condition is small,
which is usually the case in microarray experiments.Similar to Method II, Method V [13] creates a
list of significant genes based on the moderated t-statistic. Extra information
is borrowed, on the basis of the hierarchical model, from the ensemble of genes
which can assist in inference about each gene individually. However, users need
to subjectively specify the mixing proportions before the algorithm is applied.
The authors pointed out that the estimations of mixing proportions are somewhat
unstable and suggested using 0.01 as a universal prior. In real mice data, the
proportion of differentially expressed genes is about 0.002, which is
overestimated by limma. This may explain why limma detected too many false
positives both in simulation study and mice data analysis.Method VI [15] is quite similar to the
proposed method except that the observed differential expression between the
control and the treatment is used for cluster analysis instead of using the
expected differential expression (parameter) as the basis for clustering.
Again, the error structure of the expression profile is not properly
incorporated into the model when the normalization-like procedure is used. This
explains why the proposed new method outperforms Method VI. In addition, Method
VI is based on the EM algorithm, while the proposed method is based on the MCMC
algorithm. The EM algorithm sometimes may stuck in the so-called local optimal,
while the MCMC has reasonable chance to jump out of the “trap.” This may
explain why Method VI failed in the situation where the borders of clusters are
fuzzy with the effects of differentially expressed genes varying within a wide
span, rather than a narrow span (data not shown).In a quantitative trait-associated microarray data
analysis, Jia and Xu [25] classified genes into several clusters based on the
association of gene expression and the phenotypic value of a quantitative
trait. The Bayesian method developed here can be used for such a quantitative
trait-associated microarray data analysis. Recall that, in the differential
expression analysis, the design matrix X is a variable
indicating whether a gene comes from the control or the treatment. To make such
an extension, the design matrix X in the
differential expression analysis is simply replaced by the phenotypic values of
the trait in the quantitative trait-associated microarray analysis. The method
developed here has the following extra features compared to that of Jia and Xu [25]:
(a) Bayesian method implemented via the MCMC algorithm, (b) constraints on the
cluster means, (c) an imbedded normalization step, and (d) a fixed number of
clusters.In differential gene expression analysis, we usually
deal with two conditions, control and treatment, with a purpose of identifying
genes whose expression levels have altered between the two conditions. In
time-course and dose-response [26] microarray experiments, however,
gene expression are measured from multiple conditions. The Bayesian mixture
model analysis developed here may be extended to detect differentially
expressed genes across multiple conditions. To make such an extension, the
dimensions of X and γi in model (2)
need to be changed to reflect the multiple conditions. Let d+1 be the number
of conditions. The modified dimensions of X and γi should be N×d and d×1, respectively. The change of γi from a scaler
into a vector leads to the following consequences: (1) the clusters of down-
and upregulated genes are not explicitly defined, although the neutral cluster
is still defined as that with βk=0, and (2) the variance of γi in the
experiment with two conditions becomes a variance-covariance matrix in the
experiment with multiple conditions. The first problem may be solved via the
following approaches. Let C be the total
number of clusters from which these genes are sampled. Let the first cluster be
the one consisting of all the neutral genes and thus β1=0 is enforced for
this cluster. The means of the remaining clusters β2−βC are not
restricted, that is, they are sampled freely from their posterior distributions
(multivariate normal). The number of clusters C are found based
on the value of Bayesian factor or BIC. Genes not classified into the neutral
cluster are claimed to be differentially expressed. The actual number of
neutral genes may not be sensitive to the choice of the total number of clusters.
Therefore, one may simply choose any reasonable number of clusters for the
analysis. This number  may be a function of d, say C=3d. Note that C=3 for a single
dimension (one regression coefficient) as done in this study. For two regression
coefficients, C=3×3=9. Therefore, for d regression
coefficients, the number of clusters becomes C=3d. The second problem may be tackled as follows. The
scaled inverse chi-square distribution chosen for the variance may be replaced
by the generalization of the scaled inverse chi-square distribution for the
variance-covariance matrix, called the inverse-Wishart distribution
[21]. This prior distribution for the variance-covariance matrix is
conjugate, and thus the standard sampling algorithm for a random vector from an
inverse-Wishart distribution applies.The interceptαi of gene i in model (2) is
irrelevant to the differential expression analysis and thus it is removed from
the model via a special normalization process, called linear contrasting. In
fact, all effects not related to differential expression can be removed via
such a normalization process. Deleting these irrelevant effects (e.g., dye
effect, block effect, etc.) may avoid tedious estimating procedure for
unspecified parameters used in Zhang et al. [1] and save substantial time in
calculation. To remove these effects, matrices W and αi need to be
customized to reflect the general nature of the normalization. Let h be the number
of irrelevant effects. The dimensions of W and αi should be N×h and h×1, respectively. The coefficient matrix for the linear
contrasts L remains the
same as that defined by model (6). However, matrix L* takes the first N−r(W) eigenvectors of
matrix L. This generalized approach of normalization is
conceptually similar to the ANOVA analysis proposed by Kerr et al. [27], where
two steps are involved in the ANOVA. In the first step, αi is estimated
under model Yi=Wαi+ϵi and then Yi*=Yi−Wα^i, the residuals adjusted by the irrelevant effects, is
used as the raw data for further differential expression analysis. In the
second step, Yi* is described by
model Yi*=Xγi+ϵi and then γi is estimated
and tested to make an inference about the significance of gene i. Three issues need to be emphasized for the
comparison: (1) the generalized approach of normalization proposed in this
study removes αi using no
explicit estimate of αi, (2) the reduced degrees of freedom after adjusting
for the irrelevant effects are used for the new method, and (3) appropriate
covariance structure for the residuals is used for the new method.


## 5. Conclusions

In this paper,
we propose a Bayesian mixture model approach to detect genes that
differentially expressed in control-treatment microarray experiments. Genes are
assigned into one of three constrained clusters, corresponding to clusters of
downregulated, neutral, and upregulated genes, respectively. The Bayes method
is implemented via the Markov chain Monte Carlo (MCMC) algorithm. Genes that
have been assigned into nonneutral clusters are the target genes which we would
like to disclose. Using simulated data as well as data from real microarray
experiments, we demonstrate that the new method outperforms the methods
commonly used in differential expression analysis. Although the new method was
demonstrated using data generated from laboratory animals, it is able to
generalize to genome studies for plants.



---

*Source: 892927-2008-03-25.xml*