Elvumgar.fea.st

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.

Source: http://elvumgar.fea.st/papers/1608.pdf