We proposed a novel method (BE) for analysis of microarray data from experiments where only up- or down-regulated genes were expected or relevant. Its merit originates from an empirical Bayes framework and foundation in the EVT, as well as its good accuracy, precision and stability demonstrated in the simulated datasets. The main limitations of the BE method pertain to its sensitivity to the assumed pDE and the higher computational cost of its underlying numerical approximation through the MC integration method.
The main asset of empirical Bayes methods in analysis of differential expression is their use of information from hundreds or thousands of simultaneously tested genes to support testing at the gene level. However, empirical Bayes methods do make distributional assumptions, particularly when the choice of priors for the parameters is limited to conjugate priors (Lonnstedt, 2001
). We stepped out of the conjugate prior class, and in experiments where only up- or only down-regulated genes are expected or relevant, assumed that the mean ratios of expressions of DE genes follow the EVD. In microarrays expecting a single direction of regulation, the validity of this assumption is based on the fact that the distribution of the true genes’ mean ratios of expressions (the log of it) cannot be symmetric around zero. In the experiments where only a certain direction of regulation is relevant, it is advantageous to use the BE over the BN method even if the gene expression data are more or less symmetric around zero. That is because, being restricted to one direction of expression, due to the design of the experiment or the interest of an investigator, the genes characterized by such expressions are actually ‘selected’ beyond a threshold expression, and so, they show extreme behaviour.
The three families in EVD (Gumbel, Fréchet and Weibull) have distinct forms of tail behaviour, with finite limit distributions in its upper end-point for the Weibull distribution, and infinite ones for the Gumbel and Fréchet distributions (Coles, 2001
). Furthermore, the density in the upper end-point decays exponentially for the Gumbel distribution and polynomially for the Fréchet distribution (Coles, 2001
). While, consequently, the three families give quite different representations of the extreme value behaviour, it may not be obvious which best represents the data and the problem at hand. It is thus convenient to use a single family GEV rather than ‘guessing’ which of the three original EVD families is proper, because the data themselves (through inference about the shape parameter) determine the most appropriate type of tail behaviour (Coles, 2001
To estimate the hyperparameters of the EVD, we applied the MLE procedure. While many techniques have been proposed for estimation of hyperparameters in extreme value models (including graphical techniques and moment-based and likelihood-based methods), the MLE approach is particularly appealing due to its all-around utility and adaptability (Coles, 2001
). Indeed, our analyses showed that BE works well even when hyperparameters a
are estimated from a small sample of genes (~10) related to smaller pDE and/or a smaller number of tested genes G
. However, caution is needed when MLE is used to estimate means and SEs of EVD hyperparameters. Because the end-points of the GEV distributions are functions of the parameter values, regularity conditions required for the usual asymptotic properties associated with the maximum likelihood (ML) estimator are not satisfied by the GEV model (Coles, 2001
). Particularly, in the experiments involving detection of down-regulated genes, if the value of the shape parameter (c
) lies between −1 and −0.5, ML estimators are generally obtainable, but do not have the standard asymptotic properties, while when c
<−1, ML estimators are unlikely to be obtainable (Smith, 1985
). Conversely, when c
> −0.5, ML estimators are regular (have the usual asymptotic properties) (Smith, 1985
). Therefore, in experiments designed to detect up-regulated genes, where c
>0, the usual asymptotic properties immediately apply. For detection of down-regulated genes, it is reasonable to negate the contrast estimators from gene-specific linear models before estimation of the MLE hyperparameters and subsequent approximation of BE statistics, i.e. to switch from modelling minima to modelling maxima.
While only one direction of differential expression is often expected in the wild-type vs
. mutant experiments, it is possible to observe a few genes with their expression in an opposite direction. For example, while the sigma factor σB
is a positive regulator of gene expression, some genes show higher transcript levels in a sigB
null mutant, likely representing indirect effects, which may, however, be biologically relevant. For example, in the bacterium Staphylococcus aureus
, genes encoding a number of exo-enzymes and toxins show higher transcript levels in a sigB
null mutant strain as compared with the parent strain with an intact sigB
gene, suggesting indirect negative regulation of these genes by σB
(Bischoff et al., 2004
). To test differential expression of a suspected biologically relevant gene regulated in the direction opposite to the direction of regulation of the majority of the genes in the experiment, its contrast estimator should be negated prior to applying the BE method. Otherwise, the BE method would rank this gene low, indicating non-differential expression.
Because there is no conjugate structure between the normal likelihood of a gene being DE and the EVD prior, it is impossible to estimate BE analytically. Therefore, the MC integration approximation technique was applied. While this allowed us to base the choice of the prior on the nature of the problem, rather than on the available conjugate priors, it came at the price of a higher, but still reasonable (e.g. 16 min for 50 000 iterations in a dataset with 10 000 genes), computational cost related to a large number of tested genes and (usually) a large number of iterations required to achieve a reasonable accuracy of the BE statistics (because the accuracy increases only as the square root of the number of iterations). The MCSEs of the estimated BE statistics differed for various genes, being the lowest for genes with a BE statistic around zero. That is as expected because genes with a high mean expression ratio, corresponding to the upper tail of the EVD prior, will likely have a high BE statistic. Similarly, genes with a high variance, corresponding to the upper tail of the variance prior, will likely have a low BE statistic. At each iteration, the probability of drawing a sample from a tail is lower than from a body of the prior distribution, so these genes require more iterations to achieve sufficiently low MCSEs. If the BE statistic is used to identify DE genes, as in the present study, an even lower number of iterations may be sufficient. However, when BE is used with the purpose of precise ranking of genes by their strength of differential expression, a larger number of iterations may be necessary.
As discussed by Smyth (2004)
, it is possible to estimate pDE from the data, for example, from posterior odds through MLE or based on the P
-values from t
-statistics (moderated or ordinary) (such as in Langaas et al., 2005
and Wu et al., 2006
). Nevertheless, estimation of pDE is unstable (related to collapse in estimation of posterior odds caused by boundary values of pDE=1 and pDE=0 having positive probability), and may be sensitive to the particular prior distribution assumed for the contrast estimators from linear models and to dependence between the genes (Smyth, 2004
). To bypass these problems, Smyth (2004)
and Lonnstedt & Speed, (2002)
fixed pDE and then estimated the remaining hyperparameter. This same approach has been adopted in the BE method. A wrong choice of pDE would not change the shape of the BE statistic vs
. fold change plot (i.e. ranking of the genes would stay intact), but it would move the BE up and down on the y
-axis (as in Lonnstedt & Speed, 2002
). Also, through influence of the pDE on the hyperparameters of EVD, it would vary fold change corresponding to BE=0. Accepting that BE=0 is a natural cutoff separating DE from non-DE, the fold change corresponding to this point represents the ‘meaningful fold change cutoff’ above which a gene could be considered meaningfully DE. This ‘meaningful fold change cutoff’ could give some indication of whether the chosen pDE is wrong. For example, in experiments expecting only up-regulated genes, the ‘meaningful fold change cutoff’ ≤1 would indicate that the selected pDE is too high.
The accuracy and precision of BN, BL and BE were assessed based on the achieved FDR and FNR over 100 simulated datasets generated under the EVD model. In an experiment, a researcher could choose a tolerable FDR and try to minimize FNR (i.e. maximize power) or choose a tolerable FNR and try to minimize FDR. In the simulated datasets, we mimicked both situations, and in both, the accuracy of the BE method was at least as good as or better than the other two tested methods. In addition to being more accurate, BE was also more precise than BN and BL in the simulated datasets. We also tested the sensitivity of BE to several key characteristics of a microarray experiment, namely, the pDE, the number of arrays (usually representing the number of biological replicates) and the number of technical replicates, as well as its performance in experiments with a high number of tested genes but with only a few DE genes. Overall, BE performed the best in all four simulated model-datasets and under all the tested simulated scenarios. When applied to real datasets, BE performed well, with most of the findings being validated by other independent studies (elaborated in the Appendix
). This ‘ reality check’ is a valuable addition to simulation studies because datasets used in simulation analyses have been replicated from the same model as the BE is built upon, thus limiting ability of the simulation studies to give a completely unbiased evaluation of BE performance.