 Methodology article
 Open Access
 Published:
Empirical Bayes estimation of posterior probabilities of enrichment: A comparative study of five estimators of the local false discovery rate
BMC Bioinformatics volume 14, Article number: 87 (2013)
Abstract
Background
In investigating differentially expressed genes or other selected features, researchers conduct hypothesis tests to determine which biological categories, such as those of the Gene Ontology (GO), are enriched for the selected features. Multiple comparison procedures (MCPs) are commonly used to prevent excessive false positive rates. Traditional MCPs, e.g., the Bonferroni method, go to the opposite extreme: strictly controlling a familywise error rate, resulting in excessive false negative rates. Researchers generally prefer the more balanced approach of instead controlling the false discovery rate (FDR). However, the qvalues that methods of FDR control assign to biological categories tend to be too low to reliably estimate the probability that a biological category is not enriched for the preselected features. Thus, we study an application of the other estimators of that probability, which is called the local FDR (LFDR).
Results
We considered five LFDR estimators for detecting enriched GO terms: a binomialbased estimator (BBE), a maximum likelihood estimator (MLE), a normalized MLE (NMLE), a histogrambased estimator assuming a theoretical null hypothesis (HBE), and a histogrambased estimator assuming an empirical null hypothesis (HBEEN). Since NMLE depends not only on the data but also on the specified value of Π_{0}, the proportion of nonenriched GO terms, it is only advantageous when either Π_{0} is already known with sufficient accuracy or there are data for only 1 GO term. By contrast, the other estimators work without specifying Π_{0} but require data for at least 2 GO terms. Our simulation studies yielded the following summaries of the relative performance of each of those four estimators. HBE and HBEEN produced larger biases for 2, 4, 8, 32, and 100 GO terms than BBE and MLE. BBE has the lowest bias if Π_{0} is 1 and if the number of GO terms is between 2 and 32. The bias of MLE is no worse than that of BBE for 100 GO terms even when the ideal number of components in its underlying mixture model is unknown, but has high bias when the number of GO terms is small compared to the number of estimated parameters. For unknown values of Π_{0}, BBE has the lowest bias for a small number of GO terms (232 GO terms), and MLE has the lowest bias for a medium number of GO terms (100 GO terms).
Conclusions
For enrichment detection, we recommend estimating the LFDR by MLE given at least a medium number of GO terms, by BBE given a small number of GO terms, and by NMLE given either only 1 GO term or precise knowledge of Π_{0}.
Background
The development of microarray techniques and highthroughput genomic, proteomic, and bioinformatics scanning approaches (such as microarray gene expression profiling, mass spectrometry, and ChIPonchip) has enabled researchers to simultaneously study tens of thousands of biological features (e.g., genes, proteins, singlenucleotide polymorphisms [SNPs], etc.), and to identify a set of features for further investigation. However, there remains the challenge of interpreting these features biologically. For a given set of features, the determination of whether some biological information terms are enriched (i.e., differentially represented), compared to the reference feature set, is termed the feature enrichment problem. The biological information term may be, for instance, a Gene Ontology (GO) term [1, 2] or a pathway in the Kyoto Encyclopedia of Genes and Genomes (KEGG) [3]. We call this problem the feature enrichment problem.
This problem has been addressed using a number of highthroughput enrichment tools, including DAVID, MAPPFinder, OntoExpress and GoMiner [47]. Huang et al. [8] reviewed 68 distinct feature enrichment analysis tools. These authors further classified feature enrichment analysis tools into 3 categories: singular enrichment analysis (SEA), gene set enrichment analysis (GSEA), and modular enrichment analysis (MEA). In this article, we propose empirical Bayes solutions to the SEA problem using genes as archetypal features. Without loss of generality, we consider whether some specific biological categories are enriched for differentially expressed genes with respect to the reference genes.
Indeed, like other enrichmentdetection methods, our methods apply much more broadly. They can assess enrichment given any sublist of features selected for future study, not just a list of genes considered differentially expressed. An anonymous referee pointed out these examples of such lists of candidate features that arise in the context of whole genome sequencing:

genes with SNPs

genes with copy number variations

genes with loss of heterozygosity
These examples and those of our first paragraph do not exhaust contemporary applications, and the feature enrichment problem may occur in unforeseen domains of study. Thus, our illustrative use of differential gene expression as a running example should not be interpreted as a limitation.
Existing enrichment tools mainly address the feature enrichment problem using a pvalue obtained from an exact or approximate statistical test (e.g., Fisher’s exact test, the hypergeometric test, binomial test, or the χ^{2} test). For each GO term or other biological category, the null hypothesis tested and its alternative hypothesis are as follows:
Here and in the remainder of the paper, we use GO terms as concrete examples of biological categories without excluding applications of the methods to categories from other relevant databases. The general process begins as follows:

For each GO term, construct Table 1 based on the preselected genes (e.g., differentially expressed (DE) genes) and reference genes (e.g., all genes measured in a microarray experiment).

Compute the pvalue for each GO term using a statistical test that can detect enrichment for the preselected genes.
Multiple comparison procedures (MCPs) are then applied to the resulting pvalues to prevent excessive false positive rates. The false discovery rate (FDR) [9] is frequently used to control the expected proportion of incorrectly rejected null hypotheses in gene enrichment studies [1012] because it has lower false negative rates than Bonferroni correction and other methods of controlling the familywise error rate. Methods of FDR control assign qvalues [13] to biological categories, but qvalues are too low to reliably estimate the probability that the biological category is not enriched for the preselected features. Thus, we study application of better estimators of that probability, which is technically known as the local FDR (LFDR). Hong et al. [14] used an LFDR estimator to solve a GSEA problem and pointed out that this was less biased than the qvalue for estimating the LFDR, the posterior probability that the null hypothesis is true.
Efron [15, 16] devised reliable LFDR estimators for a range of applications in microarray gene expression analysis and other problems of largescale inference. However, whereas microarray gene expression analysis takes into account tens of thousands of genes, the feature enrichment problem typically concerns a much smaller number of GO terms. While these methods are appropriate for microarrayscale inference, they are less reliable for enrichmentscale inference [1719]. Thus, we will specifically adapt LFDR estimators that are appropriate for smallerscale inference to address the SEA problem. Again, we will focus on genes and GO terms for the sake of concreteness. Nevertheless, the estimators used can be applied to other features and to other biological terms (e.g., metabolic pathways).
The sections of this paper are arranged as follows. We first introduce some preliminary concepts in the feature enrichment problem. Next, two previous LFDR estimators and three new LFDR estimators are described. Following this, we compare the LFDR estimators by means of a simulation study and an application to breast cancer data. Finally, we draw conclusions and make recommendations on the basis of our results.
Preliminary concepts
The feature enrichment problem described in the Background section is stated here more formally for the application of LFDR methods in the next section.
Likelihood functions
In Table 1, x_{1} and x_{2} are the observed numbers DE genes and EE genes in a given GO category, respectively. Whereas n is the total number of DE genes, N is the total number of reference genes. Thus, N − n is the total number of EE genes. The columns gives the numbers of DE genes and EE genes, and the rows give the numbers of genes in the GO category and outside the GO category.
Let x_{1} and x_{2}, respectively, denote the random numbers of DE and EE genes in a GO category. The observed values x_{1} and x_{2} are modeled as realizations of x_{1} and x_{2}. x_{1} and x_{2} follow binomial distributions, namely, X_{1} ∼ Binomial(n, π_{1}) and X_{2} ∼ Binomial(N − n,π_{2}), where π_{1} and π_{2} are the probabilities that a gene is DE and EE, respectively, given that it is in the GO category. Under the assumption that x_{1} and x_{2} are independent, the unconditional likelihood is
where 0 ≤ x_{1} ≤ n, 0 ≤ x_{2} ≤ N − n, and 0 ≤ π_{ i } ≤ 1, i=1,2.
If we define
and
then θ is the parameter of interest, representing the log odds ratio of the GO term, and λ is a nuisance parameter. Under the new parametrization, the unconditional likelihood function (2) is
where 0 ≤ x_{1} ≤ n and 0 ≤ x_{2} ≤ N − n.
In equation (5), we take the interest parameter θ and also the nuisance parameter λ into consideration. Consider statistics T and S, functions of x_{1} and x_{2}, such that T(X_{1},X_{2}) = X_{1} and S(X_{1},X_{2}) = X_{1}+X_{2}. Thus, T represents the number of DE genes in a GO category, and S represents the number of total genes in a GO category. Let t and s be the observed values of T and S. The probability mass function of T(x_{1},x_{2}) = t evaluated at S(x_{1},x_{2}) = x_{1}+x_{2} = s, say Pr(T = tS = s;θ,λ,N,n), does not depend on the nuisance parameter λ[19]. See also Example 8.47 of Severini[20]. Thus, we derive the conditional probability mass function
understood as a function of t.
By eliminating the nuisance parameter λ, we can reduce the original data x_{1} and x_{2} by considering the statistic T = t. However, the use of the conditional probability mass function requires some justification because of concerns about losing information during the conditioning process. Unfortunately, in the presence of the nuisance parameter, the statistic S(X_{1},X_{2}) = X_{1}+X_{2} is not an ancillary statistic for the parameter of interest. In other words, the probability mass function of the conditional variable S(X_{1},X_{2}) may contain some information about parameter θ[20]. However, following the explanation of BarndorNielsen and Cox ([21], $\S 2.5$), the expectation value of statistic S(X_{1},X_{2}) equals the nuisance parameter. Hence, from the observation of S(X_{1},X_{2}) alone, the distribution of S(X_{1},X_{2}) contains little information about θ[21]. S(X_{1},X_{2}) satisfies the other 3 conditions of an ancillary statistic defined by BarndorNielsen and Cox [21]: parameters θ and λ are variation independent; (T(X_{1},X_{2}),S(X_{1},X_{2})) is the minimal sufficient statistic; and the distribution of T(X_{1},X_{2}), given S(X_{1},X_{2}) = s, is independent of the parameter of interest, θ, given the nuisance parameter λ. Therefore, the probability mass function of S(X_{1},X_{2}) contains little information about the value of θ.
Hypotheses and LFDRs
Considering GO term i, we denote the T, S, t, s, and θ used in equation (6) as T_{ i }, S_{ i }, T_{ i }, S_{ i }, and θ_{ i }. From Table 1, hypothesis comparison (1) of GO term i is equivalent to
Let S = 〈S_{1},S_{2},⋯,S_{ m }〉 and S = 〈S_{1},S_{2},⋯,S_{ m }〉. Let BF_{ i } denote the Bayes factor of GO term i:
It is called the Bayes factor because it yields posterior odds when multiplied by prior odds. More precisely, the posterior odds of the alternative hypothesis corresponding to GO term i is
where Π_{0} is the prior conditional probability that a GO term is not enriched for the preselected genes given s, i.e., Π_{0} = Pr(θ_{ i } = 0S = s). Thus, (1 − Π_{0})/Π_{0} is the prior odds of the alternative hypothesis of enrichment. According to Bayes’ theorem, the LFDR of GO term i is
where ω_{ i } is defined in equation (9).
Methods
This section is divided into two parts:

1.
Previous LFDR estimators. While not unique to this paper, these methods are included for comparison.

2.
New LFDR estimators. Our main methodological innovations are the uses of a conditional probability mass function and of normalized maximum likelihood for LFDR estimation.
The other original contributions of this paper are the estimator comparisons of the next section. The comparisons are made by simulation and by a case study.
Previous LFDR estimators
Binomialbased LFDR estimator
The version of the FDR that generalizes the LFDR is the nonlocal FDR, which is defined as the ratio of the expected number of false discoveries to the expected total number of discoveries [17]. In our running example, a discovery of enrichment is a rejection of the null hypothesis of nonenrichment at some significance level α, and a false discovery of enrichment is a discovery of enrichment corresponding to a case of no actual enrichment. (This FDR has been called the “Bayesian FDR” [22] to distinguish it from the FDR of Benjamini and Hochberg [9]).
Let α denote any significance level chosen to be between 0 and 1. For all GO terms of interest, the nonlocal FDR may be estimated by
where m is the number of GO terms, p_{ j } is the pvalue of GO term j, and ${\mathbf{1}}_{\{{p}_{j}\le \alpha \}}$ is the indicator such that ${\mathbf{1}}_{\{{p}_{j}\le \alpha \}}=1$ if p_{ j } ≤ α is true and ${\mathbf{1}}_{\{{p}_{j}\le \alpha \}}=0$ otherwise. Thus, $\sum _{j=1}^{m}{\mathbf{1}}_{\{{p}_{j}\le \alpha \}}$ represents the number of discoveries of enriched GO terms, and m α estimates the number of such discoveries that are false.
Let r_{ i } be the rank of the pvalue of GO term i, e.g., r_{ i } = 1 if the pvalue of GO term i is the smallest among all pvalues of m GO terms. Based on a modification of equation (11), the binomialbased estimator(BBE) of LFDR of the GO term i is
It is conservative in the sense that it tends to overestimate LFDR [17].
Histogrambased LFDR estimator
Efron [15, 16] devised reliable histogrambased LFDR estimators for a range of applications in microarray gene expression analysis and other problems of largescale inference. Let z_{ i } = Φ^{− 1}(p_{ i }) be the ztransformed statistic of GO term i, where Φ is the standard normal cumulative distribution function (cdf) and p_{ i } is the 2sided pvalue of GO term i. For each GO term, the density is a mixture of the form
where f_{0} is the density function of z for the nonenriched GO terms, f_{1} is that for the enriched GO terms, and Π_{0} is the probability that a GO term is nonenriched. The histogrambased LFDR of GO term i is estimated by equation (14):
where $\widehat{f}$ is the estimator of f that is estimated by a nonparametric Poisson regression method [15, 16]. We call ${\hat{\text{LFDR}}}_{i}$the histogrambased estimator (HBE) if the density function f_{0} is assumed to be standard normal, N(0,1), and the histogrambased estimator with empirical null (HBEEN) if the density function f_{0} is estimated based on the truncated maximum likelihood technique of [16]. Dalmasso et al. [23] compared the precursor of HBEEN [15] to other LFDR estimators.
New LFDR estimators
Type II maximum likelihood estimator
Bickel [17] follows Good [24] in calling the maximization of likelihood over a hyperparameter Type II maximum likelihood to distinguish it from the usual Type I maximum likelihood, which pertains only to models that lack random parameters. Type II maximum likelihood has been applied to parametric mixture models (PMMs) for the analysis of microarray data [25, 26], proteomics data [18], and genetic association data [27]. In this section, we adapt the approach to the feature enrichment problem by using the conditional probability mass function defined above. The particular models we use in this framework correspond to new methods of enrichment analysis.
Let $\mathcal{G}\left(\mathbf{s}\right)=\left\{{g}_{\theta}\right(\bullet \left\mathbf{s}\right);\theta \ge 0\}$ be a parametric family of probability mass functions with
where f_{ θ }(∙s) is defined in equation (6). We define the kcomponent PMM as
where θ_{0} = 0 and θ_{ j }≠θ_{ J } for all $j,J\in \left\{0,\dots ,k1\right\}$ such that j≠J.
Let T = 〈T_{1},T_{2},⋯,T_{ m }〉 and T = 〈T_{1},T_{2},⋯,T_{ m }〉 be vectors of the T_{ i }s and T_{ i }s used in equation (8). Assuming T_{ i } is independent of T_{ j } and S_{ j } for all $i,j\in \left\{1,\dots ,m\right\}$ such that j≠J. i≠j, the joint probability mass function is
where S_{ i } is the observed value of S_{ i } for GO term i, and S=〈S_{1},S_{2},⋯,S_{ m }〉.
Moreover, we assume that for a given number of genes in GO term i, T_{ i }(i = 1,…,m) satisfies the kcomponent PMM shown in equation (16). In other words, we assume that the possible log odds ratios of GO term i are the θ_{0},θ_{1},θ_{2},…,θ_{k − 1} of equation (16) if the alternative hypothesis H_{1} in hypothesis comparison (7) is true.
Therefore, the loglikelihood function under the kcomponent PMM for all GO terms is
The LFDR of GO term i is estimated by
where ${\hat{\theta}}_{1},\dots ,{\hat{\theta}}_{k1}$ and ${\hat{\Pi}}_{0},\dots ,{\hat{\Pi}}_{k1}$ are maximum likelihood estimates of θ_{1},…,θ_{k − 1} and Π_{0},…,Π_{k − 1} in equation (18). We call ${\hat{\text{LFDR}}}_{i}^{\left(k\right)}$ the k component maximum likelihood estimator (MLEk). Our LFDRenrich and LFDRhat software suites of R functions that implement MLE2 and MLE3 are now available at http://www.statomics.com. Moreover, ${\hat{\theta}}_{i}$ (i = 1,…,k − 1;k = 2,3) and ${\hat{\Pi}}_{j}$ (j = 0,…,k − 1;k = 2,3), also in LFDRenrich and LFDRhat, are maximum likelihood estimators of θ_{ i } (i = 1,…,k − 1;k = 2,3) and Π_{ j } (j = 0,…,k − 1;k = 2,3) under given constraints.
LFDR estimator based on the normalized maximum likelihood
Combining equations (9)(10), we obtain
Therefore, given a guessed value of Π_{0}, we may use an estimator of the Bayes factor to estimate the LFDR of a GO term.
We now develop such an estimator of the Bayes factor. For GO category i, let ${\mathcal{E}}_{i}$ stand for the set of all probability mass functions defined on {0,1,…,s_{ i }}, the set of all possible values of T_{ i }. Based on hypothesis comparison (7), the set of log odds ratios, denoted as $\mathcal{\mathscr{H}}$, is {0} under the null hypothesis and is $\mathbb{R}\setminus \left\{0\right\}=\left\{\theta \in \mathbb{R}:\theta \ne 0\right\}$, the set of all real values except 0, under the alternative hypothesis. With the assumption that random variable T_{ i } is independent of random variable S_{ j } for any i≠j, the regret of a predictive mass function $\stackrel{\u0304}{f}\in {\mathcal{E}}_{i}$ is a measure of how well it predicts the observed value t_{ i }∈{0,1,…,s_{ i }}. The regret is defined as
where ${\widehat{\theta}}_{i}\left({t}_{i}\right{s}_{i})$ is a Type I MLE with respect to $\mathcal{\mathscr{H}}$ under observed values T_{ i } given S_{ i }[28, 29].
For all members of ${\mathcal{E}}_{i}$, the optimal predictive conditional probability mass function of GO category i and the hypothesis that ${\theta}_{i}\in \mathcal{\mathscr{H}}$ is denoted by ${f}_{i}^{\u2020}(\bullet {s}_{i};\mathcal{\mathscr{H}})$. It minimizes the maximal regret in sample space {0,1,…,s_{ i }} in the sense that it satisfies
It is well known [28] that the predictive probability mass function that satisfies equation (22) is
where f_{ θ }(t_{ i }s_{ i }) is the conditional probability mass function defined in equation (6), and ${\mathcal{K}}_{i}^{\u2020}\left(\mathcal{\mathscr{H}}\right)$ is a constant defined as
where
We call ${f}_{i}^{\u2020}\left({t}_{i}\right{s}_{i};\mathcal{\mathscr{H}})$ the normalized maximum likelihood (NML) associated with the hypothesis that ${\theta}_{i}\in \mathcal{\mathscr{H}}$.
Thus, BF_{ i } is approximated by
which we call the NML ratio. (More generally, the logarithm of an NML ratio is interpreted as a measure of the evidential support for the alternative hypothesis over the null hypothesis [29, 30]). Therefore, by combining equations (8) and (9), if we guess the prior probability Π_{0}, the LFDR estimate of GO category i in the hypothesis comparison (7) is
where ${\hat{\text{BF}}}_{i}^{\u2020}$ is defined in equation (26). We call this LFDR estimator the NML estimator (NMLE).
To assess the reliability of NML ratio ${\hat{\text{BF}}}_{i}^{\u2020}$ for a particular data set, it will be compared to an empirical Bayes estimate of the Bayes factor that unlike NML, simultaneously takes all GO terms into account. Equations (19) and (20) suggest
as the empirical Bayes estimator of BF_{ i }.
Results and discussion
In this section, we compared the LFDR estimators using simulation data and breast cancer data.
For each GO category, the pvalue used in BBE to estimate LFDR is computed based on the 2sided Fisher’s exact test. In computing MLEk (k = 2,3), θ_{ i } (i = 1,…,k − 1) in equation (18) was constrained to lie between 0 and 10, whereas Π_{ i } (i = 0,…,k − 1) in equation (18) was allowed to take any value between 0 and 1 such that $\sum _{i=0}^{k1}{\Pi}_{i}=1$.
Simulation studies
The aim of the following simulation studies is to compare the LFDR estimation biases of BBE, MLE2, MLE3, HBE, and HBEEN. NMLE is not taken into account because its performance depends not only on the data, but also on the specified prior probability Π_{0}.
The simulation setting involves 10,000 genes in a microarray with 200 genes identified as DE and 100 GO terms. We conducted a separate simulation study using each of these values of Π_{0}: 50%, 60%, 70%, 80%, 90%, and 94%.
Since the PMM behind MLE is optimal when the number of enriched GO terms is equal to the number nonenriched GO terms, we assessed the sensitivity of MLE to that symmetry assumption by using strongly asymmetric log odds ratios and by using symmetric ones. For each GO term, two configurations were used in this simulation to choose log odds ratios: the asymmetric configuration shown in equation (29) and the symmetric configuration shown in equation (30). We used these values of odds ratio of the i th GO term:
Considering the log odds ratios of all GO terms in each simulation study, we generated Table 1 for GO term i and for each of the 20 simulated data sets as follows:

x_{1} is generated from a binomial distribution with parameter π_{1} used in equation (2); π_{1} is a real value randomly picked from 0 to 1.

x_{2} is obtained from a binomial distribution with parameter ${\pi}_{2}={\left[\frac{(1{\pi}_{1})\times {2}^{{\varphi}_{i}}}{{\pi}_{1}}+1\right]}^{1}$, obtained by solving
$${\varphi}_{i}=\underset{2}{log}[{\pi}_{1}/(1{\pi}_{1}\left)\right]\underset{2}{log}[{\pi}_{2}/(1{\pi}_{2}\left)\right].$$(31)
Thus, according to equation (4), we obtain ϕ_{ i } = θ_{ i } log2e for GO term i. Each of those artificial data sets represents what might have been a real data set such as that of the next subsection.
The pvalue of each GO term used in BBE, HBE, and HBEEN is obtained from the 2sided Fisher’s exact test. The kcomponent PMM (k = 2 or k = 3) used in MLE is shown in equation (16) with ${\Pi}_{j}=\left(1{\Pi}_{0}\right)/k\phantom{\rule{0.3em}{0ex}}\left[j=1,\dots ,k\right]$ and ${g}_{{\theta}_{i}}\left({t}_{i}\right{s}_{i})$ defined in equation (15). For every log odds ratio sequence, we estimated the LFDR for each GO term and each data set using BBE, MLE2, MLE3, HBE, and HBEEN. We compared the performances of the 5 estimators by means of computing the absolute value of the estimated LFDR bias. The true LFDR is computed by equation (10), where
and f_{1}(t_{ i }) is computed by
where f_{ θ }(ts) is defined in equation (6).
Figure 1 shows the performance comparisons of the 5 LFDR estimators (i.e., BBE, MLE2, MLE3, HBE, and HBEEN) for simulation data obtained from asymmetric and symmetric log odds ratios. The absolute LFDR biases estimated by BBE, MLE2, MLE3, and HBEEN are similar. The absolute bias of LFDR estimated by HBE on symmetric log odds ratios is a little higher than that on asymmetric log odds ratios when the proportion of nonenriched GO terms is greater than 80%. Therefore, the estimated LFDR biases of the estimators are not strongly affected by whether the log odds ratios are symmetric or asymmetric.
To assess the performance of the 5 estimators for smaller GO terms, we added simulation studies using 2, 4, 8, and 32 as the total number of GO terms. The proportion of nonenriched GO terms (Π_{0}) and log odds ratios of simulation studies are shown in Table 2. The simulation studies were otherwise the same as those for 100 GO terms. Figure 2 shows the performance of LFDR estimators by means of computing the absolute estimated LFDR bias for 2, 4, 8, and 32 GO terms with log odds ratios based on formulas shown in Table 2.
Considering every (m,Π_{0}) pair of the simulation studies with symmetric log odds ratios for the case of 100 GO terms, we recorded the LFDR estimator with the lowest absolute estimated LFDR bias among the 5 LFDR estimators (BBE, MLE2, MLE3, HBE, and HBEEN). Moreover, we determined the maximum absolute LFDR bias over the proportion of nonenriched GO terms (Π_{0}) in order to evaluate the worstcase bias of each estimator at each value of m. Figure 3 shows the results.
Breast cancer data analysis
The singlechannel microarray data set used here to illustrate our new methods is from an experiment applying an estrogen treatment to cells of a human breast cancer cell line [31]. The Affymetrix human genome U95Av2 genechip data are from four samples from an estrogen receptor positive breast cancer cell line. Two of the samples were exposed to estrogen and then harvested after 10 hours. The remaining two samples were left untreated and then harvested after 10 hours. For simplicity of terminology, we call probes in the microarray experiment “genes.” The relevant data consist of measurements of gene expression across the reference class of 12,625 genes. The purpose of the study was to determine which genes are affected by the estrogen treatment. (For further information concerning the data, see Gentleman et al. [32].)
We applied the R function expresso in the affy package [33] of Bioconductor [34] to convert the raw probe intensities from the the CEL data files to logarithms of gene expression levels without background correction. In doing so, we applied the “quantiles,” “pmonly,” and “medianpolish” [35] preprocessing settings.
We selected as genes of interest those that were differentially expressed between the treatment group and the control group according to the following criterion. Using the LFDR as the probability that a gene is EE, we considered genes with LFDR estimates below 0.2 as DE. In other words, we selected as DE genes those that were differentially expressed with estimated posterior probability of at least 80%. Considering four samples of each gene in the microarray, we used the unpaired ttest with equal variances to compute the pvalue. The LFDR of every gene is estimated using the theoretical null hypothesis method of Efron [15, 16]; the empirical null hypotheses method can lead to excessive bias due to deviations from normality [36]. When we compared gene expression data for the presence and absence of estrogen after 10 hours of exposure, we obtained 74 DE genes.
Defining unrelated pairs of GO terms as those that do not share any common ancestor, we selected for analysis all unrelated GO molecular function terms with at least 1 DE gene, thereby obtaining a total of 82 GO terms of interest. Figure 4 compares the BBE to the MLEs based on the 2component (MLE2) and 3component (MLE3) PMM. Figure 5 displays the probability mass of GO:0005524 under the null and alternative hypotheses of hypothesis comparison (7). Figure 6 compares MLEbased estimates of the Bayes factor given by equation (28) to the NML ratios given by equation (28).
For two GO terms, opposite conclusions would be drawn about their enrichment, depending on which estimator is used. As seen in Figure 4, the estimated LFDRs of GO:0051082 and GO:0005524 using MLE2 were 100%. However, the LFDRs estimated by MLE3 were essentially 0.
Using the MLE formula shown in equation (19), and the kcomponent PMM shown in equation (16), we conclude that the sensitivity of the LFDRs of GO term i estimated by MLE2 and MLE3 depended mainly on the sensitivity of the Bayes factor, based on the number of PMM components. Comparing the probability masses of GO:0005524, based on the 2 and 3component PMMs shown in Figure 5, we found that the probability mass of GO:0005524 under the null hypothesis is larger than that under the alternative hypothesis based on the 2component PMM (left plot in Figure 5). In contrast, the probability mass under the null hypothesis is smaller than that under the alternative hypothesis based on the 3component PMM (right plot in Figure 5). Thus, the LFDR estimated by MLE is strongly dependent on the number of PMM components.
While a real data set can in that way indicate the impact of selecting an appropriate method, that impact does not in itself say which method has lowest bias. For that, we rely on the simulation study of the previous subsection.
Conclusions
As seen in Figure 1 and Figure 2, HBE and HBEEN have relative high biases for a small number and a medium number of GO terms, respectively. The performance comparison displayed in the lefthand side of Figure 3 indicates that BBE contains the lowest minimum estimated LFDR bias for a small number of GO terms (i.e., 232 GO terms) when the proportion of nonenriched GO terms is 1. Although the minimum bias of BEE is not the lowest for some Π_{0}s under a small number of GO terms, it is very close to the lowest value of bias based on plots shown in Figure 2. The righthand side of Figure 3 indicates that MLE3 has the lowest maximum absolute estimated LFDR bias in 100 GO terms. MLE exhibits bias similar to that of BBE when the number of GO terms is much larger than k except for when the proportion of nonenriched GO terms is high (close to 1). Moreover, MLE3 has lower bias than MLE2 as an LFDR estimator. Due to its conservatism and freedom from PMM, we recommend using BBE for a small number of GO terms of interest (232 GO terms) and MLE for a medium number of GO terms of interest (100 GO terms).
Finally, we recommend that NMLE be used when there is only 1 GO term of interest since none of the other estimators is able to estimate LFDR in such a case except by conservatively giving 1 as the estimate. Otherwise, unless Π_{0} is known with sufficient accuracy, NMLE is not recommended since it depends not only on the data but also on a guess of the value of Π_{0}, which in the absence of strong prior information, is often set to the default value of 50%. A closely related approach is to use the logarithm of the NML ratio as a measure of statistical support for the enrichment hypothesis [30] without guessing Π_{0}. By using 10 and 100 as thresholds of the approximate Bayes factors from equations (26) and (28) to determine whether a GO term is enriched, we reached similar conclusions with both NML and MLE (Figure 6). Thus, in our data set, the NML ratio tends to estimate the Bayes factor almost as accurately as methods that simultaneously use information across GO terms. While we do not expect the same for all data sets, we note that similar results have been found for an application of a modified NML [29] to a proteomics data set [30].
References
 1.
Altshuler D, Daly MJ, Lander ES: Genetic mapping in human disease. Science 2008, 322: 881888. 10.1126/science.1156409
 2.
Rhee SY, Wood V, Dolinski K, Draghici S: Use and misuse of the gene ontology annotations. Nat Rev Genet 2008,9(7):509515. 10.1038/nrg2363
 3.
Kanehisa M, Goto S: KEGG: Kyoto Encyclopedia of Genes and Genome. Nucleic Acids Res 2000, 28: 2730. 10.1093/nar/28.1.27
 4.
Dennis G, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, Lempicki R: DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol 2003, 4: P3. 10.1186/gb200345p3
 5.
Doniger SW, Salomonis N, Dahlquist KD, Vranizan K, Lawlor SC, Conklin BR: MAPPFinder: using gene ontology and GenMAPP to create a global geneexpression profile from microarray data. Genome Biol 2003, 4: R7. 10.1186/gb200341r7
 6.
Khatri P, Draghici S, Ostermeier G, Krawetz S: Profiling gene expression using ontoexpress. Genomics 2002, 79: 266270. 10.1006/geno.2002.6698
 7.
Zeeberg BR, Feng W, Wang G, Wang MD, Fojo AT, Sunshine M, Narasimhan S, Kane DW, Reinhold WC, Lababidi S, Bussey KJ, Riss J, Barrett JC, Weinstein JN: GoMiner: a resource for biological interpretation of genomic and proteomic data. Genome Biol 2003, 4: R28. 10.1186/gb200344r28
 8.
Huang DW, Sherman BT, Lempicki RA: Bioinformatics enrichment tools paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res 2009, 37: 113. 10.1093/nar/gkn923
 9.
Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B 1995, 57: 289300.
 10.
Min JL, Barrett A, Watts T, Pettersson FH, Lockstone HE, Lindgren CM, Taylor JM, Allen M, Zondervan KT, McCarthy MI: Variability of gene expression profiles in human blood and lymphoblastoid cell lines. BMC Genomics 2010, 11: 96. 10.1186/147121641196
 11.
Reyal F, van Vliet MH, Armstrong NJ, Horlings HM, de Visser KE, Kok M, Teschendorff AE, Mook S, van’t Veer L, Caldas C, Salmon RJ, Vijver MJVD, Wessels LFA: A comprehensive analysis of prognostic signatures reveals the high predictive capacity of the proliferation, immune response and RNA splicing modules in breast cancer. Breast Cancer Res 2008, 10: R93. 10.1186/bcr2192
 12.
Wang R, Bencic D, Lazorchak J, Villeneuve D, Ankley GT: Transcriptional regulatory dynamics of the hypothalamicpituitarygonadal axis and its peripheral pathways as impacted by the 3beta HSD inhibitor trilostane in zebrafish (Danio rerio). Ecotoxicol Environ Saf 2011, 74: 14611470. 10.1016/j.ecoenv.2011.05.001
 13.
Storey JD: The positive false discovery rate: a Bayesian interpretation and the qvalue. Ann Stat 2003, 31: 20132035. 10.1214/aos/1074290335
 14.
Hong WJ, Tibshirani R, Chu G: Local false discovery rate facilitates comparison of different microarray experiments. Nucleic Acids Res 2009, 37: 74837497. 10.1093/nar/gkp813
 15.
Efron B: Largescale simultaneous hypothesis testing: The choice of a null hypothesis. J Am Stat Assoc 2004, 99: 96104. 10.1198/016214504000000089
 16.
Efron B: LargeScale Inference: Empirical Bayes Methods for Estimation,Testing, and Prediction Cambridge. Cambridge University Press; 2010.
 17.
Bickel DR: Simple estimators of false discovery rates given as few as one or two pvalues without strong parametric assumptions. Stat Appl Genet Mol Biol in press in press
 18.
Bickel DR: Smallscale inference: empirical Bayes and confidence methods for as few as a single comparison. Tech Rep, Ottawa Inst Syst Biol; 2011:arXiv:1104.0341arXiv:1104.0341.
 19.
Padilla M, Bickel DR: Empirical Bayes methods corrected for small numbers of tests. Stat Appl Genet Mol Biol 2012,11(5):art. 4.
 20.
Severini T: Likelihood Methods in Statistics Oxford. Oxford University Press; 2000.
 21.
BarndorffNielsen OE, Cox DR: Inference and Asymptotics. London: CRC Press; 1994.
 22.
Efron B, Tibshirani R: Empirical Bayes methods and false discovery rates for microarrays. Genet Epidemiol 2002, 23: 7086. 10.1002/gepi.1124
 23.
Dalmasso C, BarHen A, Broët P: A constrained polynomial regression procedure for estimating the local false discovery rate. BMC Bioinformatics 2007, 8: 229. 10.1186/147121058229
 24.
Good IJ: How to estimate probabilities. IMA J Appl Math 1966, 2: 364383. 10.1093/imamat/2.4.364
 25.
Pawitan Y, Murthy K, Michiels S, Ploner A: Bias in the estimation of false discovery rate in microarray studies. Bioinformatics 2005, 21: 38653872. 10.1093/bioinformatics/bti626
 26.
Muralidharan O: An empirical Bayes mixture method for effect size and false discovery rate estimation. Ann Appl Stat 2010, 4: 422438.
 27.
Yang Y, Aghababazadeh FA, Bickel DR: Parametric estimation of the local false discovery rate for identifying genetic associations. IEEE/ACM Trans Comput Biol Bioinformatics 2012. online ahead of print at http://dx.doi.org/10.1109/TCBB.2012.140 online ahead of print at
 28.
Grünwald PD: The Minimum Description Length Principle. London: MIT Press; 2007.
 29.
Bickel DR: A predictive approach to measuring the strength of statistical evidence for single and multiple comparisons. Can J Stat 2011, 39: 610631. 10.1002/cjs.10109
 30.
Bickel DR: Minimaxoptimal strength of statistical evidence for a composite alternative hypothesis. Int Stat Rev 2013. in press. 2011 version available at arXiv:1101.0305 in press. 2011 version available at arXiv:1101.0305
 31.
Scholtens D, Miron A, Merchant FM, Miller A, Miron PL, Iglehart JD, Gentleman R: Analyzing factorial designed microarray experiments. J Multivariate Anal 2004, 90: 1943. 10.1016/j.jmva.2004.02.004
 32.
Gentleman R, Carey V, Huber W, Irizarry R, Dudoit S (Eds): Bioinformatics and Computational Biology Solutions Using R and Bioconductor. New York: Springer; 2005.
 33.
Gautier L, Cope L, Bolstad BM, Irizarry RA: Affy—analysis of Affymetrix Gene Chip data at the probe level. Bioinformatics 2004,20(3):307315. 10.1093/bioinformatics/btg405
 34.
Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, Hornik K, Hothorn T, Huber W, Iacus S, Irizarry R, Leisch F, Li C, Maechler M, Rossini AJ, Sawitzki G, Smith C, Smyth G, Tierney L, Yang JYH, Zhang J: Bioconductor: open software development for computational biology and bioinformatics. Genome Biol 2004, 5: R80. 10.1186/gb2004510r80
 35.
Tukey JW: Exploratory Data Analysis. Reading: AddisonWesley; 1977.
 36.
Bickel DR: Estimating the null distribution to adjust observed confidence levels for genomescale screening. Biometrics 2011, 67: 363370. 10.1111/j.15410420.2010.01491.x
 37.
Jeffreys H: Theory of Probability. London: Oxford University Press; 1948.
 38.
Bickel DR: The strength of statistical evidence for composite hypotheses: inference to the best explanation. Statistica Sinica 2012, 22: 11471198.
Acknowledgements
We thank the two anonymous reviewers for comments that led to improvements in the manuscript, particularly those that led to clearer communication or to additional simulation studies. We also thank both Editage and Donna Reeder for the detailed copyediting. We are grateful to Corey Yanofsky and Ye Yang for the useful discussions. This work was partially supported by the Natural Sciences and Engineering Research Council of Canada, by the Canada Foundation for Innovation, by the Ministry of Research and Innovation of Ontario, and by the Faculty of Medicine of the University of Ottawa.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
ZY implemented the NMLE function, codesigned and executed the simulation study, carried out the comparison among 3 LFDR estimators (BBE, MLE, and NMLE) using the breast cancer data, and drafted the manuscript. ZL implemented the 3component MLE functions. DRB suggested and guided the project, contributed to writing the paper, codesigned the simulation study, and provided the BBE function. All authors have read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Yang, Z., Li, Z. & Bickel, D.R. Empirical Bayes estimation of posterior probabilities of enrichment: A comparative study of five estimators of the local false discovery rate. BMC Bioinformatics 14, 87 (2013). https://doi.org/10.1186/147121051487
Received:
Accepted:
Published:
Keywords
 Gene Ontology
 False Discovery Rate
 Maximum Likelihood Estimator
 Differentially Express
 Differentially Express Gene