BACKGROUND
with the development of microarray technologies, researchers now can measure the relative expressions of tens of thousands of genes simultaneously.
however, the number of replicates per gene is usually small, far less than the number of genes.
many statistical methods have been developed to identify differentially expressed  genes.
the use of fold change is among the first practice.
it can be inefficient and erroneous because of the additional uncertainty induced by dividing two intensity values.
there are variants of student's t test procedure that conduct a test on each individual gene and then correct for multiple comparisons.
the problem is, with a large number of tests and a small number of replicates per gene, the statistics are very unstable.
for example, a large t statistic might arise because of an extremely small variance, even with a minor difference in the expression.
the disadvantage of fold-change approach and t test has been pointed out by a number of authors  <cit> .
there are approaches proposed to improve estimation of gene variances by borrowing strength across genes  <cit> .
despite the flaw, fold change and t test are the most intuitive approaches and they have been applied widely in practice.
to control the error rate, many researchers use fold change and t test together, hoping that the double filtering will provide extra confidence in the test results.
specifically, a gene is flagged as de only if the p-value from t test is smaller than a certain threshold and the fold change is greater than a cutoff value.
for example, in  <cit> ,  <dig> genes were found to be de with two cutoff values .
there are numerous examples in the literature that implement the double filtering procedure with fold change and t statistic  <cit> .
we argue, however, that the double filtering procedure provides higher confidence mainly because it produces a shorter list of selected genes.
given the same number of genes selected, a well constructed shrinkage test can significantly outperform the double filtering method.
fold change takes the ratio of a gene's average expression levels under two conditions.
it is usually calculated as the difference on the log <dig> scale.
let xij be the log-transformed expression measurement of the ith gene on the jth array under the control , and yik be the log-transformed expression measurement of the ith gene on the kth array under the treatment .
we define  and .
fold change is computed by
as for the traditional t test, it is usually calculated on the log <dig> scale to adjust for the skewness in the original gene expression measurements.
the t statistic is then computed by
where  is the pooled variance of xij and yik.
comparing  and , it is obvious that fold change and t statistic are based on two contradicting assumptions.
the underlying assumption of fold change is that all genes share a common variance , which is implied by the omission of the variance component in .
on the other hand, the inclusion of  in  suggests that t test assumes gene-specific variances.
in order for a gene to be flagged as de, the double filtering procedure would require the gene to have extreme test scores under the common variance assumption as well as under the gene-specific variance assumption.
it is analogous to using the intersection of the rejection regions defined by fold change and t statistic.
assuming a common variance for all the genes apparently is an oversimplification.
the assumption of gene-specific variances, however, leads to unstable estimates due to limited replicates from each gene.
a more realistic assumption might lie in between the two extremes, i.e., modeling gene variances by a mixture of two components, one being a point mass at the common variance, another being a continuous distribution for the gene-specific variances.
under this mixture variance assumption, a de gene could have a large fold change or a large t statistic, but not necessarily both.
taking intersection of the rejection regions flagged by fold change and t statistic, as is adopted by the double filtering procedure, might not be the best strategy under the mixture variance assumption.
the goal of the paper is not to propose a new testing procedure in microarray analysis, but to provide insight on the drawback of the widely used double filtering procedure with fold change and t test.
we present a theoretically most powerful likelihood ratio  test under the mixture variance assumption.
we further demonstrate that two shrinkage test statistics, one from the bayesian model  <cit>  and the other from the significance analysis of microarrays  test  <cit> , can be united as approximations to the lr test.
this association explains why those shrinkage methods can considerably outperform the double filtering procedure.
a simulation study and real microarray data analyses are then presented to compare the shrinkage tests and the double filtering procedure.
methods
a likelihood ratio test
for gene i, we use fi = pvfi <dig> + fi <dig>  a mixture of two components fi <dig> and fi <dig>  to denote the density under the null hypothesis that the gene is not de under two experiment conditions.
density fi <dig> is defined under the gene-specific variance assumption, fi <dig> is defined under the common variance assumption, and pv is the mixing probability.
similarly, we use gi = pvgi <dig> + gi <dig> to denote the density under the alternative hypothesis, with gi <dig> and gi <dig> defined in a similar fashion as fi <dig> and fi <dig>  for example, in the context of testing de genes, we can assume fi <dig> = n , fi <dig> = n, gi <dig> = n, and gi <dig> = n, where  is the assumed common variance,  is the gene-specific variance, μi is the mean expression level under the control, and Δi is the difference in the expression levels between two conditions.
under the null hypothesis h <dig> : Δi =  <dig>  the likelihood ratio test statistic, which is the most powerful among all test statistics, is
the ri statistic is a weighted sum of two ratios gi1/fi <dig> and gi2/fi <dig>  with weight wi = pvfi1/.
under the normality assumption, it is easy to show that ri = wih <dig> + h <dig>  where fci and ti are fold change and t statistic, as defined in  and .
both h <dig> and h <dig> are monotonic increasing functions.
the rejection region of the lr test is defined by ri >λr, where λr is the threshold to attain a certain test size.
in order to reject h <dig>  it requires that either |fci| is large, or |ti| is large, or both.
in this sense, the lr test rejection region is more like a union of the rejection regions defined by fold change and t statistic.
on the other hand, the double filtering procedure with fold change and t statistic would require both |fci| and |ti| to be large.
this practice is analogous to using the intersection of the two rejections regions determined by |fci| and |ti|.
compared with the lr test, the double filtering procedure will lose power.
the "loss of power" has two meanings.
first, for a given false discovery rate , the double filtering procedure produces a shorter list of identified genes for further investigation.
second, for a given number of identified genes, the list produced by the double filtering procedure has a higher fdr.
the double filtering procedure offers a false sense of confidence by producing a shorter list.
the lr test statistic ri requires one to know the true values of parameters p, μi, , , and Δi, which are usually unknown in reality.
one strategy is to estimate ri by , where the maximum likelihood estimates  of the unknown parameters are plugged into .
unfortunately, with a small number of replicates from each gene, the mle would be extremely unstable and lead to unsatisfactory testing results.
a bayesian model  <cit>  was constructed under the mixture variance assumption to detect de genes.
the inference is made based on the marginal posterior probability of a gene being de, denoted by zi = p. here x = {xij} and y = {yik} are the collection of gene expression data under the two conditions.
we will show that, like , zi is also an approximation to ri.
the difference between  and zi is that the former plugs in the point estimates  of unknown parameters, while the latter marginalizes the unknown parameters with respect to their posterior distribution.
in the bayesian inference, the uncertainty from various sources are accounted for in a probabilistic fashion.
similar to the bayesian mixture model, some existing methods also try to strike a balance between the two extreme assumptions of a common variance and gene-specific variances.
the sam statistic slightly modifies the t-statistic by adding a constant to the estimated gene-specific standard deviation in the denominator.
we will present it as being motivated by a mixture model on the variances .
furthermore, the sam statistic can be directly written as a weighted sum of t statistic and fold change.
thus both the bayesian method and the sam method are approximations to the lr test under the mixture variance assumption, and they can achieve better performance than the double filtering procedure.
the bayesian mixture model
cao et al.
<cit>  proposed a bayesian mixture model to identify de genes, which has shown comparable performance to frequentist shrinkage methods  <cit> .
with parameters  defined similarly as in the lr test, gene expression measurements xij and yij are modeled by normal distributions with a mixture structure on the variances,
a latent variable ri is used to model the expression status of the ith gene,,
where ri = 0/ <dig> indicates that gene i is non-de/de and it is modeled by a bernoulli distribution: ri | pr ~ bernoulli.
for  and , it is assumed that  and  ~ ig.
here ig denotes an inverse gamma distribution with mean b/.
the other hyper-priors include, μi ~ n, pr ~ u, and pv ~ u. more details can be found in  <cit> .
we make inference based on zi = p = p, the marginal posterior probability that gene i is de.
a gene is flagged as de if zi >λz, where λz is a certain cutoff.
we argue that the bayesian rejection region defined by zi >λz is an approximation to the lr test rejection region defined by ri >λr.
first we have
here p is the joint posterior distribution of , marginalized with respect to other random parameters .
it is easy to show that
given parameters , p is an increasing function of ri.
rejecting h <dig> for ri >λr is equivalent to rejecting for p >λz, with λz = λr/.
thus the two test statistics, p and ri, are equivalent.
expression  demonstrates that zi is obtained from p by integrating with respect to the unknown parameters under the joint posterior distribution.
if the integral does not have a closed form, we can conduct numerical integration to calculate zi through the gibbs sampling algorithm  <cit> .
the uncertainty from those unknown parameters are accounted for in a probabilistic fashion.
it is in this sense that we consider zi a good approximation to the lr test statistic ri.
the sam test
the sam statistic  <cit>  is defined as
where si is the gene-specific standard deviation, and s <dig> is a constant that minimizes the coefficient of variation.
although it might not be the original intention of the authors  <cit> , a test statistic like di can be motivated by a model with a mixture structure on gene standard deviations.
we begin with a simple case where xij ~ n and yik ~ n, and the null hypothesis is h <dig> : Δi =  <dig>  given δi, the lr test statistic is
we assume a mixture structure on gene standard deviations, where δi = σi with probability pv and δi = σ <dig> with probability  <dig> - pv.
we can then approximate  by
replacing σi with si and  with s <dig>  we can see that di and  only differ by a factor of 1/p, which does not change the ordering of test statistics.
the above derivation suggests that the sam statistic can also be considered an approximation to the lr test statistic under the mixture variance  assumption.
we can also write di as a weighted sum of ti and fci:
recall that under the mixture variance assumption, the lr test statistic is ri = wih <dig> + h <dig>  where h <dig> and h <dig> are both monotonic increasing functions.
both di and ri define rejection regions that are analogous to the union of the rejection regions defined by t test and fold change.
in other words, the sam procedure rejects h <dig> for large |ti|, or large |fci|, or both.
the sam statistic is a better approximation to the lr test statistic than the double filtering procedure.
as a side note, cui et al.
<cit>  proposed a shrunken t test procedure based on a variance estimator that borrow information across genes using the james-stein-lindley shrinkage concept.
this variance estimator shrinks individual variances toward a common value, which conceptually serves the same purpose as the mixture variance model.
from this perspective, we also consider the shrunken t statistic an approximation to the lr test statistic.
RESULTS
simulation study
we conducted a simulation study to compare the double filtering procedure to the shrinkage methods.
the simulation truth is specified as follows.
we tested  <dig> genes with  <dig> genes being truly de.
without loss of generality, we set μi =  <dig>  we further assumed
and
three scenarios were considered.
scenario 1: 90% of the genes with gene-specific variances and 10% of the genes with a common variance, and  <dig> replicates per gene under each condition.
scenario 2: same as scenario  <dig>  but with  <dig> replicates per gene under each condition.
scenario 3: all the genes having a gene-specific variance, and  <dig> replicates per gene under each condition.
for each scenario we repeated the simulation  <dig> times.
for the bayesian mixture model, we specified noninformative priors so that the posterior inference is dominated by information from data.
we let  =  =  <dig>  where  <dig>  is sufficiently large for expression levels on the logarithm scale.
to specify the hyper-parameters for the inverse gamma priors, first we set aσ = a <dig> =  <dig>  so that the inverse gamma priors have an infinite variance.
then we set the prior means,  and , equal to the average of the sample variances to solve for bσ and b <dig>  finally, we chose ar = br = av = bv =  <dig>  which corresponds to uniform priors for pr and pv.
five test statistics were compared: the marginal posterior probability  of a gene being de based on the bayesian mixture model, the sam statistic, the shrunken t statistic, the t statistic, and the double filtering with t statistic and fold change greater than  <dig>  the first three graphs in figure  <dig> plot the fdr versus the total number of selected genes under the above three scenarios.
the shrinkage methods  have comparable performance.
the double filtering procedure performs better than the traditional t statistic, but it is obviously outperformed by the three shrinkage methods.
we have tried different fold change cutoff values for the double filtering procedure  and the results did not change materially.
given the same number of selected genes, the shrinkage methods can identify more truly de genes than the double filtering procedure.
note that under the gene-specific variance assumption , the t test, which theoretically is the most powerful likelihood ratio test, still performs the poorest.
this result indicates the usefulness of shrinkage in microarray studies, where only a small number of replicates are available for each gene.
in short, the simulation study shows that for a given number of selected genes, well constructed shrinkage methods can outperform the double filtering procedure.
in scenario  <dig> and  <dig> of the simulation study, the true variance distribution is specified as the mixture of a point mass and an inverse gamma distribution, which might lead to a result that is biased in favor of a shrinkage method.
here we conduct another simulation study with a "real" variance distribution, denoted as scenario  <dig>  specifically, let xij  and yik  be the observed expression levels from a real microarray study.
define the residual vector ei =' by
then ei can be considered a set of random errors sampled based on the true variance distribution.
we simulate  <dig> data sets according to the following steps.
for iteration s  and gene i ,
<dig>  obtain a random permutation of , denoted by ;
<dig>  generate  as described in the previous simulation scenarios;
<dig>  for j =  <dig>  ... , m0i, compute  = , and for k =  <dig>  ... , m1i, compute , where  is the jth element of .
the real data comes from a microarray study comparing the gene expressions of breast cancer tumors with brca <dig> mutations, brca <dig> mutations, and sporadic tumors  <cit> .
the data set is available at http://research.nhgri.nih.gov/microarray/nejm_supplement.
here we only consider the brca <dig> group and the brca <dig> group.
there are  <dig> genes, with  <dig> arrays in the brca <dig> group and  <dig> arrays in the brca <dig> group.
we analyzed the data on the log <dig> scale.
following storey and tibshirani  <cit> , we eliminated genes with aberrantly large expression values , which left us with measurements on n =  <dig> genes.
the fourth graph in figure  <dig> compares the different methods under scenario  <dig>  where the residual vector ei was constructed based on the breast cancer data.
we kept the same replicate number in the experiment, with  <dig> replicates per gene in one group and  <dig> replicates in the other group.
the relative performance of the five methods remains unchanged as in the other scenarios.
in current microarray studies, the number of replicates per gene can be easily  <dig> or more due to the low cost of array and the easiness to collect patients.
so we considered two scenarios with a relatively large number of replicates.
scenario 5: 90% of the genes with gene-specific variances and 10% of the genes with a common variance, and  <dig> replicates per gene under each condition.
scenario 6: all the genes having a gene-specific variance, and  <dig> replicates per gene under each condition.
in each of the two scenarios, we assume there are  <dig> genes with  <dig> genes being truly de.
the two graphs in the bottom panel of figure  <dig> plot the fdr versus the total number of selected genes for the five test statistics under scenario  <dig> and scenario  <dig>  respectively.
the comparison demonstrates that when the replicate number is large, the performance of the traditional t test becomes comparable to the performance of the shrinkage methods, thanks to the more reliable estimate of gene variance component.
more importantly, the drawback of the double filtering procedure becomes more obvious, which has substantially worse performance compared to the other methods, including the t test.
experimental datasets
in this section we compared the shrinkage methods with the double filtering procedure based on two microarray datasets.
the first is the golden spike data  <cit>  where the identities of truly de genes are known.
the golden spike dataset includes two conditions, with  <dig> replicates per condition.
each array has  <dig>  probesets, among which  <dig>  have non-spiked-in rnas,  <dig>  have equal concentrations of rnas, and  <dig>  are spiked-in at different fold-change levels, ranging from  <dig>  to 4-fold.
compared with other spike datasets, the golden spike dataset has a larger number of probsets that are known to be de, making it popular for comparing performance among different methods.
irizarry et al.
<cit>  pointed out that "the feature intensities for genes spiked-in to be at 1: <dig> ratios behave very differently from the features from non-spiked-in genes".
following opgen-rhein and strimmer  <cit> , we removed the  <dig>  probe sets for spike-ins with ratio 1: <dig> from the original data, leaving in total  <dig>  genes and  <dig>  known de genes.
figure  <dig> plots the fdr under each testing procedures versus the total number of rejections.
for the double filtering procedure, the fold change cutoff was set at  <dig>  because only  <dig> genes have a fold change greater than  <dig> .
the figure indicates that the shrinkage methods  have similar performance, and they outperform the double filtering procedure and t test.
the second is the breast cancer dataset  <cit>  described in the simulation study.
with the identities of truly de genes unknown, we estimated the fdr for the sam test, the shrunken t test, the t test, and the double filtering procedure, through the permutation approach described in  <cit> .
for bayesian methods, newton et al.
<cit>  proposed to compute the bayesian fdr, which is the posterior proportion of false positives relative to the total number of rejections.
however, the bayesian fdr is incomparable to the permutation-based fdr estimate employed by frequentist methods  <cit> .
cao and zhang  <cit>  developed a generic approach to estimating the fdr for bayesian methods under the permutation-based framework.
a computationally efficient algorithm was developed to approximate the null distribution of the bayesian test statistic, the posterior probability.
the approach can provide an unbiased estimate of the true fdr.
constructed under the same permutation-based framework, the resulting fdr estimate allows a fair comparison between full bayesian methods with other testing procedures.
we adopted the approach in  <cit>  to estimate the fdr of the bayesian mixture model .
figure  <dig> plots the permutation-based fdr estimates under each testing procedure versus the total number of rejections.
it shows that the shrinkage methods can considerably outperform the double filtering procedure.
CONCLUSIONS
it has been a common practice in microarray analysis to use fold change and t statistic to double filter de genes.
in this paper, we provided a close examination on the drawback of the double filtering procedure, where fold change and t statistic are based on contradicting assumptions.
we further demonstrated that several shrinkage methods  can be united under the mixture gene variance assumption.
based on the theoretical derivation, the simulation study, and the real data analysis, we showed compelling evidence that well constructed shrinkage methods can outperform the double filtering procedure in identifying de genes.
with publicly available softwares, these methods are as easy to implement as the double filtering procedure.
we acknowledge some researchers' argument that the double filtering procedure might work well because it filters out the genes that show relatively small differences between conditions, which are sometimes considered to be less biologically meaningful.
this argument, however, is based on the criterion of so called "biological meaningfulness" instead of testing power.
although many biologists refer to fold change in terms of "biological meaningfulness", there is in fact no clear cut-off for it, and 2-fold is often invoked merely based on convenience.
in addition, different normalization methods can differ quite drastically in terms of the fold changes they produce.
so a particular cut-off in fold change could mean one thing using one method and quite another using a different method.
taken together, even if researchers decide to employ the double filtering procedure based on the rationale of "biological meaningfulness", it is still helpful to understand the potential loss in power.
authors' contributions
sz and jc conceived the study, conducted the examination on the double filtering procedure, analyzed the data, and drafted the paper.
all authors read and approved the final manuscript.

 with the development of microarray technologies, researchers now can measure the relative expressions of tens of thousands of genes simultaneously. however, the number of replicates per gene is usually small, far less than the number of genes. many statistical methods have been developed to identify differentially expressed genes. the use of fold change is among the first practice. it can be inefficient and erroneous because of the additional uncertainty induced by dividing two intensity values. there are variants of student's t test procedure that conduct a test on each individual gene and then correct for multiple comparisons. the problem is, with a large number of tests and a small number of replicates per gene, the statistics are very unstable. for example, a large t statistic might arise because of an extremely small variance, even with a minor difference in the expression. the disadvantage of fold-change approach and t test has been pointed out by a number of authors. there are approaches proposed to improve estimation of gene variances by borrowing strength across genes. despite the flaw, fold change and t test are the most intuitive approaches and they have been applied widely in practice. to control the error rate, many researchers use fold change and t test together, hoping that the double filtering will provide extra confidence in the test results. specifically, a gene is flagged as de only if the p-value from t test is smaller than a certain threshold and the fold change is greater than a cutoff value. for example, in, genes were found to be de with two cutoff values. there are numerous examples in the literature that implement the double filtering procedure with fold change and t statistic. we argue, however, that the double filtering procedure provides higher confidence mainly because it produces a shorter list of selected genes. given the same number of genes selected, a well constructed shrinkage test can significantly outperform the double filtering method. fold change takes the ratio of a gene's average expression levels under two conditions. it is usually calculated as the difference on the log scale. let xij be the log-transformed expression measurement of the ith gene on the jth array under the control, and yik be the log-transformed expression measurement of the ith gene on the kth array under the treatment. we define and. fold change is computed by as for the traditional t test, it is usually calculated on the log scale to adjust for the skewness in the original gene expression measurements. the t statistic is then computed by where is the pooled variance of xij and yik. comparing and, it is obvious that fold change and t statistic are based on two contradicting assumptions. the underlying assumption of fold change is that all genes share a common variance, which is implied by the omission of the variance component in. on the other hand, the inclusion of in suggests that t test assumes gene-specific variances. in order for a gene to be flagged as de, the double filtering procedure would require the gene to have extreme test scores under the common variance assumption as well as under the gene-specific variance assumption. it is analogous to using the intersection of the rejection regions defined by fold change and t statistic. assuming a common variance for all the genes apparently is an oversimplification. the assumption of gene-specific variances, however, leads to unstable estimates due to limited replicates from each gene. a more realistic assumption might lie in between the two extremes, i.e., modeling gene variances by a mixture of two components, one being a point mass at the common variance, another being a continuous distribution for the gene-specific variances. under this mixture variance assumption, a de gene could have a large fold change or a large t statistic, but not necessarily both. taking intersection of the rejection regions flagged by fold change and t statistic, as is adopted by the double filtering procedure, might not be the best strategy under the mixture variance assumption. the goal of the paper is not to propose a new testing procedure in microarray analysis, but to provide insight on the drawback of the widely used double filtering procedure with fold change and t test. we present a theoretically most powerful likelihood ratio test under the mixture variance assumption. we further demonstrate that two shrinkage test statistics, one from the bayesian model and the other from the significance analysis of microarrays test, can be united as approximations to the lr test. this association explains why those shrinkage methods can considerably outperform the double filtering procedure. a simulation study and real microarray data analyses are then presented to compare the shrinkage tests and the double filtering procedure. methods a likelihood ratio test for gene i, we use fi = pvfi + fi a mixture of two components fi and fi to denote the density under the null hypothesis that the gene is not de under two experiment conditions. density fi is defined under the gene-specific variance assumption, fi is defined under the common variance assumption, and pv is the mixing probability. similarly, we use gi = pvgi + gi to denote the density under the alternative hypothesis, with gi and gi defined in a similar fashion as fi and fi for example, in the context of testing de genes, we can assume fi = n, fi = n, gi = n, and gi = n, where is the assumed common variance, is the gene-specific variance, μi is the mean expression level under the control, and Δi is the difference in the expression levels between two conditions. under the null hypothesis h: Δi = the likelihood ratio test statistic, which is the most powerful among all test statistics, is the ri statistic is a weighted sum of two ratios gi1/fi and gi2/fi with weight wi = pvfi1/. under the normality assumption, it is easy to show that ri = wih + h where fci and ti are fold change and t statistic, as defined in and. both h and h are monotonic increasing functions. the rejection region of the lr test is defined by ri >λr, where λr is the threshold to attain a certain test size. in order to reject h it requires that either |fci| is large, or |ti| is large, or both. in this sense, the lr test rejection region is more like a union of the rejection regions defined by fold change and t statistic. on the other hand, the double filtering procedure with fold change and t statistic would require both |fci| and |ti| to be large. this practice is analogous to using the intersection of the two rejections regions determined by |fci| and |ti|. compared with the lr test, the double filtering procedure will lose power. the "loss of power" has two meanings. first, for a given false discovery rate, the double filtering procedure produces a shorter list of identified genes for further investigation. second, for a given number of identified genes, the list produced by the double filtering procedure has a higher fdr. the double filtering procedure offers a false sense of confidence by producing a shorter list. the lr test statistic ri requires one to know the true values of parameters p, μi,,, and Δi, which are usually unknown in reality. one strategy is to estimate ri by, where the maximum likelihood estimates of the unknown parameters are plugged into. unfortunately, with a small number of replicates from each gene, the mle would be extremely unstable and lead to unsatisfactory testing results. a bayesian model was constructed under the mixture variance assumption to detect de genes. the inference is made based on the marginal posterior probability of a gene being de, denoted by zi = p. here x = {xij} and y = {yik} are the collection of gene expression data under the two conditions. we will show that, like, zi is also an approximation to ri. the difference between and zi is that the former plugs in the point estimates of unknown parameters, while the latter marginalizes the unknown parameters with respect to their posterior distribution. in the bayesian inference, the uncertainty from various sources are accounted for in a probabilistic fashion. similar to the bayesian mixture model, some existing methods also try to strike a balance between the two extreme assumptions of a common variance and gene-specific variances. the sam statistic slightly modifies the t-statistic by adding a constant to the estimated gene-specific standard deviation in the denominator. we will present it as being motivated by a mixture model on the variances. furthermore, the sam statistic can be directly written as a weighted sum of t statistic and fold change. thus both the bayesian method and the sam method are approximations to the lr test under the mixture variance assumption, and they can achieve better performance than the double filtering procedure. the bayesian mixture model cao et al. proposed a bayesian mixture model to identify de genes, which has shown comparable performance to frequentist shrinkage methods. with parameters defined similarly as in the lr test, gene expression measurements xij and yij are modeled by normal distributions with a mixture structure on the variances, a latent variable ri is used to model the expression status of the ith gene,, where ri = 0/ indicates that gene i is non-de/de and it is modeled by a bernoulli distribution: ri | pr ~ bernoulli. for and, it is assumed that and ~ ig. here ig denotes an inverse gamma distribution with mean b/. the other hyper-priors include, μi ~ n, pr ~ u, and pv ~ u. more details can be found in. we make inference based on zi = p = p, the marginal posterior probability that gene i is de. a gene is flagged as de if zi >λz, where λz is a certain cutoff. we argue that the bayesian rejection region defined by zi >λz is an approximation to the lr test rejection region defined by ri >λr. first we have here p is the joint posterior distribution of, marginalized with respect to other random parameters. it is easy to show that given parameters, p is an increasing function of ri. rejecting h for ri >λr is equivalent to rejecting for p >λz, with λz = λr/. thus the two test statistics, p and ri, are equivalent. expression demonstrates that zi is obtained from p by integrating with respect to the unknown parameters under the joint posterior distribution. if the integral does not have a closed form, we can conduct numerical integration to calculate zi through the gibbs sampling algorithm. the uncertainty from those unknown parameters are accounted for in a probabilistic fashion. it is in this sense that we consider zi a good approximation to the lr test statistic ri. the sam test the sam statistic is defined as where si is the gene-specific standard deviation, and s is a constant that minimizes the coefficient of variation. although it might not be the original intention of the authors, a test statistic like di can be motivated by a model with a mixture structure on gene standard deviations. we begin with a simple case where xij ~ n and yik ~ n, and the null hypothesis is h: Δi = given δi, the lr test statistic is we assume a mixture structure on gene standard deviations, where δi = σi with probability pv and δi = σ with probability - pv. we can then approximate by replacing σi with si and with s we can see that di and only differ by a factor of 1/p, which does not change the ordering of test statistics. the above derivation suggests that the sam statistic can also be considered an approximation to the lr test statistic under the mixture variance assumption. we can also write di as a weighted sum of ti and fci: recall that under the mixture variance assumption, the lr test statistic is ri = wih + h where h and h are both monotonic increasing functions. both di and ri define rejection regions that are analogous to the union of the rejection regions defined by t test and fold change. in other words, the sam procedure rejects h for large |ti|, or large |fci|, or both. the sam statistic is a better approximation to the lr test statistic than the double filtering procedure. as a side note, cui et al. proposed a shrunken t test procedure based on a variance estimator that borrow information across genes using the james-stein-lindley shrinkage concept. this variance estimator shrinks individual variances toward a common value, which conceptually serves the same purpose as the mixture variance model. from this perspective, we also consider the shrunken t statistic an approximation to the lr test statistic.  simulation study we conducted a simulation study to compare the double filtering procedure to the shrinkage methods. the simulation truth is specified as follows. we tested genes with genes being truly de. without loss of generality, we set μi = we further assumed and three scenarios were considered. scenario 1: 90% of the genes with gene-specific variances and 10% of the genes with a common variance, and replicates per gene under each condition. scenario 2: same as scenario but with replicates per gene under each condition. scenario 3: all the genes having a gene-specific variance, and replicates per gene under each condition. for each scenario we repeated the simulation times. for the bayesian mixture model, we specified noninformative priors so that the posterior inference is dominated by information from data. we let = = where is sufficiently large for expression levels on the logarithm scale. to specify the hyper-parameters for the inverse gamma priors, first we set aσ = a = so that the inverse gamma priors have an infinite variance. then we set the prior means, and, equal to the average of the sample variances to solve for bσ and b finally, we chose ar = br = av = bv = which corresponds to uniform priors for pr and pv. five test statistics were compared: the marginal posterior probability of a gene being de based on the bayesian mixture model, the sam statistic, the shrunken t statistic, the t statistic, and the double filtering with t statistic and fold change greater than the first three graphs in figure plot the fdr versus the total number of selected genes under the above three scenarios. the shrinkage methods have comparable performance. the double filtering procedure performs better than the traditional t statistic, but it is obviously outperformed by the three shrinkage methods. we have tried different fold change cutoff values for the double filtering procedure and the results did not change materially. given the same number of selected genes, the shrinkage methods can identify more truly de genes than the double filtering procedure. note that under the gene-specific variance assumption, the t test, which theoretically is the most powerful likelihood ratio test, still performs the poorest. this result indicates the usefulness of shrinkage in microarray studies, where only a small number of replicates are available for each gene. in short, the simulation study shows that for a given number of selected genes, well constructed shrinkage methods can outperform the double filtering procedure. in scenario and of the simulation study, the true variance distribution is specified as the mixture of a point mass and an inverse gamma distribution, which might lead to a result that is biased in favor of a shrinkage method. here we conduct another simulation study with a "real" variance distribution, denoted as scenario specifically, let xij and yik be the observed expression levels from a real microarray study. define the residual vector ei =' by then ei can be considered a set of random errors sampled based on the true variance distribution. we simulate data sets according to the following steps. for iteration s and gene i, obtain a random permutation of, denoted by; generate as described in the previous simulation scenarios; for j =..., m0i, compute =, and for k =..., m1i, compute, where is the jth element of. the real data comes from a microarray study comparing the gene expressions of breast cancer tumors with brca mutations, brca mutations, and sporadic tumors. the data set is available at http://research.nhgri.nih.gov/microarray/nejm_supplement. here we only consider the brca group and the brca group. there are genes, with arrays in the brca group and arrays in the brca group. we analyzed the data on the log scale. following storey and tibshirani, we eliminated genes with aberrantly large expression values, which left us with measurements on n = genes. the fourth graph in figure compares the different methods under scenario where the residual vector ei was constructed based on the breast cancer data. we kept the same replicate number in the experiment, with replicates per gene in one group and replicates in the other group. the relative performance of the five methods remains unchanged as in the other scenarios. in current microarray studies, the number of replicates per gene can be easily or more due to the low cost of array and the easiness to collect patients. so we considered two scenarios with a relatively large number of replicates. scenario 5: 90% of the genes with gene-specific variances and 10% of the genes with a common variance, and replicates per gene under each condition. scenario 6: all the genes having a gene-specific variance, and replicates per gene under each condition. in each of the two scenarios, we assume there are genes with genes being truly de. the two graphs in the bottom panel of figure plot the fdr versus the total number of selected genes for the five test statistics under scenario and scenario respectively. the comparison demonstrates that when the replicate number is large, the performance of the traditional t test becomes comparable to the performance of the shrinkage methods, thanks to the more reliable estimate of gene variance component. more importantly, the drawback of the double filtering procedure becomes more obvious, which has substantially worse performance compared to the other methods, including the t test. experimental datasets in this section we compared the shrinkage methods with the double filtering procedure based on two microarray datasets. the first is the golden spike data where the identities of truly de genes are known. the golden spike dataset includes two conditions, with replicates per condition. each array has probesets, among which have non-spiked-in rnas, have equal concentrations of rnas, and are spiked-in at different fold-change levels, ranging from to 4-fold. compared with other spike datasets, the golden spike dataset has a larger number of probsets that are known to be de, making it popular for comparing performance among different methods. irizarry et al. pointed out that "the feature intensities for genes spiked-in to be at 1: ratios behave very differently from the features from non-spiked-in genes". following opgen-rhein and strimmer, we removed the probe sets for spike-ins with ratio 1: from the original data, leaving in total genes and known de genes. figure plots the fdr under each testing procedures versus the total number of rejections. for the double filtering procedure, the fold change cutoff was set at because only genes have a fold change greater than. the figure indicates that the shrinkage methods have similar performance, and they outperform the double filtering procedure and t test. the second is the breast cancer dataset described in the simulation study. with the identities of truly de genes unknown, we estimated the fdr for the sam test, the shrunken t test, the t test, and the double filtering procedure, through the permutation approach described in. for bayesian methods, newton et al. proposed to compute the bayesian fdr, which is the posterior proportion of false positives relative to the total number of rejections. however, the bayesian fdr is incomparable to the permutation-based fdr estimate employed by frequentist methods. cao and zhang developed a generic approach to estimating the fdr for bayesian methods under the permutation-based framework. a computationally efficient algorithm was developed to approximate the null distribution of the bayesian test statistic, the posterior probability. the approach can provide an unbiased estimate of the true fdr. constructed under the same permutation-based framework, the resulting fdr estimate allows a fair comparison between full bayesian methods with other testing procedures. we adopted the approach in to estimate the fdr of the bayesian mixture model. figure plots the permutation-based fdr estimates under each testing procedure versus the total number of rejections. it shows that the shrinkage methods can considerably outperform the double filtering procedure.  it has been a common practice in microarray analysis to use fold change and t statistic to double filter de genes. in this paper, we provided a close examination on the drawback of the double filtering procedure, where fold change and t statistic are based on contradicting assumptions. we further demonstrated that several shrinkage methods can be united under the mixture gene variance assumption. based on the theoretical derivation, the simulation study, and the real data analysis, we showed compelling evidence that well constructed shrinkage methods can outperform the double filtering procedure in identifying de genes. with publicly available softwares, these methods are as easy to implement as the double filtering procedure. we acknowledge some researchers' argument that the double filtering procedure might work well because it filters out the genes that show relatively small differences between conditions, which are sometimes considered to be less biologically meaningful. this argument, however, is based on the criterion of so called "biological meaningfulness" instead of testing power. although many biologists refer to fold change in terms of "biological meaningfulness", there is in fact no clear cut-off for it, and 2-fold is often invoked merely based on convenience. in addition, different normalization methods can differ quite drastically in terms of the fold changes they produce. so a particular cut-off in fold change could mean one thing using one method and quite another using a different method. taken together, even if researchers decide to employ the double filtering procedure based on the rationale of "biological meaningfulness", it is still helpful to understand the potential loss in power. authors' contributions sz and jc conceived the study, conducted the examination on the double filtering procedure, analyzed the data, and drafted the paper. all authors read and approved the final manuscript. 
many researchers use the double filtering procedure with fold change and t test to identify differentially expressed genes, in the hope that the double filtering will provide extra confidence in the results. we further demonstrate that the posterior inference based on a bayesian mixture model and the widely used significance analysis of microarrays statistic are better approximations to the likelihood ratio test than the double filtering procedure. we show that fold change assumes all genes to have a common variance while t statistic assumes gene-specific variances. under the assumption that gene variances arise from a mixture of a common variance and gene-specific variances, we develop the theoretically most powerful likelihood ratio test statistic. we demonstrate through hypothesis testing theory, simulation studies and real data examples, that well constructed shrinkage testing methods, which can be united under the mixture gene variance assumption, can considerably outperform the double filtering procedure.
 many researchers use the double filtering procedure with fold change and t test to identify differentially expressed genes, in the hope that the double filtering will provide extra confidence in the results. due to its simplicity, the double filtering procedure has been popular with applied researchers despite the development of more sophisticated methods.  this paper, for the first time to our knowledge, provides theoretical insight on the drawback of the double filtering procedure. we show that fold change assumes all genes to have a common variance while t statistic assumes gene-specific variances. the two statistics are based on contradicting assumptions. under the assumption that gene variances arise from a mixture of a common variance and gene-specific variances, we develop the theoretically most powerful likelihood ratio test statistic. we further demonstrate that the posterior inference based on a bayesian mixture model and the widely used significance analysis of microarrays statistic are better approximations to the likelihood ratio test than the double filtering procedure.  we demonstrate through hypothesis testing theory, simulation studies and real data examples, that well constructed shrinkage testing methods, which can be united under the mixture gene variance assumption, can considerably outperform the double filtering procedure. 