Rank Difference Analysis of Microarrays (RDAM), a novel approach to statistical analysis of microarray expression profiling data
 Dietmar E Martin^{1},
 Philippe Demougin^{1},
 Michael N Hall^{1} and
 Michel Bellis^{2}Email author
DOI: 10.1186/147121055148
© Martin et al; licensee BioMed Central Ltd. 2004
Received: 22 July 2004
Accepted: 11 October 2004
Published: 11 October 2004
Abstract
Background
A key step in the analysis of microarray expression profiling data is the identification of genes that display statistically significant changes in expression signals between two biological conditions.
Results
We describe a new method, Rank Difference Analysis of Microarrays (RDAM), which estimates the total number of truly varying genes and assigns a pvalue to each signal variation. Information on a group of differentially expressed genes includes the sensitivity and the false discovery rate. We demonstrate the feasibility and efficiency of our approach by applying it to a large synthetic expression data set and to a biological data set obtained by comparing vegetativelygrowing wild type and tor2mutant yeast strains. In both cases we observed a significant improvement of the power of analysis when our method is compared to another popular nonparametric method.
Conclusions
This study provided a valuable new statistical method to analyze microarray data. We conclude that the good quality of the results obtained by RDAM is mainly due to the quasiperfect equalization of variation distribution, which is related to the standardization procedure used and to the measurement of variation by rank difference.
Background
In a typical microarray experiment, thousands of genes have their relative expression levels measured in parallel under different biological states [1, 2]. To identify differentiallyabundant genes, most published methods [3–5] progress through a similar sequence of elementary steps. First, a normalizing procedure is applied to make data sets comparable. If certain experimental conditions comprise several replicates, methods based either on parametric or nonparametric tests usually reduce the number of values generated by using their means. Then, gene variation is quantified by a statistic derived from intensity measurements. Knowledge of the null distribution of the gene variation, which is the distribution of its statistic when only random fluctuations occur, allows pvalues to be assigned to observed variations and genes to be ranked according to the significance of their variation. As the test is repeated as many times as there are genes, pvalues are corrected accordingly, and the false discovery rate (FDR, "the expected proportion of false positives among all genes declared significantly differentially expressed"[6]) is estimated.
We describe here, in detail, a new analysis method that has been used to analyze the transcriptome in yeast [7]. This method is original in several respects. First, Rank Difference Analysis of Microarrays (RDAM) replaces raw signal by its rank (R), expressed on a 0–100 scale, and we show that this simple transformation is a powerful normalizing procedure. Also, RDAM does not reduce replicated signals to their means, but instead only considers variations, expressed as rank differences (RD), between individual experimental points. An essential step is the standardization of RD observed between two replicates, permitting easy access to the empirical null distribution and allowing accurate and precise pvalues to be assigned to observed standardized RD (zRD). When dealing with replicated points, RDAM uses a random variable, the product of pvalues (ppv), for which the null distribution is straightforward to compute in a manner that is independent of the experimental conditions. Finally, RDAM estimates the total number of truly varying genes (TV), assigns a pvalue to each gene variation, characterizes the selection of a gene using the FDR and the percentage of truly varying genes included in the selection (sensitivity, S).
Analysis of synthetic data sets allowed us to specify the error distribution of all the estimators used (FDR, TV and S), and to demonstrate the strong predominance that the number of varying genes and the distribution of their variation have on the quality of the results.
We also analysed the transcriptional effects of the TOR2controlled signaling function using a genomewide microarray approach in yeast. In S. cerevisiae, TOR2 has two essential signaling functions. One, shared with TOR1, is required for translation initiation, transcription, and cell growth in response to the presence of nutrients [8–10]. The second is unique to TOR2, and functions in cellcycledependent actin polarization and possibly in transcription [8, 11]. A previous genetic screen for mutants defective in the TORshared and the TOR2unique functions identified several TOR2 temperaturesensitive alleles [12]. In this study, we compared total transcription profiles for strain SH121, which is specifically defective in the TOR2unique function, and its isogenic wild type counterpart SH100 [12].
Results
Standardization of positive variations
The simplest system to which our method can be applied comprises three experimental points, of which two are replicates, as described by the expression {Exp1A, Exp1B, Exp2A}, where the number refers to the biological condition and the final letter refers to the replicates. To identify significant variations in the comparison Exp2A vs. Exp1A, we have to first calculate the variation of gene i, VARi. From among several possibilities, we tested three different variation units: the fold change (FC), corresponding to the ratio of signals, the signal difference (SD), and the rank difference (RD). The RD uses a standardized signal measure that is independent of the scanner settings, because the signal is replaced by its rank, expressed on a 0–100 scale. This normalizing procedure consists of first calculating the absolute rank (AR) of each gene by ordering their signals from 0 to N (with the signals of all genes having a negative signal being set to zero, and N representing the number of nonnull and nonnegative signals) and then transforming the absolute rank value into a relative one (R = AR*100/N). In this way, all the signals are expressed on the same scale and are directly comparable.
We studied the variation distribution between the two replicates, i.e. Exp1A and Exp1B, reasoning that the observed empirical variation distribution would be an excellent approximation of the null distribution corresponding to the null hypothesis. The null hypothesis we have in mind states that all observed mRNA changes occurring under replicated conditions are due to a combination of biological and technological noise, and are not the result of any biologically significant process.
We tried to eliminate dependency of positive variation distribution on the signal rank by standardizing variations according to the general formula:
where VAR is to be replaced by any one of the variation units tested (SD, FC or RD). Using this expression, the sample mean and standard deviation (μVAR and stdVAR) were calculated for all genes having a rank within a given neighbourhood of R_{i}, the rank of gene i. This notation reflects the fact that the VAR distribution is not gene specific, but rank dependent.
Comparable results are obtained if the standardization is applied to the FC or the SD, but the zRD gives the best results in terms of distribution equalization: figure 2B shows, for example, that QQplots derived from zFC are more erratic than those derived from zRD, which are almost identical to the first diagonal up to the 99th percentile.
The fact that the zRD distribution is independent of rank can be explained by the fact that each gene's zRD_{i} follows the same zRD distribution. Therefore, we reasoned that the empirical cumulative frequency distribution, ecfd(zRD), approximates the distribution of zRD for any gene i under the null hypothesis, and we used Fo = 1  ecfd(zRD), based on a comparison between two replicated experiments, to assign a pvalue to any zRD calculated in a comparison between two different biological conditions (the pvalue is defined as the probability of zRD of an unchanged gene i to be equal to or greater than the observed zRD_{i} under the nulldistribution Fo). Because of the very large number of genes present on a chip, the null distribution is sampled a great number of times, generating a quasi continuous set of points that spans a wide range of values. This improves the precision of the Fo curve and allows accurate pvalues to be assigned for even large variations.
The entire procedure can then be applied to the comparison Exp2A vs. Exp1A. Because standardization curves constructed on the basis of the two replicates are used in the standardization process, the calculated zRD can be justifiably compared to the null distribution and interpolated on the Fo curve in order to assign a pvalue to each gene variation. At this step positive and negative variations can be processed together, although it is necessary to keep track of the actual type of variation, i.e. positive or negative, in order to conduct subsequent analysis.
Standardization of negative variations
False Discovery Rate (FDR), Total Variation (TV) and Sensitivity (S)
For each value of x, we estimate the False Discovery Rate,
and the sensitivity,
In the context of a transcriptome analysis, pvalues reflect how probable it is that a variation reaches or exceeds an observed value. Pvalues can always be used to rank genes, but the selection of significant variations in the context of multiple testing requires defining significance levels that are far more stringent than 0.01 or 0.05, as used in single testing. FDR and S parameters allow this difficulty to be overcome: for each c used as a potential critical value, S(c) reflects the fraction of trulyvarying genes that are selected by zRD>c, and FDR(c) estimates the fraction of selected genes that are likely to be invariant genes. It is therefore possible to plot FDR and S against c, and to construct, for positive and negative variations, what we call a selection abacus (figure 4).
Analysis of replicates
The simplest example of a replicated experimental scheme is the system {Exp1A, Exp1B, Exp2A, Exp2B}. While it would be tempting to average signals or ranks for each experimental condition and apply the method described above, this is not possible because averaging changes statistics and we have no practical way of obtaining the corresponding empirical null distribution.
In a first round of comparison, we conducted two analyses in parallel by applying RDAM to the first comparison Exp2A vs. Exp1A and to the second comparison Exp2B vs. Exp1B. Based on this first round of comparison, we obtained two pvalues for gene i: p1_{i} and p2_{i}. It could occur that gene i is detected as an increasing variation in the first comparison and as a decreasing variation in the second comparison. In this case, we apply a direction rule to decide on the final direction of variation. We consider simply that the lowest pvalue is in favor of its corresponding variation direction, and we set the pvalue of the discordant comparison to one. Once we have calculated and possibly corrected the pvalues, we construct a new random variable, the product of pvalues, ppv_{i} = p1_{i} × p2_{i}. To obtain an unbiased value for ppv, we apply the same procedure to a second round of comparison by exchanging Exp2A and Exp2B between the two comparisons, giving a second ppv value. The direction rule is applied to the two ppv before obtaining the final, averaged ppv.
The advantage of the random variable ppv is the ease of constructing its null distribution. In fact, the cfd(p_{i}H_{0}) is a uniform distribution over the interval [0,1]. Therefore, for cases in which two independent comparisons between two sets of duplicates were to be considered, we constructed two sets, U1 and U2, of 100000 points uniformly distributed over the interval [1,1], to take into account the possibility of increased and decreased variation for each point. These sets were randomized to make them independent in order to model the independence of measurement according to the null hypothesis. This hypothesis states that all variations are due to noise, and that for a particular gene all corresponding p values must be independent. Then, we apply the direction rule to the pair U1, U2 and calculate ppv for genes that are detected as increased. Thus, the F_{0} = 1  cfd(log10(1/ppv)) curve allows the significance of any value for ppv to be tested. The significance of ppv combines the significance of variation within each individual comparison and the significance of the correlation between these variations. F curve is, as usual, the observed 1  ecfd(log10(1/ppv)), and we get exactly the same kind of selector abacus, as shown in figure 4. Simulation of the ppv null distribution used exactly the same steps that the analysis process follows, i.e. application of the direction rule and construction of the product of pvalues, and resulted in a null distribution model we found appropriate seeing, both with experimental and synthetic data sets, that the observed distribution of ppv matches the null distribution when no variation occurs (data not shown).
The system {Epx1A, Exp1B, Exp2A, Exp2B} allows the construction of two sets of sandardizing curves, one from Exp1 replicates and the other from Exp2 replicates. As these curves are not equivalent, it is necessary to carry out both analyses and then use the more conservative one, whichever has the lower F curve.
The generalization of the entire procedure to more than two replicates is straightforward. For example, with three replicates {Exp1A, Exp1B, Exp1C, Exp2A, Exp2B, Exp2C}, there are 3 × 2 ways of arranging the experiments in order to obtain different sets of comparisons. Each round of comparison gives three pvalues for the gene i  p1_{i}, p2_{i} and p3_{i} – and the direction rule is applied the following way: in case the gene i is detected as an increasing variation in the first comparison and as a decreasing variation in the two other comparisons, we compare p1_{i} to p2_{i} × p3_{i} and determine the variation direction. Once we have calculated and possibly corrected the pvalues, we obtain the product of pvalues for gene i, ppv_{i} = p1_{i} × p2_{i} × p3_{i}. As the number of comparison rounds increases very rapidly with the number n of replicates (n!), we simply apply a circular permutation – the circular permutation of ABC consists in the subset of permutations ABC, BCA and CAB – to the replicates inside one of the biological condition which allows the number of rounds of comparison to be restricted to the number of replicates (n).
Generation of synthetic data
This entire procedure, composed of two successive steps, results in synthetic data sets of high quality, because the generation of data mimics the observed variation of genes. We compared the synthetic data sets Exp1C and Exp1D to the same natural data set Exp1A, and plotted the corresponding zRD of one comparison against the other, as shown in figure 5. We observed that these zRD were independent for genes that had not received an additional variation contribution, but were correlated for genes that had been changed. In general, however, this correlation was not absolute, except for some decreased genes. This phenomenon is explained by the high noise that characterizes weak signals: for such genes, there is a high probability that negative variation makes them reach the minimal rank value (zero). For high signals, there also exists a limit for the rank variation, but the noise is very small, and the truncation effect is not visible. We also observed that a small proportion of genes receiving an increased (decreased) variation contribution could be detected as increased (decreased) in one comparison, but decreased (increased) in another. All of these properties support the realistic nature of the synthetic data generated by our algorithm.
RDAM performances
We generated several synthetic data sets from the two experiments (Wt_t0a and Wt_t0b) by letting the number of increased and decreased genes equal 0, 100 or 500, and letting the maximum pvalue of extra variation equal 0.3, 0.1, 0.05 or 0.01.
Error on FDR estimation in the case of two replicates and 500 increased genes (p <= 0.1)
FDR  #1  #2  #3  #4  #5  #6  #7  #8  #9  #10 

5%  1 (12)  2 (8)  0 (3)  0 (3)  0  0 (0)  0 (0)  0 (0)  0 (0)  0 (2) 
10%  2 (73)  3 (25)  0 (14)  2 (37)  1 (2)  2 (20)  0 (0)  1 (12)  1 (20)  2 (20) 
15%  1 (160)  0 (84)  3 (68)  4 (59)  3 (5)  2 (58)  4 (41)  3 (40)  4 (63)  4 (57) 
20%  2 (217)  6 (157)  11 (130)  9 (129)  5 (65)  0 (134)  8 (152)  6 (81)  11 (146)  10 (136) 
Comparison of one step and two step procedures in the case of two replicates and 500 increased genes
Onestep analysis  Twostep analysis  

# Selected  # Target  Est FDR  True FDR  # Target  Est FDR  True FDR 
100  75  49  25  88  17  12 
200  133  58  34  167  23  17 
300  183  67  39  238  29  21 
400  229  75  43  298  35  26 
500  272  80  47  336  44  33 
Comparison with SAM on synthetic data
Comparison of RDAM and SAM in the case of two, three or five replicates and 500 increased genes
FDR50  FDR  True Positive  False Positive  

R  p var  RDAM  SAM  RDAM  SAM  RDAM  SAM  RDAM  SAM 
2R  0.3  55%  63%  33%    4    2   
0.1  25%  49%  21%    172    45    
0.05  9%  39%  11%  14%  281  23  35  4  
0.01  1%  26%  15%  23%  452  175  77  51  
3R  0.3  42%  47%  17%    55    11   
0.1  7%  28%  16%  20%  344  156  64  40  
0.05  3%  19%  13%  18%  409  246  60  54  
0.01  0%  11%  14%  22%  473  381  80  109  
5R  0.3  17%  21%  18%  25%  261  264  58  88 
0.1  2%  4%  15%  20%  455  395  81  101  
0.05  0%  2%  16%  20%  472  419  87  106  
0.01  0%  0%  19%  19%  474  443  111  105 
Analysis of the TOR experiment
Results of selection at S = 50% and FDR = 10%
Decreased  Increased  

FDR = 10%  S = 50%  FDR = 10%  S = 50%  
TV  Nb  S  Nb  FDR  TV  Nb  S  Nb  FDR  
MU vs WT t0  12  4  30%  7  13%  217  0  0%  196  45% 
MU vs WT t1  0  0  0%  0  0%  106  19  16%  73  27% 
MU vs WT t2  0  0  0%  0  0%  92  39  38%  51  11% 
t1 vs t0 WT  791  311  35%  459  14%  1164  817  63%  620  6% 
t2 vs t1 WT  231  0  0%  244  53%  233  45  17%  189  38% 
t2 vs t0 WT  1195  759  58%  641  8%  1275  1074  76%  660  3% 
t1 vs t0 MU  1158  644  50%  642  10%  1093  856  71%  564  3% 
t2 vs t1 MU  382  19  5%  351  44%  387  180  42%  226  14% 
t2 vs t0 MU  1289  725  51%  714  10%  1064  954  84%  543  2% 
Comparison with SAM on the TOR experiment
We ran our method in parallel with SAM in two situations displaying contrasting transcriptional responses. We tested a first comparison, Wtt1 vs Wtt0, which is characterized by a high number of varying genes and a good reproducibility between replicates, facilitating the detection of changes as reflected by the results of RDAM analysis which selected 620 increased genes at S = 50% and FDR = 6%. As SAM does not detect any increased genes at this selection level, we compared results obtained by the two methods at FDR = 10% and observed that among the 817 and 804 genes selected repectively by RDAM and SAM, only 426 were found in common. These results match what we found with synthetic data sets in case of high strength of variation (e.g. p <= 0.01 in figure 6D). We then considered comparisons of biological interest, i.e Mut1 vs. Wtt1 and Mut2 vs. Wtt2, and in this situation RDAM did not select any decreased genes and found only a few increased genes (see Discussion and Table 4). On the contrary SAM failed to detect any genes, either increased or increased.
Discussion
Analysis of the TOR experiment
RDAM is a method for identifying genes with changing expression levels using the userdetermined FDR and/or S selection parameters. This method was used to study the effects of a thermosensitive mutation of TOR2 in yeast. RDAM succeeded in identifying the few genes that are differentially regulated by the TOR mutation from among the entire mass of genes perturbed by the temperature shift. Recently it has been shown that TOR controls the translation of Gcn4 via the eIF4alpha kinase Gcn2 [14]. Under conditions of TOR inactivation by rapamycin, Gcn4 translation is enhanced, leading to the activation of Gcn4mediated transcription. Our data also demonstrate that TOR2 inactivation leads to enhanced transcription of Gcn4controlled target genes (biological results based on RDAM analysis are discussed in a forthcoming paper). Further experiments may show how the TOR2unique function is integrated into nutrient (or amino acid) responsive signaling pathways.
Normalizing of signal
Apart from randomlydistributed noise, microarrays are also prone to systematic effects that can bias the measurement of signal. All analysis methods are sensitive to systematic bias and include a preliminary normalizing step to make chips comparable. This is a limitation of this kind of approach, because the final result depends on the normalizing procedure used. Considering that all normalizing procedures rely on monotonous transformations that do not change the rank of raw data, we reasoned that if we used a statistics based on rank there would be no need to optimize the normalizing procedure. The rank unit we describe is similar to quantile normalization [15], but does not depend on the signal values of a particular chip as a reference: it can therefore be considered as an invariant. For example, if we focus specifically on the Affymetrix platform, we observe that the signal distribution changes with the different versions of the software: in MAS5, the 50th percentile is around 100, as compared to 1000 in MAS4. In our system, the rank of the genes at the same position in the signal distribution would not change, and would always be roughly equal to 50. This rank unit allows the drawing of plots in which all data are evenly distributed alongdimensions representing a signal. In addition, the linear density of points on the corresponding axes is constant, and the skewness of signal distribution has no effect on the graphical representation.
In our system, all values different from 0 or 100 are assigned to one and only one gene, because ordering of signals always delivers a series of contiguous rank values, even in cases of equivalent signal values. 0 is assigned to all unexpressed genes, as long as a robust method is available to detect them, and 100 to all genes for which the signal is saturated. In Affymetrix technology, especially with the scanner setup presently used, saturation is not a matter of concern and in our analysis, the value 100 is simply assigned to the highest signal. It is a complex problem to identify genes that are not expressed in a given experiment, and we decided to consider as absent only genes having a signal less than zero, as they occur in results delivered by MAS4.
Rank normalization results in transformation of the original signal distribution which is heavily skewed towards low values into a uniform distribution. As a consequence high rank variation could be assigned to small signal variations of weakly expressed genes, and it could be argued that our rank normalization method may bias variation detection towards genes with low signal. By using comparison between synthetic data sets we found no evidence of such a bias (data not shown).
Systematic usage of duplicates and standardization of variation
Our method has been developed within the framework of hypothesis testing and requires knowledge of the variation distribution for each gene when the null hypothesis is verified. The rationale of our approach considers replicated experiments as precisely representing a system in which all genes follow this hypothesis. However, it has long been recognized that the variation distribution expressed as a ratio or fold change is dependent upon the level of gene expression [16], and we show here that this property subsists when difference of signals or difference of ranks is used to measure variation. In theory it could be possible to use numerous replicates to obtain the empirical variation distribution of each gene. This is not possible for practical reasons, however, and we found that the classical centeredreduced standardization procedure can render variation distribution totally independent of gene expression level, as demonstrated by the QQplot analysis of figure 2A, and allow us to use duplicates to obtain the null variation distribution.
Independent analysis of positive and negative variations
In the algorithmic implementation of our method, we chose to proceed in two steps and deal with increased and decreased variations independently. First, this ensures that symmetrical comparisons (e.g. Exp1 vs. Exp2 and Exp2 vs. Exp1) give perfectly symmetrical results. Second, even if it were possible to devise another method that would allow one to proceed in one step, it seems more logical to consider increased and decreased variation separately. To clarify this point, a clear distinction must be made between up or downregulated mRNAs and increased and decreased variation. Regardless of the experimental points that are being compared, one always observes increased and decreased variations, but these variations have no absolute meaning because one only has to reverse the comparison to change the direction of variation. On the contrary, we can speak of up or down regulated mRNAs only when a causal effect exists, such as in a differentiation process or a kinetics or drug assay. In other words, a positive variation observed, for example, between two successive time points in a kinetic can be considered as an upregulation whatever its mechanism – gene or posttranscriptional regulation – but it is meaningless to invoke any particular form of regulation when comparing, for example, two unrelated cancer tissues. In the case of downregulated genes, the variation distribution of all weaklyexpressed genes is truncated, due to the impossibility of a decreasing signal crossing the zero line. In the case of upregulated genes, we do not observe the same effect for increasing variation of highlyexpressed genes, first because the signal distribution is heavily skewed towards low values, and second because the variance of highly expressed genes is very small (figure 1A). The observation that the reproducibility of variation was lower for downregulated genes than it was for upregulated genes [17] is partly explained by this reason, and we were able to demonstrate the statistical difference between up and downregulated genes by using synthetic data and observing that all FDR vs. S curves (figure 6) constructed with downregulated genes were lower than the corresponding curves for upregulated genes (data not shown). Moreover, we conducted a test showing that the joint analysis of increased and decreased variations degrades the quality of FDR estimation and reduces the number of true positives detected.
Replicates
We did not try to reduce the amount of raw data when using replicates, and devised a twostep method. A first statistics, the standardized rank difference zRD is constructed on each independent comparison, and pvalues are assigned by considering an empirical distribution that matches the null hypothesis. Then a second statistics, the product of pvalues ppv, is calculated and pvalues are assigned from the null distribution obtained by simulation. Simulation of the null distribution used exactly the same steps that the analysis process follows, i.e. application of the direction rule and construction of the product of pvalues, and resulted in a null distribution model we found appropriate seeing, both with experimental and synthetic data sets, that the observed distribution of ppv matches the null distribution when no variation occurs (data not shown). It turned out that this scheme is very flexible and of general applicability: because the second step is rooted in a rigorous statistical method that uses only pvalues as input data, it is possible to adapt or to improve the entire process simply by focusing on the first step of pvalue estimation. For example, to apply our method to cDNA glass arrays, the only step to be modified would be the variation standardization. Alternatively we could use the segmental approach proposed by Yang and colleagues [18], which is claimed to equalize log ratio distribution, or the variance stabilization method of Huber et al [19], which is efficient in equalizing variation distribution of transformed intensity measurements in both cDNA and oligonucleotide platforms.
FDR, Total Variation and Sensitivity Estimation
The way in which we estimated FDR is exactly the same as that suggested by B. Efron et al. in their demonstration of the equivalence of empirical Bayes and frequentist approaches (Efron B, Storey J. D. and Tibshirani R. http://wwwstat.stanford.edu/~tibs/ftp/bradfdr.pdf, see equation 3.8 and [20]). We did, however, use another heuristic approach to estimate TV because we observed that the estimator proposed by Storey et al [20] could be very difficult or impossible to calculate when the expression of a small fraction of genes changes. We demonstrated using synthetic data that our estimator was not prone to this type of instability (not shown), and that under realistic conditions (additional variation of p <= 0.10) our estimate was 60%, 65% and 80% of the true TV in the case of two, three and five replicates, respectively. The accuracy of this estimator is obviously dependent on the power of the test, which is itself under the control of the number of replicates. We have also shown that estimated sensitivity was biased by a constant factor that was mostly determined by the error made in TV estimation. Finally, it must be emphasized that the error made in TV estimation has little effect on FDR estimation, as demonstrated by forcing RDAM to use the true TV and K values during the process of synthetic data analysis (data not shown).
Synthetic data sets
Using the empirical noise distribution observed between two replicates, we devised a method for constructing synthetic data sets. Most published methods add noise to a signal that is supposed to represent the true signal of the gene. We showed here that raw signals without denoising could be used and gave excellent result as judged both by the final distribution of signals and by indirect controls such as the preservation of variation distribution and the possibility of successfully analyzing synthetic data substituted for the original data.
Synthetic data sets are well adapted for judging the respective performances of different analysis methods. To characterize the scoring procedure of a particular method, we used a new type of diagram that plots FDR vs. S, two quantities that relate to the subset of selected genes and that seem better adapted than Receiver Operator Characteristic (ROC, [21]), which relates to both selected and rejected genes (FDR vs. specificity). We proposed using FDR_{50}, the FDR at S = 50%, as a comparative index between different methods and showed that RDAM has an FDR_{50} that is 30 percentage points smaller than SAM in the case of three replicates and applied changes with p <= 0.10 (figure 6).
Conclusions
RDAM is a new statistical method whose performances have been precisely evaluated through extensive analysis of synthetic data sets. When applied to TOR experiment, our method succeeded in finding the few genes of biological interest which were concealed in the mass of varying genes induced by the temperature shift. Comparison with SAM showed that our method obtained the same (if not better) results but with a smaller consumption of chips We conclude that the good quality of the results obtained by RDAM is mostly due to the use of replicates to calibrate the noise and to the quasiperfect equalization of variation distribution, which is related to the standardization procedure used and to the measurement of variation by rank difference.
Methods
Preparation of RNA
Saccharomyces cerevisiae strains SH100 and SH121 [12] were grown overnight in yeast extract peptone glucose (YPD), diluted to an optical density measured at 600 nm of 0.05 (OD600 = 0.05), and grown for an additional 4 hours the next day. The main cultures were then inoculated in YPD medium and grown at 25°C or shifted to 37°C for 2 or 6 hours. All cultures were grown as independent duplicates and were harvested at a final OD600 of 0.8 to 0.9 to minimize the influence of differences in growth phase.
Upon harvesting by centrifugation (2 min, 3000 × g) at 4°C, cells were washed once in icecold water, centrifuged again, and the cell pellet was flash frozen in liquid nitrogen. Total RNA was extracted using a hot phenol method essentially as described by Schmitt, M.E. et al. [22].
Microarray hybridization
Affymetrix™ S98 Yeast Genome GeneChips, containing 6,400 S. cerevisiae (S288C strain) genes and 600 additional probe sets representing putative open reading frames [23], were used throughout this study. Synthesis of cDNA and in vitro transcription of biotinlabeled cRNA, as well as microarray hybridisation, washing and staining procedures, were carried out according to standard protocols as recommended by the manufacturer. Two independent preparations were used for each experimental point.
Data processing
The scanned microarray images were analysed using the algorithm implemented in MAS 5.0 (Affymetrix, Santa Clara, CA) and the generated raw data were further processed by scripts written in Matlab language (MathWorks, Natick, MA.). SAM analysis [3] of synthetic data was made using version 1.21 of the program [24] with the following default parameters: unlogged data, number of permutations set to 100 and "KNearest Neighbors Imputer" used.
Raw data files were uploaded to NCBI's GEO repository under the series number GSE1814 http://www.ncbi.nlm.nih.gov/geo/.
Abbreviations
 RD:

rank difference
 zRD:

standardized rank difference
 FDR:

false discovery rate, S, sensitivity
 TV:

total variation
Declarations
Acknowledgements
We thank Estelle Chanudet for her participation in the development of the synthetic data generator, Reinhold Koch and Gunnar Wrobel for helpful discussions, and Lydia Michaut for her comments on the manuscript. We thank anonymous reviewers for their useful comments. This work was supported by grants from the Genopole of Montpellier to MB, and from the Canton of Basel and the Swiss National Science Foundation to MNH.
Authors’ Affiliations
References
 Schena M, Shalon D, Davis RW, Brown PO: Quantitative monitoring of gene expression patterns with a complementary dna microarray. Science 1995, 270: 467–70.View ArticlePubMedGoogle Scholar
 Lipshutz RJ, Fodor SP, Gingeras TR, Lockhart DJ: High density synthetic oligonucleotide arrays. Nature Genetics 1999, 21: 20–4. 10.1038/4447View ArticlePubMedGoogle Scholar
 Tusher VG, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci USA 2001, 98: 5116–5121. 10.1073/pnas.091062498PubMed CentralView ArticlePubMedGoogle Scholar
 Pan W: A comparative review of statistical methods for discovering differentially expressed genes in replicated microarray experiments. Bioinformatics 2002, 18: 546–554. 10.1093/bioinformatics/18.4.546View ArticlePubMedGoogle Scholar
 Hatfield GW, Hung SP, Baldi P: Differential analysis of DNA microarray gene expression data. Mol Microb 2003, 47: 871–877. 10.1046/j.13652958.2003.03298.xView 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 BMethodological 1995, 57: 289–300.Google Scholar
 Williams RM, Primig M, Washburn B, Winzeler EA, Bellis M, Sarrauste de Menthière C, Davis RW, Esposito RE: The Ume6 regulon coordinates metabolic and meiotic gene expression in yeast. Proc Natl Acad Sci USA 2002, 99: 13431–13436. 10.1073/pnas.202495299PubMed CentralView ArticlePubMedGoogle Scholar
 Schmelzle T, Hall MN: TOR, a central controller of cell growth. Cell 2000, 103: 253–62. 10.1016/S00928674(00)001173View ArticlePubMedGoogle Scholar
 Barbet NC, Schneider U, Helliwell SB, Stansfield I, Tuite MF, Hall MN: TOR controls translation initiation and early G1 progression in yeast. Mol Biol Cell 1996, 7: 25–42.PubMed CentralView ArticlePubMedGoogle Scholar
 Loewith R, Jacinto E, Wullschleger S, Lorberg A, Crespo JL, Bonenfant D, Oppliger W, Jenoe P, Hall MN: Two TOR complexes, only one of which is rapamycin sensitive, have distinct roles in cell growth control. Mol Cell 2002, 10: 457–468. 10.1016/S10972765(02)006366View ArticlePubMedGoogle Scholar
 Schmidt A, Bickle M, Beck T, Hall MN: The yeast phosphatidylinositol kinase homolog TOR2 activates RHO1 and RHO2 via the exchange factor ROM2. Cell 1997, 88: 531–42. 10.1016/S00928674(00)818930View ArticlePubMedGoogle Scholar
 Helliwell SB, Howald I, Barbet N, Hall MN: TOR2 is part of two related signaling pathways coordinating cell growth in Saccharomyces cerevisiae. Genetics 1998, 148: 99–112.PubMed CentralPubMedGoogle Scholar
 Natarajan K, Meyer MR, Jackson BM, Slade D, Roberts C, Hinnebusch AG, Marton MJ: Transcriptional profiling shows that Gcn4p is a master regulator of gene expression during amino acid starvation in yeast. Mol Cell Biol 2001, 21: 347–68. 10.1128/MCB.21.13.43474368.2001View ArticleGoogle Scholar
 Cherkasova VA, Hinnebusch AG: Translational control by TOR and TAP42 through dephosphorylation of eIF2 alpha kinase GCN2. Genes Dev 2003, 17: 859–72. 10.1101/gad.1069003PubMed CentralView ArticlePubMedGoogle Scholar
 Bolstad BM, Irizarry RA, Astrand M, Speed TP: A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 2002, 19: 185–193. 10.1093/bioinformatics/19.2.185View ArticleGoogle Scholar
 Chen Y, Dougherty ER, Bittner ML: Ratiobased decisions and the quantitative analysis of cDNA microarray images. J Biomedical Optics 1997, 2: 364–374. 10.1117/1.429838View ArticlePubMedGoogle Scholar
 Cohen P, Bouaboula M, Bellis M, Baron V, Jbilo O, Poinot Chazel C, Galiegue S, Hadibi EH, Casellas P: Monitoring cellular responses to Listeria monocytogenes with oligonucleotide arrays. J Biol Chem 2000, 275: 11181–11190. 10.1074/jbc.275.15.11181View ArticlePubMedGoogle Scholar
 Yang H, Haddad H, Tomas C, Alsaker K, Papoutsakis ET: A segmental nearest neighbor normalization and gene identification method gives superior results for DNAarray analysis. Proc Natl Acad Sci USA 2003, 100: 1122–1127. 10.1073/pnas.0237337100PubMed CentralView ArticlePubMedGoogle Scholar
 Huber W, Von Heydebreck A, Sultmann H, Poustka A, Vingron M: Variance stabilization applied to microarray data calibration and to the quantification of differential expression. Bioinformatics 2002, 18(Suppl 1):S96S104.View ArticlePubMedGoogle Scholar
 Storey JD, Tibshirani R: Statistical significance for genomewide studies. Proc Natl Acad Sci USA 2003, 100: 9440–9445. 10.1073/pnas.1530509100PubMed CentralView ArticlePubMedGoogle Scholar
 Khodarev NN, Park J, Kataoka Y, Nodzenski E, Hellman S, Roizman B, Weichselbaum RR, Pelizzari CA: Receiver operating characteristic analysis: a general tool for DNA array data filtration and performance estimation. Genomics 2003, 81: 202–209. 10.1016/S08887543(02)000423View ArticlePubMedGoogle Scholar
 Schmitt ME, Brown TA, Trumpower BL: A rapid and simple method for preparation of RNA from Saccharomyces cerevisiae. Nucleic Acids Res 1990, 18: 3091–92.PubMed CentralView ArticlePubMedGoogle Scholar
 GeneChip Yeast Genome S98 Array[http://www.affymetrix.com/support/technical/datasheets/yeast_datasheet.pdf]
 SAM: Significance Analysis of Microarrays[http://wwwstat.stanford.edu/~tibs/SAM/]
 Chaudhuri P, Marron S: SiZer for exploration of structures of curves. JASA 1999, 94: 807–823.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an openaccess 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.