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 RI-plots. 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.

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. 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.

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₂(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.

Figures 3, 4, and 5 present the actual assessed Iog₂(fold change) values for 50 selected Col-0 genes in Experiment 2, along with their qRT-PCR values. These are the 50 Col-0 genes, among the 55 with qRT-PCR values, for which we have expression values for all methods. For each gene, a histogram of the Iog₂(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₂(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₂(fold change) 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.

Table 6 summarizes genotype-specific correlation results, which demonstrate that the GP model us-ing MIV input data has the highest correlation with qRT-PCR compared to the GP model and the MEV t-test using IIV_TLSF input data. The highest correlation of 0.85 is for the Col-0 results of qRT-PCR versus the GP model using MIV input. The corresponding correlations for Cvi-0, WS, and Th are 0.73, 0.66, and 0.84, respectively, which are all best among the analysis methods.

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 (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 background-corrected 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.

We used data sets from two experiments that used the Arabidopsis Oligonucleotide Microarrays 21, which include 25,712 elements, each a gene-specific 70-mer (Qiagen/Operon, Valencia, CA) for a known or putative open reading frames in Arabidopsis thaliana. There are 48 blocks per microarray, 25 rows by 24 columns (600 spots) per block, and 28,800 spots per microarray, including spots for the 25,712 gene-specific 70-mers and 302 control elements. The remaining 2,786 spots are blank.

For verification of microarray results in Experiment 2, Li, et al., 23 performed real-time quantitative reverse-transcriptase PCR (qRT-PCR) for selected genes — 55 in Col-0; 52 in Cvi-0; 59 in WS; 26 in Th. Supplementary Table 1 - see Additional file: 2 contains the annotation of the selected genes. In brief, primer pairs were selected to represent unique sequences in the Arabidopsis thaliana genome and in the Thellungiella sequences deposited in NCBI. Thellungiella actin (CX129618) cDNA primers and Arabidopsis thaliana. Ubiquitin-10 cDNA primers were used as internal controls in the qRT-PCR analyses. RT-PCR products were detected using the fluorescent dye SYBR-green (Applied Biosystems, Foster City, CA USA) and the ABI PRISM/Taqman 7900 Sequence Detection System (Applied Biosystems, Foster City, CA USA). Dissociation curves were generated for each reaction to ensure specific amplification. Three repeats were done for each gene. The averaged threshold cycle numbers were used to estimate original mRNA levels.

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-PCR. Hence, a number of the genes in Supplementary Table 1 - see Additional file: 2 are related to carbon metabolism.

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 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.

The TM4 microarray software suite consists of several components freely available from the TM4 web site 24. In this research, we employed these components: ExpressConverter, 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.

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 1314.

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.