The statistics of identifying differentially expressed genes in Expresso and TM4: a comparison

Background Analysis of DNA microarray data takes as input spot intensity measurements from scanner software and returns differential expression of genes between two conditions, together with a statistical significance assessment. This process typically consists of two steps: data normalization and identification of differentially expressed genes through statistical analysis. The Expresso microarray experiment management system implements these steps with a two-stage, log-linear ANOVA mixed model technique, tailored to individual experimental designs. The complement of tools in TM4, on the other hand, is based on a number of preset design choices that limit its flexibility. In the TM4 microarray analysis suite, normalization, filter, and analysis methods form an analysis pipeline. TM4 computes integrated intensity values (IIV) from the average intensities and spot pixel counts returned by the scanner software as input to its normalization steps. By contrast, Expresso can use either IIV data or median intensity values (MIV). Here, we compare Expresso and TM4 analysis of two experiments and assess the results against qRT-PCR data. Results The Expresso analysis using MIV data consistently identifies more genes as differentially expressed, when compared to Expresso analysis with IIV data. The typical TM4 normalization and filtering pipeline corrects systematic intensity-specific bias on a per microarray basis. Subsequent statistical analysis with Expresso or a TM4 t-test can effectively identify differentially expressed genes. The best agreement with qRT-PCR data is obtained through the use of Expresso analysis and MIV data. Conclusion The results of this research are of practical value to biologists who analyze microarray data sets. The TM4 normalization and filtering pipeline corrects microarray-specific systematic bias and complements the normalization stage in Expresso analysis. The results of Expresso using MIV data have the best agreement with qRT-PCR results. In one experiment, MIV is a better choice than IIV as input to data normalization and statistical analysis methods, as it yields as greater number of statistically significant differentially expressed genes; TM4 does not support the choice of MIV input data. Overall, the more flexible and extensive statistical models of Expresso achieve more accurate analytical results, when judged by the yardstick of qRT-PCR data, in the context of an experimental design of modest complexity.


Background
DNA microarrays are a powerful means of monitoring the expression of thousands of genes simultaneously. A variety of computational and statistical methods have been proposed to extract information from the large quantity of data generated from microarray experiments. Many methods assume, as we do here, the use of cDNA labeled with one of two fluorescent dyes to differentiate two treatments on a single microarray, implying data from two images to be analyzed. These methods include a number of data normalization techniques to reduce the effects of systematic errors and various kinds of statistical tests to identify differentially expressed genes in comparisons among different experimental conditions. There is as yet no single method that can be recommended under all circumstances for either normalization or identification of differential gene expression.
In recent years, ANOVA methods have gained popularity for identification of differential gene expression. The power of ANOVA methods derives from their flexibility in fitting and comparing different models to a given set of data [1]. One such method is the two-stage, log-linear ANOVA mixed models technique of Wolfinger, et al., [2]. Its first stage uses a normalization model designed to remove global effects across all microarrays. Its second stage uses a gene-specific model to estimate gene-treatment interactions as ratios of gene expression under control and treated conditions, along with a statistical significance. Kerr [3] notes that the global normalization model employed in this technique is conducive to combining data across genes for realistic and robust models of error, especially when random effects are included. Pan [4] compares different microarray statistical analysis methods and demonstrates that the log-linear ANOVA mixed model approach performs better than the t-test and regression approaches. The regression approach, although flexible and robust, assumes that the data is drawn from a normal distribution, while the t-test is limited due to very few degrees of freedom. Chu, et al., [5] compare two loglinear ANOVA mixed models for probe-level, oligonucleotide array data and found that both types of models capture key measurable sources of variability of oligonucleotide arrays for real and simulated data. Cui and Churchill [6] review the use of a mixed ANOVA model for analyzing a cDNA microarray experiment and conclude that such models provide a powerful way to obtain information from experiments with multiple factors or sources of variation. Rosa, et al., [7] review issues of analyzing cDNA microarrays with mixed linear models and puts such analysis in the larger context of Bayesian analysis procedures and adjustments for multiple testing.
Data normalization is the first step in analyzing microarray data; numerous data normalization methods have been proposed and investigated. While refinements of existing methods continue to appear (e.g., Futschik and Crompton [8]), naive methods, such as total intensity normalization, are still in use (e.g., Held, et al., [9]). Xie, et al., [10] did a comparative study of normalization methods and test statistics to analyze the results of a DNAprotein binding microarray experiment. Using performance and bias correction criteria, Bolstad, et al., [11] evaluate the cyclic lowess method, the contrast method, the quantile method, and baseline array scaling methods, both linear and non-linear; they demonstrate that normalization methods incorporating data from all microarrays perform better than methods employing a baseline array.
Several software tools that combine data normalization and statistical analysis are currently available. Dudoit, et al., [12] review these software tools with an emphasis on the TM4 microarray software suite, Bioconductor in R, and the BioArray Software Environment (BASE) system. Saeed, et al., [13] describe the features and capabilities of TM4, while Quackenbush [14] describes the normalization and transformation methods implemented in it. Williams, et al., [15], Zhu, et al., [16], and Khaitovich, et al., [17] have used TM4 in microarray data analysis. Another system is Expresso, an experiment management system that serves as a unifying framework to study data driven applications such as microarray experiments [18][19][20]. Expresso has adapted the two-stage ANOVA mixed models technique of Wolfinger, et al., [2] to the particular needs of individual microarray data sets. Our experience with numerous such data sets has demonstrated that modeling the underlying experiment carefully and completely is essential to obtaining meaningful and defensible results. Use of tools that require experiments to conform to their analysis methods are less than satisfactory.
In this paper, we compare the Expresso analysis methodology to the approach provided in the TM4 microarray analysis software suite [13]. Each is invoked to identify differentially expressed genes in two experimental data sets, each of which uses an Arabidopsis thaliana oligonucleotide array. Along the way, we demonstrate differences between the use of integrated intensity values (IIV) and median intensity values (MIV) as inputs. We report interactions between normalization and gene identification methods. We use quantitative reverse-transcriptase PCR (qRT-PCR) results to assess the consistency of genes reported by TM4 and Expresso as having significant differential expression.

Results and discussion
Here, we report a portion of the results obtained in our comparison of Expresso analysis and the TM4 pipeline (see Materials and Methods). Figure 1 illustrates the overall flow of the statistical analyses of microarray data that were done in this study. We began with microarray data in GPR format from Experiment 1 and Experiment 2.
Median intensity values (MIV) from the GPR files can be analyzed by the Expresso GP and GOT models directly. ExpressConverter provides integrated intensity values (IIV) for further Expresso and TM4 analysis. The MIDAS normalization and filtering pipeline executes these steps in order: total intensity normalization (subscript T), lowess normalization (subscript L), standard deviation regularization (subscript S), and low intensity filter (subscript F). MIDAS allows tapping the output of any step in the pipeline; for example, IIV TL signifies an MEV file after total intensity normalization followed by lowess normalization. The identification of genes with significant differential expression was performed on all GPR and MEV files, using the Expresso GP and GOT models and the t-test in MEV.

Normalization and low intensity filtering in TM4
Quackenbush [14] describes the use of ratio-intensity plots (RI-plots) to detect and normalize for any systematic intensity-dependent dye bias using lowess normalization (see Materials and Methods). We evaluated the effect of lowess normalization within the context of the flow in Figure 1 by creating RI-plots after each step for the second replicate microarray in Experiment 1, WT plant. Supplementary Figure 1 -see Additional file: 1 contains these RIplots. The IIV TL is indeed effective, for this data set, in correcting systematic dye bias, suggesting that preprocessing by these two normalization steps in MIDAS may be a good practice in many situations.
The normalization and filtering pipeline affects the number of genes identified as differentially expressed in both the GP and GOT models. See Table 1. For example, in Experiment 1, the GP model using IIV input data identifies 567 up-expressed genes in the WT microarrays, while it identifies only 460 WT genes as up-expressed if IIV TLSF (processed by the complete MIDAS pipeline) input data is used.
Small changes in the number of genes identified as up-or down-expressed after successive MIDAS steps may mask larger changes in the composition of sets of up-and down-expressed genes. To obtain a more precise view of the effects of MIDAS changes, we computed retention counts (RC) and retention percentages (RP) between the gene successive sets whose numbers are in Table 1. RC is the number of genes in the set before the MIDAS step that remain in the set after the step. RP is the percentage of remaining genes with respect to the number of genes in the set after the MIDAS step. Table 2 contains the RC and RP values corresponding to the counts in Table 1. For Experiment 1, there is a tremendous drop in retention during the lowess normalization that follows the total intensity normalization. There is not a drop of corresponding magnitude for Experiment 2. For both experiments, normalization has a significant effect on the sets of genes identified as differentially expressed.
In Experiment 1 results, the number of genes commonly assessed by Expresso as significantly expressed when using IIV and IIV T is high. For example, there is 95.45% retention of WT genes (545 total) assessed as up-expressed when using IIV T data in Expresso compared to that when using IIV data. Retention percentage of these genes assessed as expressed however went down after doing lowess normalization. There is only 15.20% retention of WT genes (74 total) assessed as up-expressed in the results when using IIV TS data in Expresso compared to when using IIV T . While we observe increase in the retention percentage in IIV TLS (from IIV TL ) and IIV TLSF (from IIV TLS ), there's low retention percentage in the results using IIV TLSF from IIV data. This can be traced in the low retention percentage of IIV TL from IIV T . Hence, the normalization method that affects the results in Experiment 1 the most is lowess normalization.
The results of Expresso on Experiment 2 show that low retention percentages happen after appli-cation of total intensity normalization (lowest is 59.30%) and after application of low intensity filtering (lowest is 61.19%). The low retention percentages shown in the IIV ∩ IIV TLSF column implies that the normalization pipeline also significantly affects the results in Expresso analysis of Experiment 2.

Choice of intensity signal data
The input to statistical analysis of microarray experiments is a set of real numbers that represent the measured intensity signal for each spot in a microarray. Much statistical analysis of microarray data has traditionally used median intensity values (MIV). The alternative used in TM4 is the integrated intensity value (IIV). (See Materials and Methods.) Since IIV is intended to integrate the measured intensity across the biological sample printed at a spot, one might expect IIV to be a more accurate assessment of the biological measurement than MIV data. For example, a spot having 100 pixels and a median intensity of 5,000 has the same IIV as a spot having 50 pixels and a median intensity of 10,000.
This study provides the opportunity to observe the difference that choosing MIV or IIV makes on the sets of genes ultimately identified as differentially expressed. We used the GP model to analyze unnormalized MIV and IIV data from Experiment 1, and we used the GOT model to analyze unnormalized MIV and IIV data from Experiment 2.
Overall flow of the statistical analyses of microarray data Figure 1 Overall flow of the statistical analyses of microarray data. Data input is in GPR format and provides the MIV for each spot. TM4 analysis requires ExpressConverter to generate MEV format containing the IIV for each spot. The normalization steps performed by MIDAS are T (total intensity normalization), L (lowess normalization), S (standard deviation regularization), and F (low intensity filter). Differential gene expression is obtained from the Expresso GP model, the Expresso GOT model (Experiment 2 only), and the t-test in MEV.  Table 3 reports a summary of the results. In Experiment 1, 725 WT genes are assessed as up-expressed and 774 WT genes are down-expressed when MIV data are used in Expresso. These numbers decreased to 567 up-expressed genes and 552 down-expressed genes when IIV data are used instead. A similar trend is observed in Experiment 2 results when using MIV and IIV data. These results suggest that employing IIV input data with Expresso analysis leads to more conservative results than employing MIV input data.

Comparison of statistical methods
We compared the performance of the GP model, the GOT model, and the t-test of MEV in identifying differentially expressed genes in Experiment 2. We used the IIV TLSF data of Experiment 2 as input to these methods. We also contrast these results with MIV data analyzed in GP. Table 4 reports counts for these analyses.
The plot in Figure 2a demonstrates that the estimates of Iog 2 (fold change) are the same in GP and GOT. As Figure  2b shows, the p values by GOT are smaller compared to the p values calculated by GP. Use of the MEV t-test resulted in fewer genes assessed as significantly expressed when compared to the numbers for the GP model. The results obtained when MIV data was used as input to GP, is closest to the results when using IIV TLSF .
To compare the effectiveness of Expresso and TM4 in identifying gene differential expression, we compared the identified direction of differential expression of a select set of genes per genotype in Experiment 2 with results obtained by qRT-PCR. See Table 5. The lowest overall percentage (71.9%) of agreement is between the qRT-PCR results and the MEV t-test results using IIV TLSF . The log(fold change) estimates of the GP model has 77.1% percentage agreement with the qRT-PCR results, which is slightly higher than the percentage for the MEV t-test. The results of the GP model using MIV data demonstrated the greatest agreement, 90.1%, with the qRT-PCR results. For each gene, a histogram of the Iog 2 (fold change) estimates of qRT-PCR, the GP model using MIV, the GP model using IIV TLSF , and the MEV t-test is given. The 50 histograms are spread over three figures to enhance readability and are in increasing order by qRT-PCR estimated change. In general, the log 2 (fold change) estimates of the GP model and of the MEV t-test, all IIV TLSF input data, are approximately the same, while being slightly different from estimates of the GP model using MIV input data. As might be expected, disagreement between qRT-PCR and microarray results are more prevalent for small estimated log 2 (fold change)

Numbers of identified differentially expressed genes in the GP and GOT models after preprocessing by 0 or more MIDAS pipeline steps -IIV, IIV T , IIV TL , IIV TLS , or IIV TSLF . For Experiment 1, up-expression (+) and down-expression (-) numbers are given for both WT and antiPLD. For Experiment 2, + and -1 numbers are given for all 4 genotypes separately. The ALL entries correspond to the H1 hypotheses of the GOT model (see Materials and Methods).
Genotype IIV  IIV T  IIV TL  IIV TLS  IIV TLSF GP values. The histograms for genes AT4G09020 (Figure 3), AT1G35580 (Figure 4), and AT3G29360 ( Figure 4) show that the direction of log(fold change) estimate of qRT-PCR matches the direction of the GP model using MIV input data, while differing from the direction of the estimates of the GP model and the MEV t-test using IIV TLSF input data.

Conclusion
Our integration and comparison of Expresso analysis and the capabilities of TM4 has highlighted successes in microarray analysis, some similarities, and some differences. The success of microarray analysis is demonstrated by considerable agreement between qRT-PCR results and the results of all the examined microarray analysis methods. The greatest agreement was found when median intensity value (MIV) inputs were analyzed with the Expresso GP analysis model. We also found that the use of integrated intensity value (IIV) inputs for Expresso analysis consistently resulted in fewer genes identified as differentially expressed when compared to results from MIV inputs. This suggests that the use of IIV inputs is more conservative than the use of MIV inputs, while MIV inputs may give greater agreement to qRT-PCR results than IIV inputs.
Our results demonstrate that the MIDAS normalization and filtering pipeline corrects systematic intensitydependent dye bias on a per microarray basis. The normalization stage in Expresso analysis removes global effects across all microarrays and complements the per microarray normalization methods of MIDAS. The generally better agreement of Expresso analysis with qRT-PCR results when compared to the MEV t-test suggests that it would be desirable for MEV to have an ANOVA test that has the greater flexibility of the Expresso gene model.

Median and integrated intensity values
This research considers two ways of measuring spot intensity, one or both of which are reported by typical microarray image processing software. The median intensity value  IIV T  IIV TL  IIV TLS  IIV TLSF  IIV  (MIV) of a spot is the median value of all the pixels identified as part of the spot. The integrated intensity value (IIV) of a spot is the total value of all the pixels identified as part of the spot. In this research, both are backgroundcorrected values. If the IIV data is unavailable, but the radius of the bounding circle of the spot and its average intensity value are available, then the IIV data can be estimated as the product of the average intensity and the number of pixels in the circle. This is the estimate used by ExpressConverter (below) when it converts GPR format data to MEV format data.

Microarray data sets
We used data sets from two experiments that used the Arabidopsis Oligonucleotide Microarrays [21], which include 25   of an antisense plant for phospholipase D α (antiPLD) in Columbia background to drought stress [22]. Plants were harvested at a single time point, and two biological repli-  Table 4 (a) This is a scatter plot of the estimated log 2 (fold change) values from the GP and GOT models; these values are essentially identical, (b) This is a scatter plot of the -log 10 (pvalue), again for the GP and GOT models. The dotted lines correspond to -log 10 (0.05), as our significance cutoff is 0.05. For the GP model, it is the points to the right of the dotted line that are significant. For the GOT model, it is the points above the dotted line that are significant.    -2

GP
The microarray results of Experiment 2 suggested that exposure to elevated CO2 resulted in changes in expression of many genes associated with carbon metabolism, and those associated with pho-tosynthetic carbon metabolism in particular. This included genes encoding proteins that transport pho-tosynthate out of the chloroplast, where carbon fixation takes place, for export to other parts of the cell, and also genes encoding transport proteins that export carbon skeletons out of the cell to other tissues where growth is taking place. Because of this finding, it was important to validate the results obtained for gene expression associated with carbon metabolism with qRT-

Expresso analysis
Expresso analysis employs a general and flexible method to identify differentially expressed genes that is adapted from the two-stage analysis method of Wolfinger, et al., [2]. In general, Expresso analysis consists of two log-linear ANOVA mixed models, called the normalization model and gene model. The first estimates and removes the experiment-wise systematic errors, while the second estimates and removes the gene specific errors. The residual that remains is the log-ratio estimate for each gene. In particular, the Tukey-Kramer multiple comparison of treatment effects on each gene is performed to estimate its expression level and the significance of (confidence in) that expression level. Expresso analysis is implemented for and executed on SAS (SAS/STAT version 8.2, SAS Institute Inc., Gary, NC USA).
The original model of Wolfinger, et al., [2] includes the treatment and the array as the main effects. In previous Expresso analysis, we have extended that model to experiment-appropriate models that include additional fixed and random effects. Here, the design of the two-dye oligonucleotide microarray used in Experiment 1 and Experiment 2 includes various controls strategically positioned in different blocks of the microarray. This makes it possible to estimate the random block effect in each microarray. Furthermore, the dye effect is included in the normalization model to estimate and remove the global dye bias.
For this research, we developed two Expresso models, one whose gene model assesses the gene-(plant sample) effect (the GP model) and the other whose gene model assesses the gene-genotype-treatment effect (the GOT model). The GP model is much like previous Expresso models and is applicable to both Experiment 1 and Experiment 2. However, the GOT model is specific to analyzing Experiment 2. In both experiments, we used the GP model to estimate Comparison of qRT-PCR results, Part III qRT-PCR GP model -MIV GP model -IIV TLSF MEV t-test -IIV TLSF the differences in response of individual genotypes to treatment (drought stress versus control or ozone stress versus ambient ozone). However, the GOT model was used to estimate the effects of treatment (ozone stress), aside from the effect of individual genotypes.

Expresso GP model
The normalization model is This differs from the normalization model in Wolfinger, et al., [2], in that it incorporates dye and block effects. It is a refinement of the Expresso normalization model in Watkinson, et al., [18], in that it has no printing pin effect, which is specific to the analysis in the earlier paper, and includes the block effect.
The second stage of the analysis uses the residual values r spdab computed in the first stage to estimate the interaction between an individual gene g and each plant sample p at a significance level ≤ α = 0.05. Index g is added to the residual values r spdab resulting to r gspdab . The value of g is determined using the mapping of s index values to g index values. The gene model is Here, g is a spot that represents a gene (not a blank or control) within the dye d image in block b of array a. The value G g is the mean of residual values for all spots that represent gene g in all images. The interactions (G × P) gp of gene g with plant sample p and of (G × D) gd of gene g with dye d are the fixed effects. The interaction (G × A) ga of gene g with mi-croarray a is a random effect. The λ gspdab values are stochastic errors. This differs from the gene model in Wolfinger, et al., [2], in that it incorporates interactions between gene and dye and between gene and array. It refines the Expresso gene model in Watkinson, et al., [18] to include the interaction between gene and array. The difference (G × P) gi -(G × P) gj is the estimate of the log 2 (fold change) of gene g in the experimental comparison P i versus P j . The analysis also yields a p-value for the statistical confidence in each difference.

Expresso GOT model
We wanted to estimate the gene-treatment effects separately from gene-genotype-treatment interaction effects using just one model. To do this, we designed the genegenotype-treatment model, an alternative set of log-linear ANOVA mixed models, for the elevated CO 2 experiment where we unfold the genotype information from the plant sample factor in the GP model. This resulted in a normalization model that includes the genotype effect (O o ) with 4 levels (Col-0, Cvi-0, WS, and Th) and basic treatment effect (T t ) with 2 levels (test and control). The random array (A a ) effect however needs to be removed from the model since it confounds the genotype effect.
The normalization model is Each y sotdab value is the log 2 -transformed intensity of spot s for genotype o and treatment t within the dye d image in block b of array a. We have that μ is as in the GP model.  [2], in that it incorporates genotype (organism), dye, and block effects. It is a refinement of the Expresso normalization model in Watkinson, et al., [18], in that it has no printing pin effect, and includes the genotype and block effects.
The second stage of the analysis uses the residual values r sotdab computed in the first stage to estimate the interaction among an individual gene g, each genotype o, and each treatment t, at a significance level ≤ α = 0.05. Index g is added to the residual values r sotdab resulting to r gsotdab . The value of g is determined using the mapping of s index values to g index values. The gene model is Here, G g are as for the GP model.  [2], in that it incorporates interactions between gene and genotype, between gene and genotype with treatment, and between gene and array. It refines the Expresso gene model in Watkinson, et al., [18] to include interactions between gene and genotype and between gene and genotype with treatment.
We do pairwise comparisons as in the GP model, but there are two plausible classes of null hypotheses to test within the GOT gene model. If τ is the number of treatments, then the class 1 null hypothesis for gene g and comparison i,j, where 1 ≤ i <j ≤ τ, is The difference (G × T) gi -(G × T) gj is the estimate of the log 2 (fold change) of gene g in the T i versus T j comparison. This particular comparison looks for gene-treatment effects that are independent of genotype.
We can still estimate the expression level of each gene with respect to a specific genotype level by computing the pair-wise least square differences of gene-genotype-treatment interaction effects. The class 1 null hypothesis for gene g and comparison i, j, where 1 ≤ i <j ≤ τ, is The difference (G × O × T) goi -(G × O × T) goj is the estimate of the log 2 (fold change) of gene g of a specific genotype o for the T i versus T j comparison. The analysis also yields a p-value for the statistical confidence in each difference.
The GOT model applies to Experiment 2 in a straightforward way. There are 12 arrays, two treatments, and four genotypes, namely, Col-0, Cvi-0, WS, and Th.

TM4 microarray software suite
The TM4 microarray software suite consists of several components freely available from the TM4 web site [24].
In this research, we employed these components: Express-Converter, Microarray Data Analysis Software (MIDAS), and Microarray Experiment Viewer (MEV). ExpressConverter converts microarray data from various data formats, such as the GenePix Results (GPR) format, to the MEV format, which is used by MIDAS and MEV. MEV format includes only integrated intensity values (IIV), which is the kind of intensity values expected of all TM4 components.

MIDAS data normalization methods and filters
Low intensity and saturated spots are marked by quantitation programs. These spots are filtered out from the data before doing any further normalization or statistical analysis. Data normalization methods proceed from the assumption that only a relatively small proportion of the genes change significantly in expression level between the two hybridized mRNA samples. This assumption is reasonable for these data sets since the hybridizations and subsequent analysis address nearly all Arabidopsis thaliana genes. The MIDAS component of TM4 provides a number of data normalization methods and filters and supports applying them in a pipelined fashion [13,14].

Total intensity normalization
While our assumption implies that the average measured intensities of the two channels of a cDNA or oligonucleotide microarray should be almost the same, these averages are often significantly different, due to differences in the inherent fluorescence of the two dyes. The total intensity normalization step in TM4 is a straightforward means to eliminate this global dye bias. For each spot i, where 1 ≤ i ≤ n, let R i and G i be the measured intensities of the spot in the two channels. The normalized intensity data for spot i is = κG i and = R i , where κ is the normalization factor Quackenbush [14] discusses this normalization in the context of sev-eral variations that are possible to address differing channel intensities.

Lowess normalization
Beyond the global dye bias, there is dye bias that is dependent on the measured spot intensities [25,26]. TM4 constructs a scatter plot, called an RI-plot, of the points (x i , y i ), where 1 ≤ i ≤ n, given by x i = log 10 (R i G i ) and y i = log 2 (R i /G i ). Under our assumption, the RI-plot should be very nearly symmetric with respect to the line y = 0. In lowess normalization, TM4 applies the lowess method of Cleveland [27] to fit a locally weighted regression curve to the RI-plot; TM4 then adjusts spot intensities to eliminate any systematic intensity-dependent bias. Additional details on correcting intensity-dependent bias is found in [14].

Standard deviation regularization
After total intensity and lowess normalizations eliminate dye bias on a global (per microarray) scale, TM4 employs standard deviation regularization to ensure that the perblock variances of log(R i /G i ) values are the same [25,28]. Quackenbush [14] provides the formulas for this normalization step.

Low intensity filtering
Since the relative error in the log(R i /G i ) values increases if R i or G i is close to background levels, spots with low intensities are filtered out. Quackenbush [14] provides additional details, which essentially require that both R i and G i intensities be above two standard deviations of the respective backgrounds.

The MIDAS pipeline
We applied a MIDAS pipeline consisting of total intensity normalization, lowess normalization, standard deviation regularization, and low intensity filtering to both microarray data sets. MIDAS default parameters were used throughout; the default low intensity filter cut-off is R i G i < 10,000.

TM4 MEV analysis
The Multi Experiment Viewer (MEV) component of TM4 provides a number of statistical analyses and clustering algorithms to identify differentially expressed genes. We report results from the one-class t-test analysis applied to output of the MIDAS pipeline. This test assumes that the paired distribution of treated and control groups is normally distributed. Since the intensities measured from the same spot are correlated, we can apply the one-class t-test for the two-group comparison.