BACKGROUND
microarray technology has become a standard experimental method in bio-medical research.
in the analysis of microarray data, one of the most fundamental tasks is the identification of differentially expressed genes while controlling false positives and minimizing false negatives.
this is a multiple hypothesis test problem which analyzes thousands or tens of thousands of genes simultaneously.
in these tests we often need to control the false discovery among the rejected hypotheses under a pre-specified level while maintaining maximal power.
thus, there is a trade off in the control of the type-i error between rejecting true null hypotheses  versus accepting true alternative hypotheses .
traditional bonferroni correction procedures are designed to control the family wise error rate , which guards against making one or more type i errors among a family of hypothesis tests.
however, these procedures may be excessively conservative for microarray analysis where the number of hypotheses is very large and a substantial fraction of the genes are differentially expressed  <cit> .
a more appropriate approach is to control the false discovery rate , which is the proportion of type i errors among all rejected hypotheses  <cit> .
this approach is particularly useful in exploratory analyses, where the objective is to maximize the discovery of true positives, rather than guarding against one or more false positive results.
a number of methods have been proposed to control the fdr given a population of hypothesis tests.
these methods usually assume that the distribution of the test statistics, f, can be modeled by a mixture of two components  <cit> :
f = π0f <dig> + f <dig>     π <dig> = m0/m
where f <dig> is the distribution of the test statistics under h <dig>  which by definition equals to  <dig> when using p-values when tests are independent, f <dig> is the distribution of the test statistics under h <dig>  m <dig> is the number of true h <dig>  m is the total number of hypotheses under consideration, and π <dig> is the proportion of true h <dig>  the methods proposed by benjamini et al  <cit>  to control fdr do not estimate π0; therefore, they provide the strongest controls on fdr but have the lowest power compared to other methods that do so.
in many actual applications where a considerable number of genes are differentially expressed, assuming π <dig> =  <dig> may be too conservative causing loss of power.
several alternative methods, such as nonparametric empirical bayesian pfdr criterion and its p-value equivalent called q-value method  <cit> , bin-wise model  <cit> , local fdr method  <cit> , parametric beta-uniform mixture models  <cit> , the lowest slope estimator   <cit> , the spacing loess histogram  method  <cit> , the nonparametric mle method  <cit> , the moment generating function approach  <cit> , and the poisson regression approach  <cit> , have all been proposed to estimate π <dig> by pooling test statistics and controlling fdr based on the estimated π <dig>
in these methods, one of the critical steps is estimating the proportion of null hypotheses, π <dig>  when using p-values, these estimations usually depend on the assumption that f <dig> follows a uniform distribution.
this assumption, which is of critical importance for the methods of statistical inference that employ pooling test statistics across genes  <cit> , is valid when all test hypotheses are independent and identically distributed.
furthermore, when there are only weak correlations, or "clumpy" correlations , the uniform assumption is not strongly violated and the method remains adequate.
however, in datasets with large scale strong correlations, the joint distribution of the test statistics will no longer be the product of marginal distributions, and the observed f <dig> will severely deviate from uniform, causing the current π <dig> estimation methods to become very unstable.
increased variation and bias of π <dig>  as well as fdr, was also observed by wu et al  <cit>  in datasets with strong local correlations.
the effect of correlation on simultaneous significance tests was previously discussed theoretically  <cit> , and a number of permutation based fdr control methods were proposed, such as sam  <cit> , dchip  <cit> , ge et al  <cit> , meinshausen et al  <cit>  and efron  <cit> .
in these methods, the distribution of f <dig> was modeled empirically through permutations, which naturally considered the correlation.
however, like benjamini et al  <cit> , these methods don't estimate π0; therefore, in datasets with a large number of differentially expressed genes, the fdr control may be overly conservative with a loss of power.
therefore we proposed  <dig> re-sampling schemes, similar to model averaging in bagging methods, to reduce the variation in estimating π <dig> in datasets with strong correlation between gene expression values.
our methods produced a more stable and conservative estimation of π <dig> and, therefore, provided stronger control of false discovery rate with only a minor sacrifice of power.
implementation
creating simulated data set
to test the performance of various algorithms in estimating π <dig>  we generated  <dig> types of simulated datasets.
both datasets had strong correlation between subsets of genes and a known proportion of true null hypotheses, to represent the log transformed microarray expression data.
data-b
the first simulation method is similar to the method used in qiu et al  <cit>  and wu et al  <cit> .
assume n samples and m genes, with n/ <dig> samples per class.
the m genes were divided into k blocks with each block consisting of m/k genes.
assume independence between blocks and constant correlation coefficient between genes within each block.
for block l and sample j, we first created a block center vector
blj = dl·xj + δlj, l =  <dig> ...,k, j =  <dig> ...,n
where dl was the mean difference between groups, it equals to  <dig>  with probability π <dig>  or was generated from beta distribution with parameter  otherwise; xj was a group indicator; and δlj was i.i.d.
n to represent sample specific noise.
then the expression value of gene i in sample j in that block, ylji, was generated by
ylji = ρ·blj + ·elji, i =  <dig> ...,m/k
where ρ was the correlation constant which takes value between  <cit>  determining the correlation coefficient between genes within the block, and elji was i.i.d.
n. the k blocks of genes were generated independently of each other and then pooled to form the whole dataset.
we call this type of dataset which contains blocks of highly correlated genes data-b in our experiments.
data-m
the above data-b model is over-simplified in many aspects, and is still a "clumpy" structure, although the clumpiness can be pronounced.
to mimic more realistic situations, we generated a second type of simulated data based on an actual human breast cancer microarray dataset  <cit>  obtained with affymetrix u <dig> plus  <dig>  microarrays.
the dataset contains  <dig> estrogen receptor positive  cases and  <dig> estrogen receptor negative  cases.
the data were normalized by the gcrma algorithm <cit> , and the gene  expression levels were log2-transformed.
according to published literatures  <cit> , the er status is one of the most predominant partitioning factors for molecular classification of breast cancer.
we therefore took some of the genes differentially expressed between the two classes as "truly h1" genes.
we selected  <dig> genes with differences in mean of log transformed expression levels between the two classes greater than  <dig>  .
from these genes, each time we randomly chose  <dig> to form our simulated dataset and then randomly picked π <dig> proportion of the  <dig> genes to establish h <dig> genes by scrambling these genes together between samples.
thus, we obtained a simulated dataset with a known number of null hypotheses while the correlations among both the differentially and non-differentially expressed genes were maintained.
we call this type of dataset data-m in our experiments  <cit> .
estimating π <dig> by re-sampling strategy
to get a better estimation of π <dig> in datasets with strong correlations, we proposed  <dig> re-sampling methods to replace the original π <dig> estimation step in q-value method which estimates π <dig> directly from the p-value distribution  <cit> .
the first method, termed samps, re-samples without replacement 2/ <dig> of the samples from each class, calculates p-values for each gene, and then estimates the π <dig> from the p-values.
for each dataset, we performed this procedure  <dig> times and used the upper quartile to replace the π <dig> estimated by the q-value method.
the second method, termed sampg, re-samples without replacement 2/ <dig> of the values from each class for each gene independently, calculates the p-values for each gene, and then estimates the π <dig>  we also performed this procedure  <dig> times for each dataset and used the upper quartile to replace the π <dig> estimated by the q-value method.
after each re-sampling step, we have to feed the p-values into a π <dig> estimator.
this estimator could be the π <dig> estimation by q-value method, the mgf method  <cit> , or any other unbiased π <dig> estimators.
in this paper, we tried to use both q-value and mgf as the plug-in estimator, and called them samps.q and sampg.q when using q-value as the plug-in estimator, or samps.m and sampg.m when using mgf as the plug-in estimator.
RESULTS
variation of the π <dig> estimation when genes are correlated
to evaluate the impact of correlation structure in large datasets on the estimation of the proportion of true null hypotheses , we first evaluated current methods on a published microarray dataset  <cit> .
this dataset consists of  <dig> estrogen receptor positive  and  <dig> estrogen receptor negative  breast cancers.
the gene expression levels were normalized by the gcrma algorithm  <cit>  and log <dig> transformed.
expression values were filtered to eliminate low expressing genes with mean expression below  <dig> and constant expressing genes with coefficient of variation below  <dig> .
a total of  <dig>  genes passed this filtering criterion.
we bootstrapped  <dig> datasets from this microarray data and used the q-value  <cit>  and twilight  <cit>  methods to estimate the proportion of null hypotheses.
the π <dig> estimates were similar by the  <dig> methods on each bootstrapped data set; however, both methods showed a large range of π <dig> among the bootstrapped datasets that varied from  <dig>  to  <dig>  .
to further investigate the influence of gene correlation on the estimation of π <dig>  we generated simulated dataset data-b which contained blocks of highly correlated genes, but all genes were non-differentially expressed.
in our simulation we set the correlation constant ρ equal to  <dig> , and created datasets containing  <dig>   <dig>  or  <dig>  genes with  <dig>   <dig>   <dig> or  <dig>  genes per block.
after calculating the p-value for each gene, we calculated the coefficient of variation  of the bar heights in the histogram of p-values by splitting the p-values into  <dig> bins between  <cit>  with the width of each bin being  <dig> .
figure  <dig> shows the histogram of p-values in one of the simulations having  <dig>  genes all independent with each other, and figure  <dig> shows the histogram of p-values in another simulation having  <dig>  genes with  <dig>  genes per block.
comparing figure  <dig> and  <dig> we can see that with highly correlated genes, the distribution of p-value deviated significantly from uniform, although none of the genes were differentially expressed.
we created  <dig> simulations for each type of data, calculated their cv of the bar height in the p-value histogram, and plotted the box plots of the cvs in figure  <dig>
in this simulated study, since all genes were non-differentially expressed, the p-values should follow a uniform distribution, and the histogram of p-values should be flat if genes were independent of each other.
when the number of correlated genes in each block was small, for example,  <dig>  or  <dig>  genes per block, the distribution of p-values approximated a uniform distribution and the cv of the histogram of p-values were close to expected under the independent assumption.
however, the cv of the histogram of p-values increased significantly with the growth of correlation structure.
in other words, although none of the genes were differentially expressed, the distribution of p-value deviated increasingly from a uniform distribution when large groups of genes were correlated.
we also calculated the cv of the histogram of p-values in our microarray dataset by randomly permuting the class labels, which makes all genes non-differentially expressed but still correlated, as well as randomly permuting all expression values across genes and samples which makes genes non-differentially expressed and also independent.
the results of  <dig> permutations showed a dramatically higher cv for the p-value histogram of datasets with only class labels permuted but with gene-gene correlations intact .
comparing figure  <dig>  figure  <dig> and figure  <dig>  it is apparent that increased correlation among genes greatly increased the deviation of p-value distribution from uniform.
therefore, strong correlation structures will increase the variation in estimated π <dig>  and inevitably subsequent statistical inferences and false discovery control procedures would be influenced by this unstable π <dig> estimation.
improving the π <dig> estimation by re-sampling strategies
to improve the estimation of π <dig> in datasets with strong correlations, we proposed two methods, termed samps and sampg, to replace the original π <dig> estimation step in the q-value method.
in the samps algorithm, we used a model averaging strategy.
we repeatedly sampled 2/ <dig> of the data from each class without replacement, calculated the p-values for genes, estimated π <dig> from the p-value distribution and finally used the upper quartile of the re-sampling estimations in the subsequent statistical inferences.
in the sampg algorithm, we further broke down the correlations between genes and stabilized the π <dig> estimation.
in this algorithm, we repeatedly sampled without replacement 2/ <dig> of the expression values from each class independently for each gene, calculated their p-values, estimate π <dig>  and then used the upper quartile of the re-sampling estimations in subsequent analysis.
for the choice of the plug-in π <dig> estimator, we tried to use both q-value  <cit>  and moment generating function  <cit>  methods, and called them sampg.q, samps.q and sampg.m, samps.m, respectively.
to test the variation in π <dig> estimation induced by strong correlation structures, and the performance of our proposed samps and sampg methods, we created simulated datasets data-b and data-m, with true π <dig> being  <dig> ,  <dig> ,  <dig> ,  <dig> , and the correlation constant ρ being  <dig> ,  <dig>  and  <dig> , respectively.
we created  <dig> datasets for each combination of these control parameters.
in data-b, we created  <dig> genes forming  <dig> blocks with  <dig> genes per block to simulate large scale correlation between genes.
in data-m we randomly selected  <dig> genes with differences in the mean of log transformed expression levels greater than  <dig> , then randomly scrambled 90%, 80%, 70% or 60% of them to create datasets with known proportion of true null hypothesis and strong gene-gene correlations.
we applied the samps and sampg strategies to estimate π <dig> for each of the simulated datasets, and compared our results to a number of other methods.
the π <dig> estimation methods we used and their corresponding r functions are:
<dig>  the lowest slope estimator   <cit>  with parameter determined via bootstrap   <cit> ; function fdr.estimate.eta0
<dig>  q-value method with tuning parameter chosen by smoother method   <cit> ; function qvalue
<dig>  the spacing loess histogram   <cit> ; function splosh
<dig>  the beta-uniform mixture model   <cit> ; function bum.mle
<dig>  the nonparametric mle method  <cit> ; function convest
<dig>  the poisson regression approach   <cit> ; function p <dig> mom
<dig>  the moment generating function   <cit> ; function p <dig> mom
<dig>  the lowest slope estimator   <cit>  with parameter determined adaptively   <cit> ; function fdr.estimate.eta0
<dig>  sampg method, with q-value plug in, 2nd quartile
<dig>  sampg method, with q-value plug in, 3rd quartile
<dig>  samps method, with q-value plug in, 2nd quartile
<dig>  samps method, with q-value plug in, 3rd quartile
<dig>  sampg method, with mgf plug in, 2nd quartile
<dig>  sampg method, with mgf plug in, 3rd quartile
<dig>  samps method, with mgf plug in, 2nd quartile
<dig>  samps method, with mgf plug in, 3rd quartile
the π <dig> estimations were shown as boxplots in figure  <dig>  and the mean square error  were listed in table  <dig>
we then selected genes with fdr controlled under  <dig>  level based on the π <dig> estimated by our method, calculated the actual false discovery rate and power, and compared our results to that of fdr controlling method proposed by benjamini et al.
with correlations considered   <cit> , the permutation based fdr controlling method   <cit> , and q-value method  <cit> .
the boxplot of actual false discovery rate and power on the  <dig> types of data-m simulations were shown in figure  <dig>  from figure  <dig> it can be seen that both the by and howmany methods provided strong control of fdr, their actual fdr level was much lower than  <dig> , and their power to detect true alternatives were much lower than methods where the proportion of true null hypotheses was estimated and used in the fdr control.
comparing the q-value method and our proposed samps and sampg method, for majority of the cases the actual fdr level were still controlled below  <dig>  level, although some outliers exist with actual fdr up to  <dig> ~ <dig> , due to the unstable pi_ <dig> estimations.
the samps and sampg methods, especially when using mgf as the plug-in π <dig> estimator, tend to have lower fdr on those outlier cases compared to q-value method.
the differences in power between q-value method and our proposed methods were very minor, compared to the difference between q-value and by method.
we also tested on the data-b simulations, and obtained the same result.
since we used the 3rd quartile of π <dig> estimated by samps and sampg, our estimations were biased but conservative.
with the same fdr control level our method would make smaller numbers of rejections than the q-value method, therefore the actual fdr and power of the genes selected were lower than that of the q-value method.
this was shown by the p-value in pair-wise comparison of actual fdr and power between the re-sampling based methods and q-value method in these simulations, where all fdr and most power comparisons were significant .
interestingly, comparing the mean of fdr and power, as shown in table  <dig> and  <dig>  we found that the samps and sampg methods, compared to the q-value method, can reduce the average fdr up to 40%, with a decrease of average power in most cases less than 1%.
in fact, with the highest correlation constant we tested, and in most cases for sampg.m method, the decrease of power was not even significant from the q-value method.
in contrast, the most conservative by method reduced the average fdr by more than 90% compared to the q-value method, but also reduced the average power by approximately 10% in cases with strong correlation.
the actual fdr in gene selections across the  <dig> simulations under each configuration by these methods was compared to that of q-value method using one-tail paired t-statistics.
the actual power in gene selections across the  <dig> simulations under each configuration by these methods was compared to that of q-value method using one-tail paired t-statistics.
a /fdrqvalue
b /fdrqvalue
c /fdrqvalue
d /fdrqvalue
e /fdrqvalue
a /powerqvalue
b /powerqvalue
c /powerqvalue
e /powerqvalue
f /powerqvalue
comparing the re-sampling strategy to bootstrap p-values directly
when there is strong correlation between genes, and thus also between p-values, bootstrapping p-values does not change the correlation structure and therefore the estimations are still unstable.
in contrast, re-sampling of samples or gene expression values as we proposed could address the variations induced by the correlation structure and therefore smooth the estimation.
for example, for each of the simulated datasets, we bootstrapped the p-values  <dig> times and then estimated π <dig> by the q-value method.
figure  <dig> shows the scatter plots of π <dig> estimated by the q-value method versus the 1st, 2nd and 3rd quartiles of the samps and sampg methods, as well as the estimations by bootstrapping p-values on the  <dig> data-m simulations with true π <dig> equals to  <dig> .
for the bootstrap method, as expected, the scatter plot showed that the median estimation of π <dig> was close to the estimation using original p-values for all cases.
whereas for the samps and sampg methods, since they smoothed the variation induced by correlation between p-values, the variation was much smaller than the q-value method.
especially for cases where the q-value method severely under- or over-estimated the π <dig>  the estimation by the samps and sampg methods were closer to the true value.
this was also observed in all simulated datasets that we have tested .
discussion
in actual microarray datasets, genes expression is often correlated due to co-regulation, sharing of transcription factor binding motifs, or technical reasons such as sequence similarity, cross-hybridization or signal leak during hybridization.
this is of critical importance for statistical inferences that rely on pooling of test statistics across genes  <cit> .
the distribution of p-values of these correlated genes can be viewed under a mixture model where groups of highly correlated genes share similar p-values and the whole distribution is actually a mixture of components corresponding to the groups of highly correlated genes.
the effect of strong and large scale correlations is equivalent to reducing the total number of independent components in this mixture model.
storey  <cit>  argued that subsets of genes fall into small but highly correlated groups due to co-regulation or cross-hybridization, but these groups are small in size and nearly independent with each other , therefore the uniformity of p-value distribution of true null genes would not be strongly affected.
however, other researchers, such as qiu et al  <cit> , have found that the impact of correlation may be quite strong.
it is also true that in our permutation study on a breast cancer microarray dataset, the distribution of p-values could be extremely far from uniform due to gene-gene correlation.
systems biology research has shown that biological gene networks have a scale free  <cit> , hierarchical structure  <cit> , where most of the genes are connected to a small number of other genes forming small groups of complexes, while some "hub" genes may be connected to large number of peripheral genes.
the distribution of connectivity degree  decreases with a power-law, which is much slower than the exponential decay expected in a random network  <cit> .
these gene-gene interactions may be reflected by co-regulation or correlation in expression under certain conditions, and the possible scale of gene interaction is unlimited given the scale free structure.
therefore, large scale correlation of gene expression levels is not surprising in microarray studies.
we have shown in our simulated study that with the growth of correlation structures, the p-value distribution of h <dig> genes increasingly deviates from a typical uniform distribution.
this may influence the estimation of π <dig> and the following statistical inferences.
the effect of strong correlation was also discussed by other researchers  <cit> .
to solve this problem, we proposed the samps and sampg methods.
these algorithms replaced the unbiased but unstable π <dig> estimation step in the q-value method with a model averaging procedure of re-sampling samples or furthermore re-sampling independently for each gene to partially break the correlations between genes.
strong correlation between genes will inevitably increase the variation of the π <dig> estimation, even though the variation could be partially smoothed by the re-sampling strategies proposed in this paper.
therefore, it is necessary to compromise between safety and efficiency; in this case, we would like to shrink the estimation toward  <dig> from an unbiased estimation to guarantee a strong control of fdr.
that is why we used the upper quartile instead of the median of the re-sampling estimations, although medians had a smaller mse in estimating π <dig> in our simulated studies.
we showed in our simulations that these plug-in revisions, compared to the q-value method, can greatly reduce the variation of the π <dig> estimation under strong gene-gene correlations, and enhance the performance of fdr control by reducing false discovery rate up to 40% with a reduction of power less than 1% compared to q-value method.
in our study, to create datasets with a known proportion of true null hypotheses while still having a similar correlation structure to that in actual microarray datasets, we developed  <dig> strategies to generate simulated datasets.
the first one, data-b, is simply a block-wise structure with arbitrary block size and intra-block correlation, but independent between blocks.
when the block size was small, this was similar to the "clumpy" hypothesis  <cit> .
when the block size became bigger, we showed that the correlations influenced the π <dig> estimation, and the re-sampling strategies we proposed improved the performance of gene selection by significantly reducing the fdr with a minor reduction of power.
to mimic a more realistic scenario, we also developed the data-m strategy to generate simulated data from actual microarray datasets by arbitrarily permuting a given proportion of the genes.
this permutation breaks the correlation between arbitrarily assigned differentially and non-differentially expressed genes, but maintains the correlation within the  <dig> groups of genes.
CONCLUSIONS
the samps and sampg methods we proposed and tested in this paper are simple revisions, but they greatly improved the π <dig> estimations.
the same approach using independent re-samples of expression values to estimate the π <dig> and then using the upper quartile of the re-sampling estimations in fdr control, could be applied to other fdr algorithms in data where strong correlation between hypotheses exists.
in this paper we considered the typical 2-sample comparison problem with a reasonable number of independent replicates.
for more complex problems, such as time-course studies, the sampg method  may not be able to be applied directly without a reasonable number of replicates in each time point, but the samps method  may be applicable if there is sound reason to assume the existence of stationary time patterns in the biological system under investigation.
availability and requirements
the r code of sampg and samps methods, and r code to generate simulated data sets data-b and data-m, is included as additional file  <dig>
requirements: r environment.
operating systems: windows xp or linux.
license: free.
authors' contributions
xl developed the basic strategies of sampg and samps methods, developed the methods to generate data-b and data-m simulation data, and did all simulation studies.
dpl helped xl on analyzing the experiments, comparing the methods, and improving the algorithms.
xl and dpl executed the writing together.
supplementary material
additional file 1
the file contains the r source code of functions to create data-b, data-m simulated data sets, and sampg, samps methods to estimate the π <dig>
click here for file
acknowledgements
the authors want to thank dr. ping ma, charles berry and anthony gamst for valuable advice and proofreading of the manuscript.
the comments by an anonymous referee and the assistant editor are also gratefully acknowledged.
the work was partially supported by nci grant ro <dig> ai <dig>  acs-irg #70- <dig>  and nci special cancer center support grant ca23100- <dig>

 microarray technology has become a standard experimental method in bio-medical research. in the analysis of microarray data, one of the most fundamental tasks is the identification of differentially expressed genes while controlling false positives and minimizing false negatives. this is a multiple hypothesis test problem which analyzes thousands or tens of thousands of genes simultaneously. in these tests we often need to control the false discovery among the rejected hypotheses under a pre-specified level while maintaining maximal power. thus, there is a trade off in the control of the type-i error between rejecting true null hypotheses versus accepting true alternative hypotheses. traditional bonferroni correction procedures are designed to control the family wise error rate, which guards against making one or more type i errors among a family of hypothesis tests. however, these procedures may be excessively conservative for microarray analysis where the number of hypotheses is very large and a substantial fraction of the genes are differentially expressed. a more appropriate approach is to control the false discovery rate, which is the proportion of type i errors among all rejected hypotheses. this approach is particularly useful in exploratory analyses, where the objective is to maximize the discovery of true positives, rather than guarding against one or more false positive results. a number of methods have been proposed to control the fdr given a population of hypothesis tests. these methods usually assume that the distribution of the test statistics, f, can be modeled by a mixture of two components: f = π0f + f     π = m0/m where f is the distribution of the test statistics under h which by definition equals to when using p-values when tests are independent, f is the distribution of the test statistics under h m is the number of true h m is the total number of hypotheses under consideration, and π is the proportion of true h the methods proposed by benjamini et al to control fdr do not estimate π0; therefore, they provide the strongest controls on fdr but have the lowest power compared to other methods that do so. in many actual applications where a considerable number of genes are differentially expressed, assuming π = may be too conservative causing loss of power. several alternative methods, such as nonparametric empirical bayesian pfdr criterion and its p-value equivalent called q-value method, bin-wise model, local fdr method, parametric beta-uniform mixture models, the lowest slope estimator, the spacing loess histogram method, the nonparametric mle method, the moment generating function approach, and the poisson regression approach, have all been proposed to estimate π by pooling test statistics and controlling fdr based on the estimated π in these methods, one of the critical steps is estimating the proportion of null hypotheses, π when using p-values, these estimations usually depend on the assumption that f follows a uniform distribution. this assumption, which is of critical importance for the methods of statistical inference that employ pooling test statistics across genes, is valid when all test hypotheses are independent and identically distributed. furthermore, when there are only weak correlations, or "clumpy" correlations, the uniform assumption is not strongly violated and the method remains adequate. however, in datasets with large scale strong correlations, the joint distribution of the test statistics will no longer be the product of marginal distributions, and the observed f will severely deviate from uniform, causing the current π estimation methods to become very unstable. increased variation and bias of π as well as fdr, was also observed by wu et al in datasets with strong local correlations. the effect of correlation on simultaneous significance tests was previously discussed theoretically, and a number of permutation based fdr control methods were proposed, such as sam, dchip, ge et al, meinshausen et al and efron. in these methods, the distribution of f was modeled empirically through permutations, which naturally considered the correlation. however, like benjamini et al, these methods don't estimate π0; therefore, in datasets with a large number of differentially expressed genes, the fdr control may be overly conservative with a loss of power. therefore we proposed re-sampling schemes, similar to model averaging in bagging methods, to reduce the variation in estimating π in datasets with strong correlation between gene expression values. our methods produced a more stable and conservative estimation of π and, therefore, provided stronger control of false discovery rate with only a minor sacrifice of power. implementation creating simulated data set to test the performance of various algorithms in estimating π we generated types of simulated datasets. both datasets had strong correlation between subsets of genes and a known proportion of true null hypotheses, to represent the log transformed microarray expression data. data-b the first simulation method is similar to the method used in qiu et al and wu et al. assume n samples and m genes, with n/ samples per class. the m genes were divided into k blocks with each block consisting of m/k genes. assume independence between blocks and constant correlation coefficient between genes within each block. for block l and sample j, we first created a block center vector blj = dl·xj + δlj, l =...,k, j =...,n where dl was the mean difference between groups, it equals to with probability π or was generated from beta distribution with parameter otherwise; xj was a group indicator; and δlj was i.i.d. n to represent sample specific noise. then the expression value of gene i in sample j in that block, ylji, was generated by ylji = ρ·blj + ·elji, i =...,m/k where ρ was the correlation constant which takes value between determining the correlation coefficient between genes within the block, and elji was i.i.d. n. the k blocks of genes were generated independently of each other and then pooled to form the whole dataset. we call this type of dataset which contains blocks of highly correlated genes data-b in our experiments. data-m the above data-b model is over-simplified in many aspects, and is still a "clumpy" structure, although the clumpiness can be pronounced. to mimic more realistic situations, we generated a second type of simulated data based on an actual human breast cancer microarray dataset obtained with affymetrix u plus microarrays. the dataset contains estrogen receptor positive cases and estrogen receptor negative cases. the data were normalized by the gcrma algorithm, and the gene expression levels were log2-transformed. according to published literatures, the er status is one of the most predominant partitioning factors for molecular classification of breast cancer. we therefore took some of the genes differentially expressed between the two classes as "truly h1" genes. we selected genes with differences in mean of log transformed expression levels between the two classes greater than. from these genes, each time we randomly chose to form our simulated dataset and then randomly picked π proportion of the genes to establish h genes by scrambling these genes together between samples. thus, we obtained a simulated dataset with a known number of null hypotheses while the correlations among both the differentially and non-differentially expressed genes were maintained. we call this type of dataset data-m in our experiments. estimating π by re-sampling strategy to get a better estimation of π in datasets with strong correlations, we proposed re-sampling methods to replace the original π estimation step in q-value method which estimates π directly from the p-value distribution. the first method, termed samps, re-samples without replacement 2/ of the samples from each class, calculates p-values for each gene, and then estimates the π from the p-values. for each dataset, we performed this procedure times and used the upper quartile to replace the π estimated by the q-value method. the second method, termed sampg, re-samples without replacement 2/ of the values from each class for each gene independently, calculates the p-values for each gene, and then estimates the π we also performed this procedure times for each dataset and used the upper quartile to replace the π estimated by the q-value method. after each re-sampling step, we have to feed the p-values into a π estimator. this estimator could be the π estimation by q-value method, the mgf method, or any other unbiased π estimators. in this paper, we tried to use both q-value and mgf as the plug-in estimator, and called them samps.q and sampg.q when using q-value as the plug-in estimator, or samps.m and sampg.m when using mgf as the plug-in estimator.  variation of the π estimation when genes are correlated to evaluate the impact of correlation structure in large datasets on the estimation of the proportion of true null hypotheses, we first evaluated current methods on a published microarray dataset. this dataset consists of estrogen receptor positive and estrogen receptor negative breast cancers. the gene expression levels were normalized by the gcrma algorithm and log transformed. expression values were filtered to eliminate low expressing genes with mean expression below and constant expressing genes with coefficient of variation below. a total of genes passed this filtering criterion. we bootstrapped datasets from this microarray data and used the q-value and twilight methods to estimate the proportion of null hypotheses. the π estimates were similar by the methods on each bootstrapped data set; however, both methods showed a large range of π among the bootstrapped datasets that varied from to. to further investigate the influence of gene correlation on the estimation of π we generated simulated dataset data-b which contained blocks of highly correlated genes, but all genes were non-differentially expressed. in our simulation we set the correlation constant ρ equal to, and created datasets containing or genes with or genes per block. after calculating the p-value for each gene, we calculated the coefficient of variation of the bar heights in the histogram of p-values by splitting the p-values into bins between with the width of each bin being. figure shows the histogram of p-values in one of the simulations having genes all independent with each other, and figure shows the histogram of p-values in another simulation having genes with genes per block. comparing figure and we can see that with highly correlated genes, the distribution of p-value deviated significantly from uniform, although none of the genes were differentially expressed. we created simulations for each type of data, calculated their cv of the bar height in the p-value histogram, and plotted the box plots of the cvs in figure in this simulated study, since all genes were non-differentially expressed, the p-values should follow a uniform distribution, and the histogram of p-values should be flat if genes were independent of each other. when the number of correlated genes in each block was small, for example, or genes per block, the distribution of p-values approximated a uniform distribution and the cv of the histogram of p-values were close to expected under the independent assumption. however, the cv of the histogram of p-values increased significantly with the growth of correlation structure. in other words, although none of the genes were differentially expressed, the distribution of p-value deviated increasingly from a uniform distribution when large groups of genes were correlated. we also calculated the cv of the histogram of p-values in our microarray dataset by randomly permuting the class labels, which makes all genes non-differentially expressed but still correlated, as well as randomly permuting all expression values across genes and samples which makes genes non-differentially expressed and also independent. the results of permutations showed a dramatically higher cv for the p-value histogram of datasets with only class labels permuted but with gene-gene correlations intact. comparing figure figure and figure it is apparent that increased correlation among genes greatly increased the deviation of p-value distribution from uniform. therefore, strong correlation structures will increase the variation in estimated π and inevitably subsequent statistical inferences and false discovery control procedures would be influenced by this unstable π estimation. improving the π estimation by re-sampling strategies to improve the estimation of π in datasets with strong correlations, we proposed two methods, termed samps and sampg, to replace the original π estimation step in the q-value method. in the samps algorithm, we used a model averaging strategy. we repeatedly sampled 2/ of the data from each class without replacement, calculated the p-values for genes, estimated π from the p-value distribution and finally used the upper quartile of the re-sampling estimations in the subsequent statistical inferences. in the sampg algorithm, we further broke down the correlations between genes and stabilized the π estimation. in this algorithm, we repeatedly sampled without replacement 2/ of the expression values from each class independently for each gene, calculated their p-values, estimate π and then used the upper quartile of the re-sampling estimations in subsequent analysis. for the choice of the plug-in π estimator, we tried to use both q-value and moment generating function methods, and called them sampg.q, samps.q and sampg.m, samps.m, respectively. to test the variation in π estimation induced by strong correlation structures, and the performance of our proposed samps and sampg methods, we created simulated datasets data-b and data-m, with true π being,,,, and the correlation constant ρ being, and, respectively. we created datasets for each combination of these control parameters. in data-b, we created genes forming blocks with genes per block to simulate large scale correlation between genes. in data-m we randomly selected genes with differences in the mean of log transformed expression levels greater than, then randomly scrambled 90%, 80%, 70% or 60% of them to create datasets with known proportion of true null hypothesis and strong gene-gene correlations. we applied the samps and sampg strategies to estimate π for each of the simulated datasets, and compared our results to a number of other methods. the π estimation methods we used and their corresponding r functions are: the lowest slope estimator with parameter determined via bootstrap; function fdr.estimate.eta0 q-value method with tuning parameter chosen by smoother method; function qvalue the spacing loess histogram; function splosh the beta-uniform mixture model; function bum.mle the nonparametric mle method; function convest the poisson regression approach; function p mom the moment generating function; function p mom the lowest slope estimator with parameter determined adaptively; function fdr.estimate.eta0 sampg method, with q-value plug in, 2nd quartile sampg method, with q-value plug in, 3rd quartile samps method, with q-value plug in, 2nd quartile samps method, with q-value plug in, 3rd quartile sampg method, with mgf plug in, 2nd quartile sampg method, with mgf plug in, 3rd quartile samps method, with mgf plug in, 2nd quartile samps method, with mgf plug in, 3rd quartile the π estimations were shown as boxplots in figure and the mean square error were listed in table we then selected genes with fdr controlled under level based on the π estimated by our method, calculated the actual false discovery rate and power, and compared our results to that of fdr controlling method proposed by benjamini et al. with correlations considered, the permutation based fdr controlling method, and q-value method. the boxplot of actual false discovery rate and power on the types of data-m simulations were shown in figure from figure it can be seen that both the by and howmany methods provided strong control of fdr, their actual fdr level was much lower than, and their power to detect true alternatives were much lower than methods where the proportion of true null hypotheses was estimated and used in the fdr control. comparing the q-value method and our proposed samps and sampg method, for majority of the cases the actual fdr level were still controlled below level, although some outliers exist with actual fdr up to ~, due to the unstable pi_ estimations. the samps and sampg methods, especially when using mgf as the plug-in π estimator, tend to have lower fdr on those outlier cases compared to q-value method. the differences in power between q-value method and our proposed methods were very minor, compared to the difference between q-value and by method. we also tested on the data-b simulations, and obtained the same result. since we used the 3rd quartile of π estimated by samps and sampg, our estimations were biased but conservative. with the same fdr control level our method would make smaller numbers of rejections than the q-value method, therefore the actual fdr and power of the genes selected were lower than that of the q-value method. this was shown by the p-value in pair-wise comparison of actual fdr and power between the re-sampling based methods and q-value method in these simulations, where all fdr and most power comparisons were significant. interestingly, comparing the mean of fdr and power, as shown in table and we found that the samps and sampg methods, compared to the q-value method, can reduce the average fdr up to 40%, with a decrease of average power in most cases less than 1%. in fact, with the highest correlation constant we tested, and in most cases for sampg.m method, the decrease of power was not even significant from the q-value method. in contrast, the most conservative by method reduced the average fdr by more than 90% compared to the q-value method, but also reduced the average power by approximately 10% in cases with strong correlation. the actual fdr in gene selections across the simulations under each configuration by these methods was compared to that of q-value method using one-tail paired t-statistics. the actual power in gene selections across the simulations under each configuration by these methods was compared to that of q-value method using one-tail paired t-statistics. a /fdrqvalue b /fdrqvalue c /fdrqvalue d /fdrqvalue e /fdrqvalue a /powerqvalue b /powerqvalue c /powerqvalue e /powerqvalue f /powerqvalue comparing the re-sampling strategy to bootstrap p-values directly when there is strong correlation between genes, and thus also between p-values, bootstrapping p-values does not change the correlation structure and therefore the estimations are still unstable. in contrast, re-sampling of samples or gene expression values as we proposed could address the variations induced by the correlation structure and therefore smooth the estimation. for example, for each of the simulated datasets, we bootstrapped the p-values times and then estimated π by the q-value method. figure shows the scatter plots of π estimated by the q-value method versus the 1st, 2nd and 3rd quartiles of the samps and sampg methods, as well as the estimations by bootstrapping p-values on the data-m simulations with true π equals to. for the bootstrap method, as expected, the scatter plot showed that the median estimation of π was close to the estimation using original p-values for all cases. whereas for the samps and sampg methods, since they smoothed the variation induced by correlation between p-values, the variation was much smaller than the q-value method. especially for cases where the q-value method severely under- or over-estimated the π the estimation by the samps and sampg methods were closer to the true value. this was also observed in all simulated datasets that we have tested. discussion in actual microarray datasets, genes expression is often correlated due to co-regulation, sharing of transcription factor binding motifs, or technical reasons such as sequence similarity, cross-hybridization or signal leak during hybridization. this is of critical importance for statistical inferences that rely on pooling of test statistics across genes. the distribution of p-values of these correlated genes can be viewed under a mixture model where groups of highly correlated genes share similar p-values and the whole distribution is actually a mixture of components corresponding to the groups of highly correlated genes. the effect of strong and large scale correlations is equivalent to reducing the total number of independent components in this mixture model. storey argued that subsets of genes fall into small but highly correlated groups due to co-regulation or cross-hybridization, but these groups are small in size and nearly independent with each other, therefore the uniformity of p-value distribution of true null genes would not be strongly affected. however, other researchers, such as qiu et al, have found that the impact of correlation may be quite strong. it is also true that in our permutation study on a breast cancer microarray dataset, the distribution of p-values could be extremely far from uniform due to gene-gene correlation. systems biology research has shown that biological gene networks have a scale free, hierarchical structure, where most of the genes are connected to a small number of other genes forming small groups of complexes, while some "hub" genes may be connected to large number of peripheral genes. the distribution of connectivity degree decreases with a power-law, which is much slower than the exponential decay expected in a random network. these gene-gene interactions may be reflected by co-regulation or correlation in expression under certain conditions, and the possible scale of gene interaction is unlimited given the scale free structure. therefore, large scale correlation of gene expression levels is not surprising in microarray studies. we have shown in our simulated study that with the growth of correlation structures, the p-value distribution of h genes increasingly deviates from a typical uniform distribution. this may influence the estimation of π and the following statistical inferences. the effect of strong correlation was also discussed by other researchers. to solve this problem, we proposed the samps and sampg methods. these algorithms replaced the unbiased but unstable π estimation step in the q-value method with a model averaging procedure of re-sampling samples or furthermore re-sampling independently for each gene to partially break the correlations between genes. strong correlation between genes will inevitably increase the variation of the π estimation, even though the variation could be partially smoothed by the re-sampling strategies proposed in this paper. therefore, it is necessary to compromise between safety and efficiency; in this case, we would like to shrink the estimation toward from an unbiased estimation to guarantee a strong control of fdr. that is why we used the upper quartile instead of the median of the re-sampling estimations, although medians had a smaller mse in estimating π in our simulated studies. we showed in our simulations that these plug-in revisions, compared to the q-value method, can greatly reduce the variation of the π estimation under strong gene-gene correlations, and enhance the performance of fdr control by reducing false discovery rate up to 40% with a reduction of power less than 1% compared to q-value method. in our study, to create datasets with a known proportion of true null hypotheses while still having a similar correlation structure to that in actual microarray datasets, we developed strategies to generate simulated datasets. the first one, data-b, is simply a block-wise structure with arbitrary block size and intra-block correlation, but independent between blocks. when the block size was small, this was similar to the "clumpy" hypothesis. when the block size became bigger, we showed that the correlations influenced the π estimation, and the re-sampling strategies we proposed improved the performance of gene selection by significantly reducing the fdr with a minor reduction of power. to mimic a more realistic scenario, we also developed the data-m strategy to generate simulated data from actual microarray datasets by arbitrarily permuting a given proportion of the genes. this permutation breaks the correlation between arbitrarily assigned differentially and non-differentially expressed genes, but maintains the correlation within the groups of genes.  the samps and sampg methods we proposed and tested in this paper are simple revisions, but they greatly improved the π estimations. the same approach using independent re-samples of expression values to estimate the π and then using the upper quartile of the re-sampling estimations in fdr control, could be applied to other fdr algorithms in data where strong correlation between hypotheses exists. in this paper we considered the typical 2-sample comparison problem with a reasonable number of independent replicates. for more complex problems, such as time-course studies, the sampg method may not be able to be applied directly without a reasonable number of replicates in each time point, but the samps method may be applicable if there is sound reason to assume the existence of stationary time patterns in the biological system under investigation. availability and requirements the r code of sampg and samps methods, and r code to generate simulated data sets data-b and data-m, is included as additional file requirements: r environment. operating systems: windows xp or linux. license: free. authors' contributions xl developed the basic strategies of sampg and samps methods, developed the methods to generate data-b and data-m simulation data, and did all simulation studies. dpl helped xl on analyzing the experiments, comparing the methods, and improving the algorithms. xl and dpl executed the writing together. supplementary material additional file 1 the file contains the r source code of functions to create data-b, data-m simulated data sets, and sampg, samps methods to estimate the π click here for file acknowledgements the authors want to thank dr. ping ma, charles berry and anthony gamst for valuable advice and proofreading of the manuscript. the comments by an anonymous referee and the assistant editor are also gratefully acknowledged. the work was partially supported by nci grant ro ai acs-irg #70- and nci special cancer center support grant ca23100- 
when conducting multiple hypothesis tests, it is important to control the number of false positives, or the false discovery rate. therefore, we developed a re-sampling strategy to reduce the variation by breaking the correlations between gene expression values, then using a conservative strategy of selecting the upper quartile of the re-sampling estimations to obtain a strong control of fdr. with simulation studies and perturbations on actual microarray datasets, our method, compared to competing methods such as q-value, generated slightly biased estimates on the proportion of null hypotheses but with lower mean square errors. when selecting genes with controlling the same fdr level, our methods have on average a significantly lower false discovery rate in exchange for a minor reduction in the power. our objective was to develop methods to control the fdr while maintaining a greater level of power in highly correlated datasets by improving the estimation of the proportion of null hypotheses. we showed that when strong correlation exists among the data, which is common in microarray datasets, the estimation of the proportion of null hypotheses could be highly variable resulting in a high level of variation in the fdr. these methods usually depend on the assumption that the test statistics are independent. several methods have been proposed, such as the q-value method, to estimate the proportion of true null hypothesis among the tested hypotheses, and use this estimation in the control of fdr.
 when conducting multiple hypothesis tests, it is important to control the number of false positives, or the false discovery rate. however, there is a tradeoff between controlling fdr and maximizing power. several methods have been proposed, such as the q-value method, to estimate the proportion of true null hypothesis among the tested hypotheses, and use this estimation in the control of fdr. these methods usually depend on the assumption that the test statistics are independent. however, many types of data, for example microarray data, often contain large scale correlation structures. our objective was to develop methods to control the fdr while maintaining a greater level of power in highly correlated datasets by improving the estimation of the proportion of null hypotheses.  we showed that when strong correlation exists among the data, which is common in microarray datasets, the estimation of the proportion of null hypotheses could be highly variable resulting in a high level of variation in the fdr. therefore, we developed a re-sampling strategy to reduce the variation by breaking the correlations between gene expression values, then using a conservative strategy of selecting the upper quartile of the re-sampling estimations to obtain a strong control of fdr.  with simulation studies and perturbations on actual microarray datasets, our method, compared to competing methods such as q-value, generated slightly biased estimates on the proportion of null hypotheses but with lower mean square errors. when selecting genes with controlling the same fdr level, our methods have on average a significantly lower false discovery rate in exchange for a minor reduction in the power. 