A simple implementation of a normal mixture approachto differential gene expression in multiclass microarraysG.J. McLachlan1,2,Ã, R.W. Bean2 and L. Ben-Tovim Jones21Department of Mathematics, University of Queensland and 2ARC Centre in Bioinformatics, Institute forMolecular Bioscience, University of Queensland, St Lucia, Brisbane 4072, Australia
Received on January 25, 2006; revised on April 10, 2006; accepted on April 12, 2006Advance Access publication April 21, 2006Associate Editor: Martin Bishop
suggested in the case where the null density is taken to be standard
Motivation: An important problem in microarray experiments is the
normal (the theoretical null distribution). We also consider the
detection of genes that are differentially expressed in a given number
provision of an initial partition of the genes into two groups for
of classes. We provide a straightforward and easily implemented
the application of the EM algorithm in the case where the adoption
method for estimating the posterior probability that an individual
of the theoretical null distribution would appear not to be appro-
gene is null. The problem can be expressed in a two-component mixture
priate and an empirical null distribution needs to be used. We
framework, using an empirical Bayes approach. Current methods of
demonstrate our approach on three datasets that have been analyzed
implementing this approach either have some limitations due to
previously in the bioinformatics literature.
the minimal assumptions made or with more specific assumptionsare computationally intensive.
Results: By converting to a z-score the value of the test statistic used to
test the significance of each gene, we propose a simple two-componentnormal mixture that models adequately the distribution of this score.
Although biological experiments vary considerably in their design,
The usefulness of our approach is demonstrated on three real datasets.
the data generated by microarrays can be viewed as a matrix of
Availability: An R-program for implementing the approach is freely
expression levels. For M microarray experiments (corresponding to
M tissue samples), where we measure the expression levels of N
genes in each experiment, the results can be represented by the
N · M matrix. Typically, M is no more than 100 (usually much
less in the present context), while the number of genes N is of theorder of 104. The M tissue samples on the N available genes areclassified with respect to g different classes, and it is assumed that
the (logged) expression levels have been preprocessed with adjust-ment for array effects.
Often the first step, and indeed the major goal for many microarraystudies, is the detection of genes that are differentially expressed in a
known number of classes, C1, . . . , Cg. Statistical significance ofdifferential expression can be tested by performing a test for
Differential expression of a gene means that the (class-conditional)
each gene. When many hypotheses are tested, the probability
distribution of its expression levels is not the same for all g classes.
that a type I error (a false positive error) is committed increases
These distributions can differ in any possible way, but the statistics
sharply with the number of hypotheses. In this paper, we focus on
usually adopted are designed to be sensitive to primarily a differ-
the use of a two-component mixture model to handle the multiplic-
ence in the means; for example, the one-way analysis of variance
ity issue. This model is becoming more widely adopted in the
(ANOVA) F-statistic. Even so, the gene hypotheses being tested are
context of microarrays, where one component density corresponds
of equality of distributions across the g classes, which allows the use
to that of the test statistic for genes that are not differentially
of permutation methods to estimate P-values if necessary.
expressed, and the other component density to that of the test statis-
In the special case of g ¼ 2 classes, the one-way ANOVA
tic for genes that are differentially expressed. For the adopted test
F-statistic reduces to the square of the classical (pooled) t-statistic.
statistic, we propose that its values be transformed to z-scores,
Various refinements of the t-statistic have been suggested; see, for
whose null and non-null distributions can be represented by a single
example, the procedure of Tusher et al. (2001).
normal each. We show how this two-component normal mixturemodel can be fitted very quickly via the EM algorithm started from a
point that is completely determined by an initial specification of
Posterior probability of non-differential
the proportion p0 of genes that are not differentially expressed.
A procedure for determining suitable initial values for p0 is
In this paper, we focus on a decision–theoretic approach to the
ÃTo whom correspondence should be addressed.
problem of finding genes that are differentially expressed. We
Ó The Author 2006. Published by Oxford University Press. All rights reserved. For Permissions, please email: [email protected]
Normal mixture approach to differential gene
use a prediction rule approach based on a two-component mixture
can be expressed in terms of the estimated posterior probabilities
model as formulated in Lee et al. (2000) and Efron et al. (2001). We
let G denote the population of genes under consideration. It can bedecomposed into two groups G
p0 0ðwjÞ= ^fðwjÞ ðj ¼ 1‚ . . . ‚ NÞ
genes that are not differentially expressed, and G1 is the comple-ment of G
is the estimated posterior probability that the j-th gene is not dif-
0; that is, G1 contains the genes that are differentially
ferentially expressed. An optimal ranking of the genes can therefore
i denote the prior probability of a gene belong-
be obtained by ranking the genes according to the ^
i (i ¼ 0, 1), and assume that the common density of the test
from smallest to largest. A short list of genes can be obtained by
j for a gene j in Gi is fi(wj). The unconditional density of
j is then given by the two-component mixture model,
taking the top No genes in the ranked list.
Using Bayes Theorem, the posterior probability that the jth gene
is not differentially expressed (i.e. belongs to G0) is given by
In this framework, the gene-specific posterior probabilities provide
Then McLachlan et al. (2004) have proposed that the false discov-
the basis for optimal statistical inference about differential expres-
ery rate (FDR) of Benjamini–Hochberg (1995) can be estimated as
sion. The posterior probability t0(wj) has been termed the local false
discovery rate (local FDR) by Efron and Tibshirani (2002). It quan-
tifies the gene-specific evidence for each gene. As noted by Efron
(2004), it can be viewed as an empirical Bayes version of theBenjamini–Hochberg (1995) methodology, using densities rather
where Nr is the number of selected genes and IA(x) is the indicator
function, which is one if x 2 A and is zero otherwise.
It can be seen from (2) that in order to use this posterior proba-
Similarly, the false non-discovery rate (FNDR) can be estimated by
bility of non-differential expression in practice, we need to be able
0, the mixture density f(wj) and the null density f0(wj),
or equivalently, the ratio of densities f0(wj)/f(wj). Efron et al. (2001)
have developed a simple empirical Bayes approach to this problemwith minimal assumptions. This problem has been studied since
We can also estimate the false positive rate (FPR), e01, and the
under more specific assumptions, including works by Newton et al.
false negative rate (FNR), e10, in a similar manner to give
(2001, 2004), Lo¨nnstedt and Speed (2002), Pan et al. (2003), Zhao
and Pan (2003), Broe¨t et al. (2004), Newton et al. (2004), Smyth
(2004), Do et al. (2005) and Gottardo et al. (2006), among many
others. The fully parametric methods that have been proposed arecomputationally intensive.
01 and e10 denote the two errors when a rule is used to assign a
gene as being differentially expressed or not, where e01 is the proba-
bility of a false positive and e10 is the probability of a false negative.
When controlling the FDR, it is important to have a guide to the
That is, the sensitivity is 1 À e10 and the specificity is 1 À e01. The
value of the associated FNR in particular, as setting the FDR too low
may result in too many false negatives in situations where the genes
of interest (related to biological pathway or target drug) are notnecessarily the top ranked genes; see, for example, Pawitan et al.
where (1 À c) is the cost of a false positive. As the risk depends only
(2005). The local FDR in the form of the posterior probability of
on the ratio of the costs of misallocation, they have been scaled to
non-differential expression of a gene has an advantage over the
add to one without loss of generality.
global measure of FDR in interpreting the data for an individual
The Bayes rule, which is the rule that minimizes the risk (3),
gene; see more details in Efron (2005b) and also Supplementary
assigns a gene to G1 if t0(wj) c; otherwise, the j-th gene is
In practice, we do not know the prior probability p0 nor the densitiesf
0(wj) and f(wj), which will have to be estimated. We shall shortly
denote the test statistic for the test of the null
discuss a simple and quick approach to this estimation problem. If
p0, ^f0(wj) and ^f1(wj) denote estimates of p0, f0(wj) and f1(wj),respectively, the gene-specific summaries of differential expression
Hj : j-th gene is not differentially expressed:
For example, as discussed above, Wj might be the t- or F-statistic,
the j-th gene is to be represented by the two-component normal
depending on whether there are two or multiple classes. Whatever
the test statistic, we proceed in a similar manner as in Efron (2004)
and transform the observed value of the test statistic to a z-score
where p1 ¼ 1 À p0. In (13), f0(zj) ¼ f(zj; 0, 1) is the (theoretical)null density of Zj, where f(z; m, s2) denotes the normal density with
mean m and variance s2, and f1(zj) is the non-null density of Zj. Itcan be approximated with arbitrary accuracy by taking q sufficiently
where Pj is the P-value for the value wj of the original test statistic
large in the normal mixture representation
Wj and F is the N(0, 1) distribution function. Thus
where F0 is the null distribution of Wj. If F0 is the true null distri-
For the datasets that we have analyzed, it has been sufficient to use
bution, then the null distribution of the test statistic Zj corresponding
just a single normal component (q ¼ 1) in (14). In such cases, we can
to zj is exactly standard normal. With this definition of zj, departures
from the null are indicated by large positive values of zj. Ourtransformation (11) is slightly different to that in Efron (2004),
f ðzjÞ ¼ p0fðzj; 0‚ 1Þ + p1fðzj; m ‚s2Þ:
as we wish that only large positive values of the z-score be con-
As pointed out in a series of papers by Efron (2004, 2005a, b), for
sistent with the alternative hypothesis; that is, we want the latter to
some microarray datasets the normal scores do not appear to have
be (upper) one-sided so that the non-null distribution of the z-score
the theoretical null distribution, which is the standard normal. In this
can be represented by a single normal distribution rather than a
case, Efron has considered the estimation of the actual null distri-
mixture in equal proportions of two normal components with
bution called the empirical null as distinct from the theoretical null.
means of opposite sign. Previously, Allison et al. (2002) had con-
As explained in Efron (2005b), the two-component mixture model
sidered mixture modelling of the P-values directly in terms of a
(1) assumes two classes, null and non-null, whereas in reality the
mixture of beta distributions with the uniform (0, 1) distribution (a
differences between the genes range smoothly from zero or near
special form of a beta distribution) as the null component. Pounds
and Morris (2003) considered a less flexible beta mixture model for
In the case where the theoretical null distribution does not appear
the P-values, being a mixture of a uniform (0, 1) distribution for
to be valid and the use of an empirical null distribution would seem
the null and a single beta distribution for the non-null component.
appropriate, we shall adopt the two-component mixture model
In the work of Broe¨t et al. (2004), they used a transformation similar
obtained by replacing the standard normal density by a normal
to the approximation of Wilson and Hilferty (1931) for the x2 dis-
with mean m0 and variance s2 to be inferred from the data. That
tribution to transform the value Fj for the F-statistic for the j-th gene
is, the density of the zj-score is modelled as
In the sequel, we shall model the density of the zj-score by (16). In
In cases where we are unwilling to assume the null distribution F0 of
the case of the theoretical N(0, 1) null being adopted, we shall set
the original test statistic Wj for use in our normal transformation (11),
we can obtain an assessment of the P-value Pj via permutationmethods. We can use just permutations of the class labels for the
gene-specific statistic Wj. This suffers from a granularity problem,
since it estimates the P-value with a resolution of only 1/B, where B isthe number of the permutations. Hence it is common to pool over all
We now consider the fitting of the two-component mixture model
N genes; see Supplementary information for details. The drawback
(15) to the zj, firstly with the theoretical N(0, 1) null adopted. In
of pooling the null statistics across the genes to assess the null dis-
order to fit the two-component normal mixture (15), we need to be
j is that one is using different distributions unless all the
lihood (ML) via the EM algorithm of Dempster et al. (1977), using
j are true. The distribution of the null values of the
differentially expressed genes is different from that of the truly null
the EMMIX program as described in McLachlan and Peel (2000);
genes, and so the tails of the true null distribution of the test statistic is
see also McLachlan and Krishnan (1997). To provide a suitable
overestimated, leading to conservative inferences; see, for example,
starting value for the EM algorithm in this task, we note that the ML
Pan (2003), Guo and Pan (2005) and Xie et al. (2005).
estimate of the parameters in a two-component mixture model sat-isfies the moment equations obtained by equating the sample meanand variance of the mixture to their population counterparts, which
We now proceed to show that by working in terms of the z
as defined by (11), we can provide a parametric version of the
two-component mixture model (1) that is easy to fit. The density
j corresponding to the use of the z-score (11) for
Normal mixture approach to differential gene
and on substituting for them in (17) and (18), we obtain
Hence with the specification of an initial value p
values for the other parameters to be estimated, m1 and s2, are
automatically obtained from (19) and (20). If there is aproblem in so finding a suitable solution for mð0Þ and sð0Þ2 , it
gives a clue that perhaps the theoretical null is inappropriate andthat consideration should be given to the use of an empirical null, asto be discussed shortly.
Following the approach of Storey and Tibshirani (2003) to the
estimation of p0, we can obtain an initial estimate pð0Þ for use in
for an appropriate value of j. There is an inherent bias-variancetrade-off in the choice of j. In most cases, as j grows smaller, thebias of pð0ÞðjÞ grows larger, but the variance becomes smaller.
In this case, we do not assume that the mean m0 and variance s2 of
Fig. 1. Breast cancer data: plot of fitted two-component normal mixture
the null distribution are zero and one, respectively, but rather they
model with theoretical N(0, 1) null and non-null components (weighted
are estimated in addition to the other parameters p
p0)) imposed on histogram of z-scores.
0, we let n0 be the greatest integer less than
of pð0Þ, as obtained from (21). For example, using (21) for j ¼ 0 and
class corresponding to the null component and the remaining N À n
À0.675 led to the initial values of 0.70 and 0.66 for p . The fit we
to the other class corresponding to the alternative component. We
obtained (corresponding to the largest local maximum) is given by
then obtain initial values for the mean and variances of the null and
alternative components by taking them equal to the means and
fitted mixture density superimposed on the histogram of zj-scores,
variances of the corresponding classes so formed. The two-
along with its two components, the theoretical N(0, 1) null density
component mixture model is then run from these starting values
and the N(1.49, 0.94) non-null density weighted by their prior
for the parameters. For the real datasets to be considered next, the
fitting of the normal mixture model takes <0.1 s. The calculation of
component normal mixture model gives a good fit to the empirical
the t-statistics in the first instance for the individual genes takes $4 s
for the Hedenfalk dataset and 2 s for the other two sets on a Pentium
In Table 1, we have listed the FDR estimated from (6) for various
levels of the threshold co in (5). It can be seen, for example, that if cois set equal to 0.1, then the estimated FDR is 0.06 and Nr ¼ 143genes would be declared to be differentially expressed. It is not
suggested that the FDR should be controlled to be around 0.05. It is
For our first example, we consider some data from the study of
just that in this example, its control at this approximate level yields a
Hedenfalk et al. (2001), which examined gene expressions in breast
number (143) of differentially expressed genes that is not too
cancer tissues from women who were carriers of the hereditary
unwieldy for a biologist to handle in subsequent confirmatory
BRCA1 or BRCA2 gene mutations, predisposing to breast cancer.
experiments; the choice of co is discussed in Efron (2005b).
The dataset comprised the measurement of N ¼ 3226 genes using
In the original paper, Hedenfalk et al. (2001) selected 176 genes
based on a modified F-test, with a p-value cut off of 0.001. Com-
tumours. We column normalized the logged expression values,
paring genes which were selected in our set of 143, we found 107 in
and ran our analysis with the aim of finding differentially expressed
common, including genes involved in DNA repair and cell death,
genes between the tumours associated with the different mutations.
which are over-expressed in BRCA1-mutation-positive tumours,
As in Efron (2004), we adopted the classical pooled t-statistic as our
such as MSH2 (DNA repair) and PDCD5 (induction of apoptosis).
Storey and Tibshirani (2003) in their analysis of this dataset,
j for each gene j and we used the t-distribution
selected 160 genes by thresholding genes with q-values less than
or equal to a ¼ 0.05 (an arbitrary cut-off value), of which there are
j in the computation of the P-value Pj from (12).
We fitted the two-component normal mixture model (15) with the
113 in common with our set of 143. Overall, 101 genes were
standard normal N(0, 1) as the theoretical null, using various values
selected in common to all three studies, with 24 genes unique to
Table 1. Estimated FDR and other error rates for various levels of thethreshold co applied to the posterior probability of non-differential expressionfor the breast cancer data, where Nr is the number of selected genes
our set. We searched publicly available databases for the biologicalfunctions of these genes, and found these included DNA repair, cellcycle control and cell death, suggesting good evidence for inclusionof these genes (see the Supplementary information for a fullerdescription).
Concerning the other type of allocation rates for the choice of
co ¼ 0.1 in (5), the estimates of the FNDR, FNR and FPR are equalto 0.32, 0.88 and 0.004, respectively. The FNR of 0.88 means thatthere would be quite a few false negatives among the genes declaredto be null (not differentially expressed).
Among other analyses of this dataset, p0 was estimated to be 0.52
by Broe¨t et al. (2004), 0.64 by Gottardo et al. (2006), 0.61 by Ploner
Fig. 2. Colon cancer data: plot of fitted two-component normal mixture
et al. (2006) and 0.47 by Storey (2002). In the fully parametric
model with theoretical N(0, 1) null and non-null components (weighted re-
Bayesian approach of Broe¨t et al. (2004), the mean of the null
p0)) imposed on histogram of z-scores.
component was fixed at zero, but the variance was allowed to befree during the estimation process for computational convenience.
In Ploner et al. (2006), 56 genes with highly extreme expressionvalues were first removed as in Storey and Tibshirani (2003).
We apply our method next to the well-known study of Alon et al.
We also considered the fitting of the two-component normal
(1999), where Affymetrix arrays were used to compare the gene
mixture model (16) with the null component mean and variance,
expressions from colon cancer tissues with normal colon tissues.
They measured expressions of over 6500 genes in n
1 and s2 . We found that this fit from using the
2 ¼ 22 normal colon tissue samples. The samples were taken
empirical null in place of the N(0, 1) theoretical null is similar
from 40 different patients, so that 22 patients supplied both a normal
to the fit in Figure 1; see Figure 1 in Supplementary information.
and a tumor tissue sample. We consider the set of N ¼ 2000 genes
Efron (2003) writes that ‘there is ample reason to distrust the theo-
with highest minimal intensity across the samples as in Alon et al.
retical null’ in the case of the Hedenfalk data. The difference in our
findings may be due to the fact that our gene expression data seems
As in the last example, we adopted the (pooled) t-statistic as our
to differ when compared with the expression data presented in Efron
test statistic Wj for each gene j, and so we used the t-distribution
and Tibshirani (2002). In other analyses of this dataset, Newton
function with 60 degrees of freedom, F60, as the null distribution of
et al. (2001), Tusher et al. (2001) and Gottardo et al. (2006) con-
Wj in the computation of the P-value Pj from (12). We fitted the
cluded that there were 375, 374 and 291 genes, respectively, dif-
two-component normal mixture model (15) with the theoretical
ferentially expressed when the FDR is controlled at the 10% level. It
N(0, 1) null to the zj-scores, which gave the fit, ^
can be seen from Table 1 that our approach gives 338 genes as being
2:21. A plot of the normal mixture density
differentially expressed when the threshold of c
and its two components is displayed in Figure 2. Our estimate of
on the posterior probability of non-differential expression for which
0.39 for p0 coincides with the empirical Bayes estimate reported in
the implied FDR is 11% and the FNR is 73%. The corresponding
values with the use of the empirical null are 12 and 76% for the FDR
If we again consider a threshold of co ¼ 0.1 in (5), then 433 genes
and FNR, respectively, with 235 genes declared to be differentially
would be declared to be differentially expressed with an implied
expressed. According to the Bayesian Information criterion (BIC),
FDR and FNR of 3 and 65%, respectively. The smooth muscle
the empirical null would not be selected in favor of the theoretical
genes provide a positive control for differential gene expression,
N(0, 1) null (see Supplementary information). The (main) reason for
as the normal tissues were later found to have included smooth
fewer genes being declared differentially expressed with the use of
muscle tissue in the biopsy samples. We find that for the smooth
the empirical than with the theoretical null is that the estimate of
muscle genes (J02854, T60155, M63391, D31885, X74295,
Normal mixture approach to differential gene
non-differential expression that is <0.0015, and these genes arefound within our top 66 ranked genes. For this dataset, Do et al. (2005) concluded that each gene in this cluster has an estimatedposterior probability (of non-differential expression) of 0.002or less.
In their original analysis, Alon et al. (1999) identified a ribosomal
gene cluster, associated with over-expression in tumor tissues rel-ative to normal tissues. Of these 29 genes (as listed in Table 1 ofAlon et al. 1999), only 10 are declared differentially expressed at athreshold of co of 0.05 (within a subset of 296 genes). Do et al. (2005) similarly identified only 13 genes with posterior probabilitiesof non-differential expression of <0.05. Their two genes with weak-est discriminatory ability were H77302 and T63484, with posteriorprobabilities of non-differential expression of 0.44 and 0.49; oursare H77302 and R85464, with posterior probabilities of 0.59 and0.65. Our results, supported by the findings in Do et al. (2005),suggest that only a minority of the ribosomal genes are actuallydifferentially expressed between the tumor and normal classes. Wealso considered the use of the empirical null instead of the theo-retical N(0, 1) null for this dataset. The plot of the normal mixturedensity is given in Figure 2 in the Supplementary information. For athreshold of co ¼ 0.1 in (5), 253 genes would be declared to bedifferentially expressed with an implied FDR and FNR of 4 and74%, respectively. The reason for fewer genes being declared dif-ferentially expressed with the use of the empirical than with thetheoretical null is that the estimate of p0 is greater ( ^
Fig. 3. HIV data: plot of fitted two-component normal mixture model with
According to the BIC, the empirical null would not be selected
empirical null and non-null components (weighted respectively by ( ^
in favor of the theoretical N(0, 1) null (see Supplementary
p0) imposed on histogram of z-scores.
considered as a means for providing improved versions of the theo-retical null rather than empirical nulls.
We considered an empirical null for this dataset by fitting the
We consider finally the HIV dataset of van’t Wout et al. (2003), and
normal mixture model (16), obtaining as our fit, ^
as discussed in Gottardo et al. (2005). van’t Wout et al. used cDNA
arrays to measure expression levels of N ¼ 7680 genes in CD4-
T-cell lines, at time t ¼ 24 h after infection with the HIV-1 virus.
obtained by Efron (2005b) for this dataset. Using a fully Bayesian
There were four control slides (pooled mRNA from uninfected
approach to model the gene expressions, Gottardo et al. (2006)
cells) and four slides run using pooled mRNA from the infected
estimated p0 to be 0.993, which is similar to our quick and easy
1 ¼ 4 and n2 ¼ 4. We column normalized the logged
expression values, using the dataset as downloaded from
On setting a threshold of co ¼ 0.01 on the posterior probability of
hajek.stat.ubc.ca/$raph. The dataset contains 12 HIV-1 genes
non-differential expression, we find that there is a subset of 15 genes
which were used as positive controls, as they are known to be
that are differentially expressed, which contains the 12 genes known
differentially expressed in infected versus uninfected cells. As in
to be differentially expressed from external information. The
Efron (2004), we adopted the (pooled) t-statistic and used the
implied FDR is 0.002. For co ¼ 0.1, we obtain that 37 genes are
t distribution with 6 degrees of freedom for F
differentially expressed with an implied FDR of 0.03.
The HIV dataset is unique in the examples discussed, as cell lines
j from (12). As noted by Efron (2004, 2005b), this
dataset is an example of one where the theoretical N(0, 1) is inap-
were used rather than tissue biopsies from patients. Additionally,
plicable due to the underdispersion of the z
each of the four microarray experiments within the class (uninfected
(Fig. 3). A fit of a single normal to all the z
or infected) used the same pooled RNA, so that each of these were
N(À0.16, 1.06) distribution, clearly showing that it is not appropri-
technical replicates rather than individual samples. These factors
ate to adopt the theoretical N(0, 1) null component in our
may explain the decreased dispersion.
two-component mixture model (15). As there are only 6 degreesavailable for the use of the parametric t assumption for the nulldistribution F
0 of the test statistic, we also computed the P-values
using the permutation estimate with B ¼ 35 permutations (pooled
In this paper, we consider the problem of detecting which genes are
over the genes). It gave a similar fit, N(À0.146, 0.997), for a single
differentially expressed in multiple classes of tissue samples, where
normal fitted to the zj-scores. As explained by Efron (2005b),
the classes represent various clinical or experimental conditions.
permutation methods in the context of this problem should be
The available data consist of the expression levels of typically a
very large number of genes for a limited number of tissues in each
model the distribution of the P-values in their analyses. However, it
class. Usually, a test statistic such as the classical t in the case of two
is advantageous to work as proposed here in terms of the zj-scores,
classes or the F in case of multiple classes is formed for a test of
which can be modelled by normal components on the real line rather
equality of the class means. The key step in this approach is to
than working in terms of the P-values.
transform the observed value of the test statistic for each gene j to a
Finally, we should mention explicitly that it has been assumed
z-score zj by using the inverse standard normal distribution function
throughout that the genes are all independently distributed. Typi-
of the implied P-value Pj, similar to its use in Efron (2004) and his
cally in practice, this independence assumption will not hold for all
subsequent papers on this problem. We demonstrate that a two-
the genes. As cautioned by Qiu et al. (2005), care is needed in
component normal mixture model is adequate for modelling the
extrapolating results valid in the case of independence to dependent
empirical distribution of the zj-scores, where the first component is
gene data. In the Supplementary information, we report the results
the standard normal, corresponding to the null distribution of the
of a few initial simulations that we have performed to investigate the
score, and the second component is a normal density with unspe-
effect of correlation between some of the genes on the distribution
cified (positive) mean and variance, corresponding to the non-null
of the z-scores. For moderate correlation, the effect was small, while
distribution of the score. This model can be used to provide a
for strong correlation it was found that the use of the empirical null
straightforward and easily implemented assessment of whether a
in place of the theoretical N(0, 1) null avoided any marked effect
gene is null (not differentially expressed) in terms of its posterior
(see Figs 3 and 4 in the Supplementary information). Another
probability of being a null gene. Estimates of this posterior proba-
approach to handle the effect of strong correlation between the
bility can be easily obtained by using the EM algorithm to fit the
genes would be first to cluster the (normalized) gene profiles on
two-component normal mixture model via ML. As there are mul-
the basis of Euclidean distance, so that highly correlated genes are
tiple local maximizers, consideration has to be given to the choice of
put in the same cluster. The genes in each of those clusters with
starting values for the algorithm. We show that the specification of
highly correlated members can then be represented by a single gene
or a linear combination of the member genes (a metagene). To
specifies a starting point for the fitting of the normal mixture model
incorporate a priori knowledge that some genes belong to the
with the theoretical choice of N(0, 1) as the null component. An
same pathway, we can take the cluster labels for these genes to
can be tried, and a guide to its endpoints is
be the same in the specification of the complete data in the EM
given by values of p0 obtained by equating the number of zj values
less than a threshold j to the expected number under the theoretical
Conflict of Interest: none declared.
N(0, 1) null component. We consider too the case where the theo-retical N(0, 1) null is not tenable and an empirical null is adoptedwith mean and variance estimated from the data. Also, the estima-
tion of the FDR and its control are considered, along with the
Allison,D.B. et al. (2002) A mixture model approach for the analysis of microarray
estimation of other relevant rates such as the FNR. Note that it
gene expression data. Comput. Statist. Data Anal., 39, 1–20.
is not valid to make claims as to the relative superiority of the
Alon,U. et al. (1999) Broad patterns of gene expression revealed by clustering analysis
two models corresponding to the theoretical and empirical nulls
of tumor and normal colon tissues probed by oligonucleotide arrays. Proc. NatlAcad. Sci. USA, 96, 6745–6750.
on the basis of these error rates, as they are only valid for the
Benjamini,Y. and Hochberg,Y. (1995) Controlling the false discovery rate: a practical
model under which they were calculated. Our approach is demon-
and powerful approach to multiple testing. J. R. Stat. Soc. B, 57, 289–300.
Broe¨t,P. et al. (2004) A mixture model-based strategy for selecting sets of genes in
Concerning the choice between the use of the theoretical N(0, 1)
multiclass response microarray experiments. Bioinformatics, 20, 2562–2571.
Dempster,A.P. et al. (1977) Maximum likelihood from incomplete data via the EM
null and an empirical null, the intent in the first instance is to use the
algorithm (with discussion). J. R. Stat. Soc. B, 39, 1–38.
former in modelling the density of the zj-scores. In some situations
Do,K-A. et al. (2005) A Bayesian mixture model for differential gene expression. Appl.
as with the HIV dataset above, it will be clear that the use of the
theoretical null is inappropriate. In other situations, an informed
Efron,B. (2004) Large-scale simultaneous hypothesis testing: the choice of a null
choice between the theoretical and empirical null components can
hypothesis. J. Am. Stat. Assoc., 99, 96–104.
Efron,B. (2005a) Selection and estimation for large-scale simultaneous inference.
be made on the basis of the increase in the log likelihood due to the
Technical Report. Department of Statistics, Stanford University, Stanford, CA,
use of an empirical null with its two extra parameters. For this
purpose we can use BIC as in the first two examples above.
Efron,B. (2005b) Local false discovery rates. Technical Report. Department of Statis-
The reliability of our approach obviously depends on how well
tics, Stanford University, Stanford, CA, False.pdf
the proposed two-component normal mixture model approximates
Efron,B. and Tibshirani,R. (2002) Empirical Bayes methods and false discovery rates
the empirical distribution of the zj-scores. In the three examples here
for microarrays. Genet. Epidemiol., 23, 70–86.
and in our analyses of other datasets not presented here, this normal
Efron,B. et al. (2001) Empirical Bayes analysis of a microarray experiment. J. Am. Stat.
mixture model provides an excellent approximation. Its fit can be
assessed either by visual inspection of a plot of the fitted normal
Gottardo,R. et al. (2006) Bayesian robust inference for differential gene expression in
cDNA microarrays with multiple samples. Biometrics, 62, 10–18.
mixture density versus a histogram of the zj-scores or, more for-
Guo,X. and Pan,W. (2005) Using weighted permutation scorse to detect differential
mally, by a likelihood ratio test for the need for an additional normal
gene expression with microarray data. J. Bioinformatics Compat. Biol., 3,
density to represent the non-null distribution of the z
similar note on the adequacy of a two-component normal mixture
Hedenfalk,I. et al. (2001) Gene-expression profiles in hereditary breast cancer. N. Engl.
model, Pounds and Morris (2003) found that a two-component
Lee,M.-L.T. et al. (2000) Importance of replication in microarray gene expression
mixture of the uniform (0,1) distribution and a single beta compo-
studies: statistical methods and evidence from repetitive cDNA hybridizations.
nent (with one unspecified unknown parameter) was adequate to
Proc. Natl Acad. Sci. USA, 97, 9834–9838.
Normal mixture approach to differential gene
Lo¨nnstedt,I. and Speed,T. (2002) Replicated microarray data. Statist. Sinica, 12, 31–46.
Qiu,X. et al. (2005) Correlation between gene expression levels and limitations of the
McLachlan,G.J. and Krishnan,T. (1997) The EM Algorithm and Extensions. Wiley,
empirical Bayes methodology for finding differentially expressed genes. Stat. Appl.
Genet. Mol. Biol., 4, No. 1, Article 34.
McLachlan,G.J. and Peel,D. (2000) Finite Mixture Models. Wiley, New York.
Smyth,G.K. (2004) Linear models and empirical Bayes methods for assessing differ-
Newton,M.A. et al. (2001) On differential variability of expression ratios: improving
ential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol., 3, No. 1,
statistical inference about gene expression changes from microarray data.
Storey,J.D. (2002) A direct approach to false discovery rates. J. R. Stat. Soc. Ser B, 64,
McLachlan,G.J. et al. (2004) Analyzing Microarray Gene Expression Data. Wiley,
Storey,J.D. and Tibshirani,R. (2003) Statistical significance for genome-wide studies.
Newton,M.A. et al. (2004) Detecting differential gene expression with a semiparamet-
Proc. Natl Acad. Sci. USA, 100, 9440–9445.
ric hierarchical mixture method. Biostatistics, 5, 155–176.
Tusher,V.G. et al. (2001) Significance analysis of microarrays applied to the ionizing
Pan,W. (2003) On the use of permutation in and the performance of a class of non-
radiation response. Proc. Natl Acad. Sci. USA, 98, 5116–5121.
parametric methods to detect differential gene expression. Bioinformatics, 19,
van’t Wout,A.B. et al. (2003) Cellular gene expression upon human immuno-
deficiency virus type 1 infection of CD4+-T-cell lines. J. Virol., 77,
Pan,W. et al. (2003) A mixture model approach to detecting differentially expressed
genes with microarray data. Functional and Integrative Genomics, 3, 117–124.
Wilson,E.B. and Hilferty,M.M. (1931) The distribution of chi-square. Proc. Nat. Acad.
Pawitan,Y. et al. (2005) False discovery rate, sensitivity and sample size for microarray
studies. Bioinformatics, 21, 3017–3024.
Xie,Y. et al. (2005) A note on using permutation-based false discovery rate estimates to
Ploner,A. et al. (2006) Multidimensional local false discovery rate for microarray
compare different analysis methods for microarray data. Bioinformatics, 21,
studies. Bioinformatics, 22, 556–565.
Pounds,S. and Morris,S.W. (2003) Estimating the occurrence of false positives and
Zhao,Y. and Pan,W. (2003) Modified nonparametric approaches to detecting differ-
false negatives in microarray studies by approximating and partitioning the empiri-
entially expressed genes in replicated microarray experiments. Bioinformatics, 19,
cal distribution of p-values. Informatics, 19, 1236–1242.
Date : 25-11-2013 Folio n° : Commune de Mellionnec REGISTRE DES DELIBERATIONS DU CONSEIL MUNICIPAL ET DES DECISIONS DU MAIRE SEANCE DU 25 NOVEMBRE 2013 L’an deux mil treize, le vingt-cinq novembre, à dix-huit heures trente. Le Conseil Municipal de cette commune, régulièrement convoqué, le quinze novembre deux mil treize s’est réuni au nombre prescrit par la loi, da
TITLE: A SAFETY GUIDE FOR SPINAL ANESTHESIA AUTHORS: Yasuhiro Shimada (1), Naofumi Ashizawa (2), Yoshio Kato (3), Reiji Wakuzawa (4), Toyohisa Arai (5), Yasuhiko Hashimoto (6), Ken Asayama (7), and Members of the Subcommittee on the Safety of Anesthesia in the Society for the Safety of Medical Practice (1) Department of Anesthesiology, Nagoya University School of Medicine (2) Department of An