Ranking analysis of F-statistics for microarray data
© Tan et al; licensee BioMed Central Ltd. 2008
Received: 23 November 2007
Accepted: 06 March 2008
Published: 06 March 2008
Microarray technology provides an efficient means for globally exploring physiological processes governed by the coordinated expression of multiple genes. However, identification of genes differentially expressed in microarray experiments is challenging because of their potentially high type I error rate. Methods for large-scale statistical analyses have been developed but most of them are applicable to two-sample or two-condition data.
We developed a large-scale multiple-group F-test based method, named ranking analysis of F-statistics (RAF), which is an extension of ranking analysis of microarray data (RAM) for two-sample t-test. In this method, we proposed a novel random splitting approach to generate the null distribution instead of using permutation, which may not be appropriate for microarray data. We also implemented a two-simulation strategy to estimate the false discovery rate. Simulation results suggested that it has higher efficiency in finding differentially expressed genes among multiple classes at a lower false discovery rate than some commonly used methods. By applying our method to the experimental data, we found 107 genes having significantly differential expressions among 4 treatments at <0.7% FDR, of which 31 belong to the expressed sequence tags (ESTs), 76 are unique genes who have known functions in the brain or central nervous system and belong to six major functional groups.
Our method is suitable to identify differentially expressed genes among multiple groups, in particular, when sample size is small.
Microarray gene expression technology, which profiles the expression of multiple genes in parallel [1, 2], affords the means for globally exploring physiological and pathological processes  regulated by the coordinated expression of thousands of genes . However, identification of genes that are differentially expressed in large-scale gene expression experiments requires global statistical methods rather than traditional statistical methods based on single hypothesis testing. A variety of multiple-testing procedures, such as the Bonferroni procedure, Holm procedure , Hochberg procedure , Benjamini-Hochberg (BH) procedure , and Benjamini-Liu (BL) procedures  have already been developed for testing a large family of null hypotheses. The first three methods bound the family-wise-error rate (FWER) that is the probability of at least one false positive over all tests and hence remain too stringent and have lower power for finding genes from the real data sets. The last two methods have an upper bound for the false discovery rate (FDR) with both strong and weak controls  and require a large sample size for valid p-values. Tusher et al.  has proposed a ranking statistic approach based on permutation for resampling. However, permutation is not a desirable approach to estimating null distribution in microarray data [10–12] because in general a microarray dataset has a large number of genes but small sample sizes  due to cost. Permutation fails to remove treatment effect and due to small sample sizes the difference of treatment effects between permutated groups may become a main component in differences between group means so that the estimated null distribution is not well approximate to the true null distribution ( and also see Appendix in Tan et al. ). For example, Xie et al.  found that the estimated null F-distribution based on permutation has a larger variance and a heavier tail compared to the true null F-distribution, which leads to a potential loss of power. Similar phenomenon was also observed in comparison of the estimated null t-distribution to the true null t-distribution . To remove the group or treatment effects on the estimated null distribution, Tan et al.  developed a random splitting (RS) approach. Since treatment effects are completely eliminated, the estimated null distribution obtained by the RS method is smooth, stable and approximate true null distribution well.
For the multi-class microarray data, the analysis of variance (ANOVA) is useful to identify differentially expressed genes . In ANOVA, the F-test is a generalization of the t-test that allows for comparison of more than two samples . However, due to small sample sizes, the classical F-test is also subject to the same problems as the t-test: bias and unstable estimates of gene-specific variances. To tackle this issue, many authors [15–19] proposed modified F-statistics. However, like the classical F-test, these modified F-tests still suffer from high false-positive rates because (i) the sample size is often so limited that the asymptotic F distribution is not accurate enough to obtain a valid p-value and (ii) they appeal to multiple-testing procedures such as the Bonferroni procedure or the BH-procedure. As mentioned above, these multiple-testing procedures have a basic requirement that sample sizes are large enough for valid p-values. In microarray data, the requirement is not realistic. Based on consideration of these problems, we here propose a novel statistical method for the analysis of multi-class gene-expression data called Ranking Analysis of F-statistics (RAF). RAF is a natural extension of our previous work, i.e., the ranking analysis of microarrary (RAM) for two-class t-tests . It works on finding genes that are differentially expressed among multiple treatment groups by comparing the ordered real F-statistics with the ordered estimated null F-statistics and implementing a two-simulation strategy to estimate the false discovery rate (FDR).
Animal model and design
Studies were performed on male stroke-resistant SHR/N (CRiv) (SHRSR) and stroke-prone SHR/A3 (Heid) (SHRSP) rats from a breeding colony maintained by the investigators as previously described . Age-matched male rats from each strain (N = 12 SHRSP and 12 SHRSR) were fed a standard rat chow and water ad libitum until age 8 weeks. Subsequently, animals from each strain were randomized to one of 2 dietary regimens (N = 6 in each strain × diet group): a "stroke-permissive diet" high in sodium (HS) (0.63% potassium, 0.37% sodium) and 1% NaCl drinking solution; a "stroke-protective diet" low in sodium and high in potassium (LS) (1.3% potassium, 0.35% sodium) and regular drinking water. All animals were housed at 23°C on a 12-hour light-dark cycle. Rats were sacrificed at 12 weeks of age, and brain tissue was collected for RNA extraction and subsequent microarray analysis. The study protocols were approved by the Animal Care Committee of the University of Texas – Houston. Since strain and dietary factor each have only two levels, we here treated them as one-way in statistical analysis instead of two-way, that is, we are neither interested in strain effects alone nor in dietary effects alone but focus on their combined contributions to gene expression. Thus, HS-SHRSPs, LS-SHRSPs, HS-SHRSRs, and LS-SHRSRs are viewed as four treatment groups for the purpose of the analyses.
Microarray analysis was performed as described by Lockhart et al. . Briefly, 10 μ g total RNA extracted from each of the 24 rats was used to synthesize cDNA, which was then used as a template to generate biotinylated cRNA. cRNA was fragmented and hybridized to a Test2 chip to verify quality and quantity of the samples. Each sample was then hybridized to a RGU34A array (Affymetrix, Santa Clara, CA). After hybridization, each array was washed and scanned, and fluorescence values were measured and normalized using the Affymetrix Microarray Suite v. 5.0 software.
where σ2 (G g ) and σ2 (e g ) are inter- and intra-group variances of the expression values of gene g, respectively. In the conventional F-tests, for example, significance of p = 0.01 in the context of the standard F distribution is for a single hypothesis to be tested; therefore, it is unsuitable to microarray data because in a microarray experiment for 10,000 genes we would expect to identify 100 genes by chance . To address this problem, an alternative approach is to rank genes by magnitude of their F values so that F1 is the largest value, F2 is the second largest value, etc., and Fg*is the g*th largest value where g* is a rank position of gene g. Thus, we have a ranking F-test where
Fg*- fg*> Δ (2)
indicates that the expression variation of gene g among multiple groups (or multiple conditions) is significant. In Eq. (2), fg*is the expectation of Fg*under the null hypothesis and Δ is a given threshold.
Estimation of fg*
Therefore, the null hypothesis is equivalent to F g = f g because σ2 (τ g ) = 0 under the null hypothesis. Note that σ2 (τ g ) = 0 means the treatment effects τ gi = ... = τ gn = μ g . In order to do a ranking F-test, it is necessary to obtain a good estimate of fg*. In the two-group scenario, Tusher et al.  employed a permutation approach to estimate the expected t-statistics. The permutation process cannot completely clear the treatment effect in the ranked d-statistics so that the estimated ranked d-statistics distribution is biased against its null distribution and unstable, in particular, when sample sizes are small (see Appendix A in Tan et al., ). Tan et al.  developed a "Randomly Splitting" (RS) approach to estimate the null distribution of t-statistics. In this study, we extended the RS approach to estimating the null distribution of F-statistics.
Note that since treatment effect is completely removed from the difference between two sub-sample means, the difference is pure noise. We rank across all g and let denote the value in ordered position g* at split J. After running M splits, we have M values of for position g*. Thus fg*in Eq. (2) can be estimated by the average of over all M splits, i.e., .
Estimation of FDR
To identify genes whose expression is significantly changed among multiple conditions, it is necessary to estimate the FDR for a given threshold [7, 22]. Here we propose a two-simulation approach for FDR estimation . Consider a series of threshold values Δ k (k = 1,...L) and let N k be the number of genes that are claimed as significant by RAF at threshold Δ k . N k comprises two parts: the number N k (t) of the true positives and the number N k (f) of the false positives, i.e., N k = N k (t) + N k (f). Thus, given a threshold Δ k , FDR is defined as λ k = N k (f)/N k . N k (f) is unknown, hence λ k must be estimated. Many approaches such as BH procedure [7, 22], BL procedure , Storey's procedure [23, 24], and Pounds and Cheng's procedure  have been proposed to estimate the FDR. These approaches, however, are based on the assumption that the tests are independent. As mentioned previously, this assumption may not be met in practice. Therefore, these methods may not be suitable to our ranking test. Based on the fact that sampling distribution fluctuates around the expected distribution via permutation, Tusher et al.  developed a permutation-based estimator to estimate FDR in the ranking tests. It has been proved, however, in theory and in simulation that when the sample sizes are small, the number of permutations is very limited so that the treatment effects cannot be removed in the permutated data . As a result, the estimator is biased for a given threshold. Here we extend the interval approach by Tan et al.  to the ranking analysis of F-statistics. In this approach, we first construct an estimated interval of the true FDR, and then we find a reasonable estimate of FDR. This interval is based on the complete and partial null distributions given by two simulations.
In simulation 1, for each gene, n samples (groups) each having r replicates are generated from normal distributions with a set of sample means () and a set of sample error variances [s2 (eg 1),..., s2 (e gn )]. Here we set and i is a randomly chosen group from the observed data, for each of a half of the genes with the null effect that the group variance is zero, i.e., σ2 (G g ) = 0 and for each of the other half with unknown effect that the group variance is larger than or equal to zero, i.e., σ2 (G g ) ≥ 0. s2 (e gi ) is set to be equal to σ2 (e gi ) where and σ2 (e gi ) are the observed values from the real microarray data set.
The second simulation for estimating FDR is carried out in the following fashion. n samples (groups) each having r replicates for each gene are generated from normal distributions with a set of sample means, and a set of sample variances s2 (eg 1) = σ2 (eg 1),..., s2 (e gn ) = σ2 (e gn ).
as the second function of threshold (see Figure 1). In particular, we let λ2k= 1 if N k = N2k= 0 because λ2k= 1 when N k = 0 and N2k> 0.
Thus, an interval for FDR estimation at threshold Δ k can be constructed between λ1kand λ2k. The third function of threshold for FDR estimation is given as
λ3k= α k λ1k+ β k λ2k
Note that as shown in the simulation result section, λ2kis an underestimate of λ k and λ1kis an overestimate of FDR when the threshold Δ k < Δ*. However, the situation is reversed when threshold Δ k > Δ*. This is because N1kbecomes very small when Δ k > Δ* so that λ1kbecomes very small whereas, from Eq. (10), λ2kslowly decreases if N k > N2kor increases if N k <N2kas threshold increases. In addition, when the microarray data have no treatment effects for all the genes detected, then λ1k= λ2k= λ3k= 1, leading to = 1
where p k = (N k - Nk+1)/(1 + N k - Nk+1) and q k = 1 - p k . Eq. (13) suggests that λk+1= λ k if N k = Nk+1. The number of false discoveries among those found to be significant at threshold Δ k in the observed data is estimated by . Figure 1 shows that the curve of agrees well with that of λ k .
Estimation of the null distribution of F-statistics
To examine if the empirical distributions obtained by the RS approach are appropriate for the analysis of the expression data, we simulated a microarray data set consisting of 3770 genes and four groups each having 6 replicates using one group mean and error variance for each gene. Thus, the simulation without treatment effect generated a set of pure noise data.
Estimation of FDR
Since it is in general unknown if a given gene expresses differently among multiple conditions, it is not necessarily best to use real data of gene expression to evaluate an FDR estimator. But simulation is a useful approach to doing such a task. Therefore, we also conducted a computer simulation for comparing expression status (significant or not) of a gene identified by a method with its real status. This simulation was also based on our real data set of 3770 genes. Treatment effect τ on expression variation was set for 30 % of the genes and assigned in 4 groups. The mean expression value of gene g was set and for the 4 groups where τ = 100U, 0 <U ≤ 1, is overall observed average for gene g, and each group has 6 replicates. Obviously, treatment effect τ on expression changes randomly with genes in our simulation, which would make it more difficult to identify differentially expressed genes than the simulations with a fixed treatment effect. Figure 1 displays a comparison between RAF estimated and true FDRs. One can see that the RAF estimate curve is very close to the true FDR curve given a series of thresholds.
Efficiencies of different methods in finding genes differentially expressed among multiple groups
Efficiencies of different methods in identifying genes differentially expressed among four groups each with 6 replicates in 30 simulated datasets
Difference between ENFP and TNFP
C(d ≥ 0)
λ = 0.05
λ = 0.05
0.04 <λ ≤ 0.05
0.03 <λ ≤ 0.04
0.02 <λ ≤ 0.03
0.01 <λ ≤ 0.02
0.00 <λ ≤ 0.01
λ = 0.00
0.04 <λ ≤ 0.05
0.03 <λ ≤ 0.04
0.02 <λ ≤ 0.03
0.01 <λ ≤ 0.02
0.00 <λ ≤ 0.01
λ = 0.00
Comparison between SAM and RAF in finding genes differentially expressed among four classes in a simulated data set of small sample size (n = 4)
Number of significances
Number of false positives
Number of significances
Number of false positive
Array findings by RAF
The results of RAF identifying genes differentially expressed among HS-SHRSPs, LS-SHRSPs, HS-SHRSRs, and LS-SHRSRs.
Number of genes called significant
Number of false discoveries
Independent verification of array findings
Fornage et al (2003) used TagMan assay to measure the relative expressions of 7 genes encoding atrial natriuretic peptide (Anp), the neurotrophin receptor protein tyrosine kinase (TrkB, short), casein kinase 2 (Ck2), complexin 2 (Cplx2), stearoyl CoA desaturase 2 (Scd2), glycerol-3-phosophate acyltransterase (Gpan), and inositol 1,4,5-triphosphate receptor (Itpr1). They found these 7 genes had significantly differentially expressed between SHRSP and SHR strains with p < 0.05. Except that genes Anp and Gpan were out of our data, genes for TrkB (short), Cplx2, and Scd2 called significant differential expressions at FDR<0.7%, and for CK2 and Itpr1 at FDR = 0.7% were found among HS-SHRSP, LS-SHRSP, HS-SHR, and LS-SHR strains. Interestingly, Tropea et al  also found the genes encoding glutamate receptor (GluR-A) and GABA receptor had significant expression difference between two groups of mice treated by dark rearing and monocular deprivation.
To our knowledge, the ranking analysis of F-statistics for identifying differentially expressed genes among multiple groups (classes) has not been reported. There are two main difficulties to be overcome: (a) estimate of the null F-distribution and (b) estimate of FDR. In conventional statistical methods, permutation is very popular to generate empirical distributions as estimates of the null distributions. However, the permutation approach may not be suitable for microarray data [10–13] because in general microarray experiments have a small sample size due to cost, as a result, treatment effect residues that cannot be removed are amplified in permutation distribution and resulting estimated null distribution has a heavier tail compared to true null distribution . This would results in two consequences: (a) the estimated null distribution is not stable, which, as seen in Table 3, leads to low conservativeness of estimate of FDR, and (b) low power. Our RAF method is successful because the -distribution obtained by applying the RS approach  does not contain treatment effects and hence is a desirable estimate of the null F-distribution, which is supported by the fact that the observed and simulated results agree well with those expected by theory. Therefore, the number (M) of splits is much smaller than that of permutations for estimate of the null F-distribution. Simulation results showed that 50 splits are enough to obtain a stable and smooth -distribution. In addition, since the -distribution is generated from all the genes detected on microarrays and does not contain treatment effect residences, impact of sample size on the -distribution is very weak. However, we also noted that the -distribution would underestimate the null F-distribution when sample sizes are smaller than 4. In this situation, Eq. (7) should be changed to .
FDR is often used to control the error rate in the BH procedure , the BL procedure , and in SAM . In practice, for a ranking test, it is necessary to obtain an accurate estimate of FDR. In SAM, FDR is estimated through the permutation approach in which fluctuations around expectation occur among permutated samples. The fluctuations would be dependent on the data itself, i.e., sample size, treatment effect, and data noise. In addition, as indicated above, permutation fails to remove the treatment effects in the data permuted from the microarray data with a small sample size so that the fluctuations are not purely due to random errors. Thus, this approach may give a biased estimate of FDR for a given threshold. The RAF estimator is based on a two-simulation strategy and hence avoids these problems of the SAM estimator, that is, its accuracy is not affected by sample size, treatment effect, and noise. As a result, the number (B) of simulations is also relatively small. Our simulation study indicates that more than 40 simulated data sets (B ≥ 40) would produce stable estimates of FDR across all given thresholds.
Our current RAF method can be readily extended to other test statistics such as Brown-Forsythe test statistic , Welch test statistic , and Cochran test statistic  by replacing F-statistic with the respective statistics.
We developed a new statistical method that is suitable for analyzing microarray data to identify differentially expressed genes among multiple groups, especially, when sample size is small.
This research was supported by grants from the U.S. National Institutes of Health (HL69126) to MF.
- DeRisi J, Penland L, Brown PO, Bittner ML, Meltzer PS, Ray M, Chen Y, Su YA, Trent JM: Use of a cDNA microarray to analyse gene expression patterns in human cancer. Nature genetics 1996, 14(4):457–460. 10.1038/ng1296-457View ArticlePubMedGoogle Scholar
- DeRisi JL, Iyer VR, Brown PO: Exploring the metabolic and genetic control of gene expression on a genomic scale. Science 1997, 278(5338):680–686. 10.1126/science.278.5338.680View ArticlePubMedGoogle Scholar
- Kim YD, Sohn NW, Kang C, Soh Y: DNA array reveals altered gene expression in response to focal cerebral ischemia. Brain research bulletin 2002, 58(5):491–498. 10.1016/S0361-9230(02)00823-7View ArticlePubMedGoogle Scholar
- Kerr MK, Martin M, Churchill GA: Analysis of variance for gene expression microarray data. J Comput Biol 2000, 7(6):819–837. 10.1089/10665270050514954View ArticlePubMedGoogle Scholar
- Holm S: A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics 1979, 6: 65–70.Google Scholar
- Hochberg Y: A sharper Bonferroni procedure for multiple tests of significance. Biometrika 1988, 75(4):800–802. 10.1093/biomet/75.4.800View ArticleGoogle Scholar
- Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B Methodological 1995, 57(1):289–300.Google Scholar
- Benjamini Y, Liu W: A step-down multiple hypotheses tesing procedure that controls the false discovery rate under independence. Journal of Statistical Planning and Inference 1999, 82: 163–170. 10.1016/S0378-3758(99)00040-3View ArticleGoogle Scholar
- Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proceedings of the National Academy of Sciences of the United States of America 2001, 98(9):5116–5121. 10.1073/pnas.091062498PubMed CentralView ArticlePubMedGoogle Scholar
- Zhao Y, Pan W: Modified nonparametric approaches to detecting differentially expressed genes in replicated microarray experiments. Bioinformatics (Oxford, England) 2003, 19(9):1046–1054. 10.1093/bioinformatics/btf879View ArticleGoogle Scholar
- Pan W: On the use of permutation in and the performance of a class of nonparametric methods to detect differential gene expression. Bioinformatics (Oxford, England) 2003, 19(11):1333–1340. 10.1093/bioinformatics/btg167View ArticleGoogle Scholar
- Xie Y, Pan W, Khodursky AB: A note on using permutation-based false discovery rate estimates to compare different analysis methods for microarray data. Bioinformatics 2005, 21(23):4280–4288. 10.1093/bioinformatics/bti685View ArticlePubMedGoogle Scholar
- Gao X: Construction of null statistics in permutation-based multiple testing for multi-factorial microarray experiments. Bioinformatics (Oxford, England) 2006, 22(12):1486–1494. 10.1093/bioinformatics/btl109View ArticleGoogle Scholar
- Tan YD, Fornage M, Fu YX: Ranking analysis of microarray data: a powerful method for identifying differentially expressed genes. Genomics 2006, 88(6):846–854. 10.1016/j.ygeno.2006.08.003PubMed CentralView ArticlePubMedGoogle Scholar
- Cui X, Churchill GA: Statistical tests for differential expression in cDNA microarray experiments. Genome biology 2003, 4(4):210. 10.1186/gb-2003-4-4-210PubMed CentralView ArticlePubMedGoogle Scholar
- Li H, Wood CL, Liu Y, Getchell TV, Getchell ML, Stromberg AJ: Identification of gene expression patterns using planned linear contrasts. BMC Bioinformatics 2006, 7: 245. 10.1186/1471-2105-7-245PubMed CentralView ArticlePubMedGoogle Scholar
- Chen D, Liu Z, Ma X, Hua D: Selecting genes by test statistics. Journal of biomedicine & biotechnology 2005, 2005(2):132–138. 10.1155/JBB.2005.132View ArticleGoogle Scholar
- Tsai PW, Lee ML: Split-plot microarray experiments: issues of design, power and sample size. Applied bioinformatics 2005, 4(3):187–194.View ArticlePubMedGoogle Scholar
- Cui X, Hwang JT, Qiu J, Blades NJ, Churchill GA: Improved statistical tests for differential gene expression by shrinking variance components estimates. Biostatistics (Oxford, England) 2005, 6(1):59–75.View ArticleGoogle Scholar
- Fornage M, Swank MW, Boerwinkle E, Doris PA: Gene expression profiling and functional proteomic analysis reveal perturbed kinase-mediated signaling in genetic stroke susceptibility. Physiological genomics 2003, 15(1):75–83.View ArticlePubMedGoogle Scholar
- Lockhart DJ, Dong H, Byrne MC, Follettie MT, Gallo MV, Chee MS, Mittmann M, Wang C, Kobayashi M, Horton H, et al.: Expression monitoring by hybridization to high-density oligonucleotide arrays. Nature biotechnology 1996, 14(13):1675–1680. 10.1038/nbt1296-1675View ArticlePubMedGoogle Scholar
- Benjamini Y, Yekutieli D: The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 2001, 29(4):1165–1188. 10.1214/aos/1013699998View ArticleGoogle Scholar
- Storey JD: A direct approach to false discovery rates. Journal of the Royal Statistical Society Series B Methodological 2002, 64(3):479–498. 10.1111/1467-9868.00346View ArticleGoogle Scholar
- Storey JD, Taylor JE, Siegmund D: Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach. Journal of the Royal Statistical Society Series B Methodological 2004, 66(1):187–205. 10.1111/j.1467-9868.2004.00439.xView ArticleGoogle Scholar
- Pounds S, Cheng C: Robust estimation of the false discovery rate. Bioinformatics (Oxford, England) 2006, 22(16):1979–1987. 10.1093/bioinformatics/btl328View ArticleGoogle Scholar
- Cho H, Lee JK: Bayesian hierarchical error model for analysis of gene expression data. Bioinformatics (Oxford, England) 2004, 20(13):2016–2025. 10.1093/bioinformatics/bth192View ArticleGoogle Scholar
- Tropea D, Kreiman G, Lyckman A, Mukherjee S, Yu H, Horng S, Sur M: Gene expression changes and molecular pathways mediating activity-dependent plasticity in visual cortex. Nature neuroscience 2006, 9(5):660–668. 10.1038/nn1689View ArticlePubMedGoogle Scholar
- Brown MB, Forsythe AB: The small sample behavior of some statistics which test the equality of several means. Technometrics 1974, 16: 129–132. 10.2307/1267501View ArticleGoogle Scholar
- Welch BL: On the comparison of several mean values: An alternative approach. Biometrika 1951, 38: 330–336.View ArticleGoogle Scholar
- Cochran WG: Problems arising in the analysis of a series of similar experiments. Journal of Royal Statistics Society Serial C Applied Statistics 1937, 4: 102–118.Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.